Skip to content

Commit 464a9d9

Browse files
Marco1410Marco1410
authored
Merge pull request #14272 from KratosMultiphysics/marco/convection_diffusion
[ConvectionDiffusionApplication] Adding crosswind stabilization to Eulerian Convection–Diffusion element
2 parents 30603e0 + 6418514 commit 464a9d9

11 files changed

Lines changed: 5385 additions & 2 deletions

‎applications/ConvectionDiffusionApplication/custom_elements/eulerian_conv_diff.cpp‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -103,6 +103,9 @@ namespace Kratos
103103
this-> GetNodalValues(Variables,rCurrentProcessInfo);
104104
double h = this->ComputeH(DN_DX);
105105

106+
array_1d<double,TDim> grad_phi_halfstep = prod(trans(DN_DX), 0.5*(Variables.phi+Variables.phi_old));
107+
const double norm_grad = norm_2(grad_phi_halfstep);
108+
106109
//Computing the divergence
107110
for (unsigned int i = 0; i < TNumNodes; i++)
108111
{
@@ -116,6 +119,11 @@ namespace Kratos
116119
BoundedMatrix<double,TNumNodes, TNumNodes> aux1 = ZeroMatrix(TNumNodes, TNumNodes); //terms multiplying dphi/dt
117120
BoundedMatrix<double,TNumNodes, TNumNodes> aux2 = ZeroMatrix(TNumNodes, TNumNodes); //terms multiplying phi
118121
bounded_matrix<double,TNumNodes, TDim> tmp;
122+
// CrossWind variables
123+
// Projected velocity
124+
array_1d<double,TDim> u_proj;
125+
// Diffusion contribution
126+
BoundedMatrix<double,TDim,TDim> Dcw;
119127

120128
// Gauss points and Number of nodes coincides in this case.
121129
for(unsigned int igauss=0; igauss<TNumNodes; igauss++)
@@ -141,6 +149,65 @@ namespace Kratos
141149
//terms which multiply the gradient of phi
142150
noalias(aux2) += (1.0+tau*Variables.beta*Variables.div_v)*outer_prod(N, a_dot_grad);
143151
noalias(aux2) += tau*outer_prod(a_dot_grad, a_dot_grad);
152+
153+
// Crosswind term according to https://doi.org/10.1016/0045-7825(93)90213-H
154+
const double norm_vel2 = norm_vel * norm_vel;
155+
if(Variables.crosswind_constant > 0.0 && norm_grad > 1e-3 && norm_vel2 > 1e-9)
156+
{
157+
// Temporal derivative
158+
const double phi_gauss = inner_prod(N, Variables.phi);
159+
const double phi_old_gauss = inner_prod(N, Variables.phi_old);
160+
const double dphi_dt = Variables.dt_inv * (phi_gauss - phi_old_gauss);
161+
162+
// Convective term
163+
const double convection = inner_prod(vel_gauss, grad_phi_halfstep);
164+
165+
// Reaction term
166+
const double reaction = Variables.beta * Variables.div_v * phi_gauss;
167+
168+
// Source termn
169+
const double source = inner_prod(N, Variables.volumetric_source);
170+
171+
// Complete residual
172+
const double residual = dphi_dt + convection + reaction - source;
173+
174+
// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
175+
// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
176+
// // Peclet's projected formulation was commented after some testing due to
177+
// // residual radial instabilities in 3D cases.
178+
// // It was decided to use a complete approach with 'crosswind_constant' without
179+
// // considering the projection.
180+
// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
181+
// // Velocity projected in the direction of the solution gradient
182+
// u_proj = (convection/std::pow(norm_grad,2)) * grad_phi_halfstep;
183+
184+
// // Peclet number projected in the direction of the solution gradient
185+
// const double Pe_proj = norm_2(u_proj) * h / (2.0 * Variables.conductivity + 1e-12);
186+
187+
// // Limiter
188+
// const double alpha_c = std::max( 0.0, Variables.crosswind_constant - 1.0/(Pe_proj + 1e-12));
189+
190+
// // Discontinuity capturing coefficient
191+
// double k_c = 0.5 * alpha_c * h * std::abs(residual / norm_grad);
192+
// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
193+
// //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
194+
195+
// Discontinuity capturing coefficient
196+
double k_c = Variables.crosswind_constant * h * std::abs(residual / norm_grad);
197+
198+
// Crosswind diffusion tensor
199+
Dcw = k_c * IdentityMatrix(TDim);
200+
201+
// Removing diffusion in streamline direcction
202+
const double k_streamline = tau * norm_vel2;
203+
const double correction = std::max(k_c - k_streamline, 0.0) - k_c;
204+
205+
Dcw += (correction / norm_vel2) * outer_prod(vel_gauss, vel_gauss);
206+
207+
// Add diffusion contribution
208+
noalias(tmp) = prod(DN_DX, Dcw);
209+
noalias(aux2) += prod(tmp, trans(DN_DX));
210+
}
144211
}
145212

146213
//adding the second and third term in the formulation
@@ -174,6 +241,7 @@ namespace Kratos
174241
{
175242
KRATOS_TRY
176243

244+
rVariables.crosswind_constant = rCurrentProcessInfo[CROSS_WIND_STABILIZATION_FACTOR];
177245
rVariables.theta = rCurrentProcessInfo[TIME_INTEGRATION_THETA]; //Variable defining the temporal scheme (0: Forward Euler, 1: Backward Euler, 0.5: Crank-Nicolson)
178246
rVariables.dyn_st_beta = rCurrentProcessInfo[DYNAMIC_TAU];
179247
const double delta_t = rCurrentProcessInfo[DELTA_TIME];

‎applications/ConvectionDiffusionApplication/custom_elements/eulerian_conv_diff.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -119,6 +119,7 @@ class KRATOS_API(CONVECTION_DIFFUSION_APPLICATION) EulerianConvectionDiffusionEl
119119
double density;
120120
double beta;
121121
double div_v;
122+
double crosswind_constant;
122123

123124
array_1d<double,TNumNodes> phi;
124125
array_1d<double,TNumNodes> phi_old;

‎applications/ConvectionDiffusionApplication/python_scripts/convection_diffusion_transient_solver.py‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,8 @@ def GetDefaultParameters(cls):
3939
"time_integration_method" : "implicit",
4040
"transient_parameters" : {
4141
"dynamic_tau": 1.0,
42-
"theta" : 0.5
42+
"theta" : 0.5,
43+
"cross_wind_stabilization_factor" : 0.0
4344
}
4445
}""")
4546
this_defaults.AddMissingParameters(super().GetDefaultParameters())
@@ -50,6 +51,7 @@ def _CreateScheme(self):
5051
# Variable defining the temporal scheme (0: Forward Euler, 1: Backward Euler, 0.5: Crank-Nicolson)
5152
self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.TIME_INTEGRATION_THETA] = self.settings["transient_parameters"]["theta"].GetDouble()
5253
self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DYNAMIC_TAU] = self.settings["transient_parameters"]["dynamic_tau"].GetDouble()
54+
self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.CROSS_WIND_STABILIZATION_FACTOR] = self.settings["transient_parameters"]["cross_wind_stabilization_factor"].GetDouble()
5355

5456
# As the time integration is managed by the element, we set a "fake" scheme to perform the solution update
5557
if not self.main_model_part.IsDistributed():

‎applications/ConvectionDiffusionApplication/tests/SourceTermTest/SourceTermTestProjectParameters.json‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,8 @@
4545
},
4646
"transient_parameters" : {
4747
"dynamic_tau": 1.0,
48-
"theta" : 1.0
48+
"theta" : 1.0,
49+
"cross_wind_stabilization_factor": 0.0
4950
}
5051
},
5152
"processes" : {

‎applications/ConvectionDiffusionApplication/tests/cpp_tests/test_eulerian_conv_diff.cpp‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -256,4 +256,60 @@ namespace Kratos::Testing
256256
KRATOS_EXPECT_MATRIX_NEAR(LHS, expected_LHS, 1.0e-4)
257257
}
258258

259+
KRATOS_TEST_CASE_IN_SUITE(EulerianConvDiff2D3NCrossWindStabilization, KratosConvectionDiffusionFastSuite)
260+
{
261+
// Create the test element
262+
Model model;
263+
auto &r_test_model_part = model.CreateModelPart("TestModelPart");
264+
ConvectionDiffusionTestingUtilities::SetEntityUnitTestModelPart(r_test_model_part);
265+
266+
// Fill the process info container
267+
auto &r_process_info = r_test_model_part.GetProcessInfo();
268+
r_process_info.SetValue(CROSS_WIND_STABILIZATION_FACTOR, 0.7);
269+
r_process_info.SetValue(TIME_INTEGRATION_THETA, 1.0);
270+
r_process_info.SetValue(DELTA_TIME, 0.1);
271+
r_process_info.SetValue(DYNAMIC_TAU, 0.001);
272+
273+
// Element creation
274+
r_test_model_part.CreateNewNode(1, 0.0, 0.0, 0.0);
275+
r_test_model_part.CreateNewNode(2, 1.0, 0.0, 0.0);
276+
r_test_model_part.CreateNewNode(3, 0.0, 1.0, 0.0);
277+
std::vector<ModelPart::IndexType> elem_nodes{1, 2, 3};
278+
r_test_model_part.CreateNewElement("EulerianConvDiff2D", 1, elem_nodes, r_test_model_part.pGetProperties(0));
279+
280+
// Set the nodal values
281+
for (auto &i_node : r_test_model_part.Nodes()) {
282+
i_node.FastGetSolutionStepValue(DENSITY) = 1.0;
283+
i_node.FastGetSolutionStepValue(CONDUCTIVITY) = 1.0;
284+
i_node.FastGetSolutionStepValue(SPECIFIC_HEAT) = 1.0;
285+
array_1d<double,3> aux_vel = ZeroVector(3);
286+
aux_vel[0] = i_node.X();
287+
aux_vel[1] = i_node.Y();
288+
i_node.FastGetSolutionStepValue(VELOCITY) = aux_vel;
289+
}
290+
291+
// Test element
292+
auto p_element = r_test_model_part.pGetElement(1);
293+
Vector RHS = ZeroVector(3);
294+
Matrix LHS = ZeroMatrix(3,3);
295+
const auto& rConstProcessInfoRef = r_test_model_part.GetProcessInfo();
296+
p_element->CalculateLocalSystem(LHS, RHS, rConstProcessInfoRef);
297+
298+
std::vector<double> expected_RHS = {0.0, 0.0, 0.0};
299+
Matrix expected_LHS(3, 3);
300+
expected_LHS(0, 0) = 1.71340;
301+
expected_LHS(0, 1) = -0.12313;
302+
expected_LHS(0, 2) = -0.12313;
303+
expected_LHS(1, 0) = -0.19003;
304+
expected_LHS(1, 1) = 1.47086;
305+
expected_LHS(1, 2) = 0.48560;
306+
expected_LHS(2, 0) = -0.19003;
307+
expected_LHS(2, 1) = 0.48560;
308+
expected_LHS(2, 2) = 1.47086;
309+
310+
KRATOS_EXPECT_VECTOR_NEAR(RHS, expected_RHS, 1.0e-4)
311+
KRATOS_EXPECT_MATRIX_NEAR(LHS, expected_LHS, 1.0e-4)
312+
313+
}
314+
259315
} // namespace Kratos::Testing.

‎applications/ConvectionDiffusionApplication/tests/discontinuity_propagation_unit_square_test/ConvectionDiffusionMaterials.json‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
{
2+
"properties" : [{
3+
"model_part_name" : "ThermalModelPart",
4+
"properties_id" : 0,
5+
"Material" : {
6+
"Variables" : []
7+
}
8+
},{
9+
"model_part_name" : "ThermalModelPart.Body",
10+
"properties_id" : 1,
11+
"Material" : {
12+
"Variables" : {
13+
"DENSITY" : 1.0,
14+
"CONDUCTIVITY" : 1e-8,
15+
"SPECIFIC_HEAT" : 1.0
16+
},
17+
"Tables" : null
18+
}
19+
}]
20+
}

‎applications/ConvectionDiffusionApplication/tests/discontinuity_propagation_unit_square_test/conv_diffusion_test_parameters.json‎

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
{
2+
"analysis_stage": "KratosMultiphysics.ConvectionDiffusionApplication.convection_diffusion_analysis",
3+
4+
"problem_data": {
5+
"problem_name": "codina_test1",
6+
"parallel_type": "OpenMP",
7+
"start_time": 0.0,
8+
"end_time": 0.5,
9+
"echo_level": 1
10+
},
11+
"modelers" : [{
12+
"name" : "Modelers.KratosMultiphysics.ImportMDPAModeler",
13+
"parameters" : {
14+
"input_filename" : "square20x20",
15+
"model_part_name" : "ThermalModelPart"
16+
}
17+
},{
18+
"name" : "Modelers.KratosMultiphysics.CreateEntitiesFromGeometriesModeler",
19+
"parameters" : {
20+
"elements_list" : [{
21+
"model_part_name" : "ThermalModelPart.Body",
22+
"element_name" : "Element2D3N"
23+
}],
24+
"conditions_list" : [{
25+
"model_part_name" : "ThermalModelPart.Top_Wall",
26+
"condition_name" : "LineCondition2D2N"
27+
},{
28+
"model_part_name" : "ThermalModelPart.Bottom_Wall",
29+
"condition_name" : "LineCondition2D2N"
30+
},{
31+
"model_part_name" : "ThermalModelPart.Left_Wall",
32+
"condition_name" : "LineCondition2D2N"
33+
},{
34+
"model_part_name" : "ThermalModelPart.Right_Wall",
35+
"condition_name" : "LineCondition2D2N"
36+
}]
37+
}
38+
}],
39+
"solver_settings": {
40+
"solver_type": "transient",
41+
"analysis_type": "linear",
42+
"model_part_name": "ThermalModelPart",
43+
"domain_size": 2,
44+
"model_import_settings": {
45+
"input_type": "use_input_model_part"
46+
},
47+
"material_import_settings": {
48+
"materials_filename": "ConvectionDiffusionMaterials.json"
49+
},
50+
"time_stepping": {
51+
"time_step": 0.01 // reduce dt for sharper resolution (0.0001)
52+
},
53+
"transient_parameters" : {
54+
"dynamic_tau": 1.0,
55+
"theta" : 0.5,
56+
// "cross_wind_stabilization_factor" : 0.0
57+
"cross_wind_stabilization_factor" : 0.7
58+
},
59+
"convection_diffusion_variables": {
60+
"unknown_variable": "TEMPERATURE",
61+
"velocity_variable": "VELOCITY",
62+
"diffusion_variable": "CONDUCTIVITY",
63+
"specific_heat_variable": "SPECIFIC_HEAT",
64+
"density_variable": "DENSITY",
65+
"volume_source_variable": "HEAT_FLUX"
66+
}
67+
},
68+
"processes": {
69+
"initial_conditions_process_list": [
70+
{
71+
"python_module": "assign_scalar_variable_process",
72+
"kratos_module": "KratosMultiphysics",
73+
"process_name": "AssignScalarVariableProcess",
74+
"Parameters": {
75+
"model_part_name": "ThermalModelPart.Body",
76+
"variable_name": "TEMPERATURE",
77+
"interval": [0.0, 0.0],
78+
"value": 0.0
79+
}
80+
}
81+
],
82+
"boundary_conditions_process_list": [],
83+
"auxiliar_process_list": [
84+
{
85+
"python_module": "assign_vector_variable_process",
86+
"kratos_module": "KratosMultiphysics",
87+
"process_name": "AssignVectorVariableProcess",
88+
"Parameters": {
89+
"model_part_name": "ThermalModelPart.Body",
90+
"variable_name": "VELOCITY",
91+
"interval": [0.0, "End"],
92+
"constrained": false,
93+
"value": [1.0, -2.0, 0.0]
94+
}
95+
}
96+
]
97+
},
98+
"output_processes": {
99+
// "gid_output" : [{
100+
// "python_module" : "gid_output_process",
101+
// "kratos_module" : "KratosMultiphysics",
102+
// "process_name" : "GiDOutputProcess",
103+
// "Parameters" : {
104+
// "model_part_name" : "ThermalModelPart",
105+
// "postprocess_parameters" : {
106+
// "result_file_configuration" : {
107+
// "gidpost_flags" : {
108+
// "GiDPostMode" : "GiD_PostBinary",
109+
// "WriteDeformedMeshFlag" : "WriteDeformed",
110+
// "WriteConditionsFlag" : "WriteConditions",
111+
// "MultiFileFlag" : "SingleFile"
112+
// },
113+
// "file_label" : "time",
114+
// "output_control_type" : "step",
115+
// "output_interval" : 1,
116+
// "body_output" : true,
117+
// "node_output" : false,
118+
// "skin_output" : false,
119+
// "plane_output" : [],
120+
// "nodal_results" : ["TEMPERATURE","VELOCITY"],
121+
// "nodal_nonhistorical_results" : [],
122+
// "gauss_point_results" : []
123+
// },
124+
// "point_data_configuration" : []
125+
// },
126+
// "output_name" : "gid_output/crosswind"
127+
// }
128+
// }]
129+
// ,
130+
// "vtk_output" : [{
131+
// "python_module" : "vtk_output_process",
132+
// "kratos_module" : "KratosMultiphysics",
133+
// "process_name" : "VtkOutputProcess",
134+
// "Parameters" : {
135+
// "model_part_name" : "ThermalModelPart",
136+
// "output_control_type" : "step",
137+
// "output_interval" : 100,
138+
// "file_format" : "ascii",
139+
// "output_precision" : 7,
140+
// "output_sub_model_parts" : false,
141+
// "output_path" : "crosswind",
142+
// "save_output_files_in_folder" : true,
143+
// "nodal_solution_step_data_variables" : ["TEMPERATURE"],
144+
// "nodal_data_value_variables" : [],
145+
// "element_data_value_variables" : [],
146+
// "condition_data_value_variables" : [],
147+
// "gauss_point_variables_extrapolated_to_nodes" : []
148+
// }
149+
// }]
150+
}
151+
}

0 commit comments

Comments
 (0)