Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,9 @@ namespace Kratos
this-> GetNodalValues(Variables,rCurrentProcessInfo);
double h = this->ComputeH(DN_DX);

array_1d<double,TDim> grad_phi_halfstep = prod(trans(DN_DX), 0.5*(Variables.phi+Variables.phi_old));
const double norm_grad = norm_2(grad_phi_halfstep);

//Computing the divergence
for (unsigned int i = 0; i < TNumNodes; i++)
{
Expand All @@ -116,6 +119,11 @@ namespace Kratos
BoundedMatrix<double,TNumNodes, TNumNodes> aux1 = ZeroMatrix(TNumNodes, TNumNodes); //terms multiplying dphi/dt
BoundedMatrix<double,TNumNodes, TNumNodes> aux2 = ZeroMatrix(TNumNodes, TNumNodes); //terms multiplying phi
bounded_matrix<double,TNumNodes, TDim> tmp;
// CrossWind variables
// Projected velocity
array_1d<double,TDim> u_proj;
// Diffusion contribution
BoundedMatrix<double,TDim,TDim> Dcw;

// Gauss points and Number of nodes coincides in this case.
for(unsigned int igauss=0; igauss<TNumNodes; igauss++)
Expand All @@ -141,6 +149,65 @@ namespace Kratos
//terms which multiply the gradient of phi
noalias(aux2) += (1.0+tau*Variables.beta*Variables.div_v)*outer_prod(N, a_dot_grad);
noalias(aux2) += tau*outer_prod(a_dot_grad, a_dot_grad);

// Crosswind term according to https://doi.org/10.1016/0045-7825(93)90213-H
const double norm_vel2 = norm_vel * norm_vel;
if(Variables.crosswind_constant > 0.0 && norm_grad > 1e-3 && norm_vel2 > 1e-9)
{
// Temporal derivative
const double phi_gauss = inner_prod(N, Variables.phi);
const double phi_old_gauss = inner_prod(N, Variables.phi_old);
const double dphi_dt = Variables.dt_inv * (phi_gauss - phi_old_gauss);

// Convective term
const double convection = inner_prod(vel_gauss, grad_phi_halfstep);

// Reaction term
const double reaction = Variables.beta * Variables.div_v * phi_gauss;

// Source termn
const double source = inner_prod(N, Variables.volumetric_source);

// Complete residual
const double residual = dphi_dt + convection + reaction - source;

// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// // Peclet's projected formulation was commented after some testing due to
// // residual radial instabilities in 3D cases.
// // It was decided to use a complete approach with 'crosswind_constant' without
// // considering the projection.
// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// // Velocity projected in the direction of the solution gradient
// u_proj = (convection/std::pow(norm_grad,2)) * grad_phi_halfstep;

// // Peclet number projected in the direction of the solution gradient
// const double Pe_proj = norm_2(u_proj) * h / (2.0 * Variables.conductivity + 1e-12);

// // Limiter
// const double alpha_c = std::max( 0.0, Variables.crosswind_constant - 1.0/(Pe_proj + 1e-12));

// // Discontinuity capturing coefficient
// double k_c = 0.5 * alpha_c * h * std::abs(residual / norm_grad);
// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

// Discontinuity capturing coefficient
double k_c = Variables.crosswind_constant * h * std::abs(residual / norm_grad);

// Crosswind diffusion tensor
Dcw = k_c * IdentityMatrix(TDim);

// Removing diffusion in streamline direcction
const double k_streamline = tau * norm_vel2;
const double correction = std::max(k_c - k_streamline, 0.0) - k_c;

Dcw += (correction / norm_vel2) * outer_prod(vel_gauss, vel_gauss);

// Add diffusion contribution
noalias(tmp) = prod(DN_DX, Dcw);
noalias(aux2) += prod(tmp, trans(DN_DX));
}
}

//adding the second and third term in the formulation
Expand Down Expand Up @@ -174,6 +241,7 @@ namespace Kratos
{
KRATOS_TRY

rVariables.crosswind_constant = rCurrentProcessInfo[CROSS_WIND_STABILIZATION_FACTOR];
rVariables.theta = rCurrentProcessInfo[TIME_INTEGRATION_THETA]; //Variable defining the temporal scheme (0: Forward Euler, 1: Backward Euler, 0.5: Crank-Nicolson)
rVariables.dyn_st_beta = rCurrentProcessInfo[DYNAMIC_TAU];
const double delta_t = rCurrentProcessInfo[DELTA_TIME];
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,7 @@ class KRATOS_API(CONVECTION_DIFFUSION_APPLICATION) EulerianConvectionDiffusionEl
double density;
double beta;
double div_v;
double crosswind_constant;

array_1d<double,TNumNodes> phi;
array_1d<double,TNumNodes> phi_old;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,8 @@ def GetDefaultParameters(cls):
"time_integration_method" : "implicit",
"transient_parameters" : {
"dynamic_tau": 1.0,
"theta" : 0.5
"theta" : 0.5,
"cross_wind_stabilization_factor" : 0.0
}
}""")
this_defaults.AddMissingParameters(super().GetDefaultParameters())
Expand All @@ -50,6 +51,7 @@ def _CreateScheme(self):
# Variable defining the temporal scheme (0: Forward Euler, 1: Backward Euler, 0.5: Crank-Nicolson)
self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.TIME_INTEGRATION_THETA] = self.settings["transient_parameters"]["theta"].GetDouble()
self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DYNAMIC_TAU] = self.settings["transient_parameters"]["dynamic_tau"].GetDouble()
self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.CROSS_WIND_STABILIZATION_FACTOR] = self.settings["transient_parameters"]["cross_wind_stabilization_factor"].GetDouble()

# As the time integration is managed by the element, we set a "fake" scheme to perform the solution update
if not self.main_model_part.IsDistributed():
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,8 @@
},
"transient_parameters" : {
"dynamic_tau": 1.0,
"theta" : 1.0
"theta" : 1.0,
"cross_wind_stabilization_factor": 0.0
Comment thread
Marco1410 marked this conversation as resolved.
}
},
"processes" : {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -256,4 +256,60 @@ namespace Kratos::Testing
KRATOS_EXPECT_MATRIX_NEAR(LHS, expected_LHS, 1.0e-4)
}

KRATOS_TEST_CASE_IN_SUITE(EulerianConvDiff2D3NCrossWindStabilization, KratosConvectionDiffusionFastSuite)
{
// Create the test element
Model model;
auto &r_test_model_part = model.CreateModelPart("TestModelPart");
ConvectionDiffusionTestingUtilities::SetEntityUnitTestModelPart(r_test_model_part);

// Fill the process info container
auto &r_process_info = r_test_model_part.GetProcessInfo();
r_process_info.SetValue(CROSS_WIND_STABILIZATION_FACTOR, 0.7);
r_process_info.SetValue(TIME_INTEGRATION_THETA, 1.0);
r_process_info.SetValue(DELTA_TIME, 0.1);
r_process_info.SetValue(DYNAMIC_TAU, 0.001);

// Element creation
r_test_model_part.CreateNewNode(1, 0.0, 0.0, 0.0);
r_test_model_part.CreateNewNode(2, 1.0, 0.0, 0.0);
r_test_model_part.CreateNewNode(3, 0.0, 1.0, 0.0);
std::vector<ModelPart::IndexType> elem_nodes{1, 2, 3};
r_test_model_part.CreateNewElement("EulerianConvDiff2D", 1, elem_nodes, r_test_model_part.pGetProperties(0));

// Set the nodal values
for (auto &i_node : r_test_model_part.Nodes()) {
i_node.FastGetSolutionStepValue(DENSITY) = 1.0;
i_node.FastGetSolutionStepValue(CONDUCTIVITY) = 1.0;
i_node.FastGetSolutionStepValue(SPECIFIC_HEAT) = 1.0;
array_1d<double,3> aux_vel = ZeroVector(3);
aux_vel[0] = i_node.X();
aux_vel[1] = i_node.Y();
i_node.FastGetSolutionStepValue(VELOCITY) = aux_vel;
}

// Test element
auto p_element = r_test_model_part.pGetElement(1);
Vector RHS = ZeroVector(3);
Matrix LHS = ZeroMatrix(3,3);
const auto& rConstProcessInfoRef = r_test_model_part.GetProcessInfo();
p_element->CalculateLocalSystem(LHS, RHS, rConstProcessInfoRef);

std::vector<double> expected_RHS = {0.0, 0.0, 0.0};
Matrix expected_LHS(3, 3);
expected_LHS(0, 0) = 1.71340;
expected_LHS(0, 1) = -0.12313;
expected_LHS(0, 2) = -0.12313;
expected_LHS(1, 0) = -0.19003;
expected_LHS(1, 1) = 1.47086;
expected_LHS(1, 2) = 0.48560;
expected_LHS(2, 0) = -0.19003;
expected_LHS(2, 1) = 0.48560;
expected_LHS(2, 2) = 1.47086;

KRATOS_EXPECT_VECTOR_NEAR(RHS, expected_RHS, 1.0e-4)
KRATOS_EXPECT_MATRIX_NEAR(LHS, expected_LHS, 1.0e-4)
Comment thread
Marco1410 marked this conversation as resolved.

}

} // namespace Kratos::Testing.
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
{
"properties" : [{
"model_part_name" : "ThermalModelPart",
"properties_id" : 0,
"Material" : {
"Variables" : []
}
},{
"model_part_name" : "ThermalModelPart.Body",
"properties_id" : 1,
"Material" : {
"Variables" : {
"DENSITY" : 1.0,
"CONDUCTIVITY" : 1e-8,
"SPECIFIC_HEAT" : 1.0
},
"Tables" : null
}
}]
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,151 @@
{
"analysis_stage": "KratosMultiphysics.ConvectionDiffusionApplication.convection_diffusion_analysis",

"problem_data": {
"problem_name": "codina_test1",
"parallel_type": "OpenMP",
"start_time": 0.0,
"end_time": 0.5,
"echo_level": 1
},
"modelers" : [{
"name" : "Modelers.KratosMultiphysics.ImportMDPAModeler",
"parameters" : {
"input_filename" : "square20x20",
"model_part_name" : "ThermalModelPart"
}
},{
"name" : "Modelers.KratosMultiphysics.CreateEntitiesFromGeometriesModeler",
"parameters" : {
"elements_list" : [{
"model_part_name" : "ThermalModelPart.Body",
"element_name" : "Element2D3N"
}],
"conditions_list" : [{
"model_part_name" : "ThermalModelPart.Top_Wall",
"condition_name" : "LineCondition2D2N"
},{
"model_part_name" : "ThermalModelPart.Bottom_Wall",
"condition_name" : "LineCondition2D2N"
},{
"model_part_name" : "ThermalModelPart.Left_Wall",
"condition_name" : "LineCondition2D2N"
},{
"model_part_name" : "ThermalModelPart.Right_Wall",
"condition_name" : "LineCondition2D2N"
}]
}
}],
"solver_settings": {
"solver_type": "transient",
"analysis_type": "linear",
"model_part_name": "ThermalModelPart",
"domain_size": 2,
"model_import_settings": {
"input_type": "use_input_model_part"
},
"material_import_settings": {
"materials_filename": "ConvectionDiffusionMaterials.json"
},
"time_stepping": {
"time_step": 0.01 // reduce dt for sharper resolution (0.0001)
},
"transient_parameters" : {
"dynamic_tau": 1.0,
"theta" : 0.5,
// "cross_wind_stabilization_factor" : 0.0
"cross_wind_stabilization_factor" : 0.7
},
"convection_diffusion_variables": {
"unknown_variable": "TEMPERATURE",
"velocity_variable": "VELOCITY",
"diffusion_variable": "CONDUCTIVITY",
"specific_heat_variable": "SPECIFIC_HEAT",
"density_variable": "DENSITY",
"volume_source_variable": "HEAT_FLUX"
}
},
"processes": {
"initial_conditions_process_list": [
{
"python_module": "assign_scalar_variable_process",
"kratos_module": "KratosMultiphysics",
"process_name": "AssignScalarVariableProcess",
"Parameters": {
"model_part_name": "ThermalModelPart.Body",
"variable_name": "TEMPERATURE",
"interval": [0.0, 0.0],
"value": 0.0
}
}
],
"boundary_conditions_process_list": [],
"auxiliar_process_list": [
{
"python_module": "assign_vector_variable_process",
"kratos_module": "KratosMultiphysics",
"process_name": "AssignVectorVariableProcess",
"Parameters": {
"model_part_name": "ThermalModelPart.Body",
"variable_name": "VELOCITY",
"interval": [0.0, "End"],
"constrained": false,
"value": [1.0, -2.0, 0.0]
}
}
]
},
"output_processes": {
// "gid_output" : [{
// "python_module" : "gid_output_process",
// "kratos_module" : "KratosMultiphysics",
// "process_name" : "GiDOutputProcess",
// "Parameters" : {
// "model_part_name" : "ThermalModelPart",
// "postprocess_parameters" : {
// "result_file_configuration" : {
// "gidpost_flags" : {
// "GiDPostMode" : "GiD_PostBinary",
// "WriteDeformedMeshFlag" : "WriteDeformed",
// "WriteConditionsFlag" : "WriteConditions",
// "MultiFileFlag" : "SingleFile"
// },
// "file_label" : "time",
// "output_control_type" : "step",
// "output_interval" : 1,
// "body_output" : true,
// "node_output" : false,
// "skin_output" : false,
// "plane_output" : [],
// "nodal_results" : ["TEMPERATURE","VELOCITY"],
// "nodal_nonhistorical_results" : [],
// "gauss_point_results" : []
// },
// "point_data_configuration" : []
// },
// "output_name" : "gid_output/crosswind"
// }
// }]
// ,
// "vtk_output" : [{
// "python_module" : "vtk_output_process",
// "kratos_module" : "KratosMultiphysics",
// "process_name" : "VtkOutputProcess",
// "Parameters" : {
// "model_part_name" : "ThermalModelPart",
// "output_control_type" : "step",
// "output_interval" : 100,
// "file_format" : "ascii",
// "output_precision" : 7,
// "output_sub_model_parts" : false,
// "output_path" : "crosswind",
// "save_output_files_in_folder" : true,
// "nodal_solution_step_data_variables" : ["TEMPERATURE"],
// "nodal_data_value_variables" : [],
// "element_data_value_variables" : [],
// "condition_data_value_variables" : [],
// "gauss_point_variables_extrapolated_to_nodes" : []
// }
// }]
}
}
Loading
Loading