/
2d_FVM_flow_around_cylinder.cpp
198 lines (191 loc) · 10.7 KB
/
2d_FVM_flow_around_cylinder.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
/**
* @file 2d_FVM_flow_around_cylinder.cpp
* @brief This is the test file for the weakly compressible viscous flow around a cylinder in FVM.
* @details We consider a flow passing by a cylinder in 2D in FVM framework.
* @author Zhentong Wang and Xiangyu Hu
*/
#include "common_weakly_compressible_FVM_classes.h"
#include "sphinxsys.h"
using namespace SPH;
//----------------------------------------------------------------------
// Basic geometry parameters and numerical setup.
//----------------------------------------------------------------------
Real DL = 50.0; /**< Channel length. */
Real DH = 30.0; /**< Channel height. */
Real resolution_ref = 1.0 / 5.0; /**< Initial reference particle spacing. */
Real DL_sponge = 2.0; /**< Sponge region to impose inflow condition. */
Real DH_sponge = 2.0; /**< Sponge region to impose inflow condition. */
Real cylinder_radius = 1.0; /**< Radius of the cylinder. */
//----------------------------------------------------------------------
// Material properties of the fluid.
//----------------------------------------------------------------------
Real rho0_f = 1.0; /**< Density. */
Real U_f = 1.0; /**< freestream velocity. */
Real c_f = 10.0 * U_f; /**< Speed of sound. */
Real Re = 100.0; /**< Reynolds number. */
Real mu_f = rho0_f * U_f * (2.0 * cylinder_radius) / Re; /**< Dynamics viscosity. */
//----------------------------------------------------------------------
// Set the file path to the data file.
//----------------------------------------------------------------------
std::string ansys_mesh_file_path = "./input/fluent_0.3.msh";
//----------------------------------------------------------------------
// Define geometries and body shapes
//----------------------------------------------------------------------
std::vector<Vecd> createWaterBlockShape()
{
std::vector<Vecd> water_block_shape;
water_block_shape.push_back(Vecd(-DL_sponge, -DH_sponge));
water_block_shape.push_back(Vecd(-DL_sponge, DH + DH_sponge));
water_block_shape.push_back(Vecd(DL, DH + DH_sponge));
water_block_shape.push_back(Vecd(DL, -DH_sponge));
water_block_shape.push_back(Vecd(-DL_sponge, -DH_sponge));
return water_block_shape;
}
class WaterBlock : public ComplexShape
{
public:
explicit WaterBlock(const std::string &shape_name) : ComplexShape(shape_name)
{
MultiPolygon water_block(createWaterBlockShape());
add<MultiPolygonShape>(water_block, "WaterBlock");
}
};
//----------------------------------------------------------------------
// Case dependent boundary condition
//----------------------------------------------------------------------
class FACBoundaryConditionSetup : public BoundaryConditionSetupInFVM
{
public:
FACBoundaryConditionSetup(BaseInnerRelationInFVM &inner_relation, GhostCreationFromMesh &ghost_creation)
: BoundaryConditionSetupInFVM(inner_relation, ghost_creation),
fluid_(DynamicCast<WeaklyCompressibleFluid>(this, particles_->getBaseMaterial())){};
virtual ~FACBoundaryConditionSetup(){};
void applyNonSlipWallBoundary(size_t ghost_index, size_t index_i) override
{
vel_[ghost_index] = -vel_[index_i];
p_[ghost_index] = p_[index_i];
rho_[ghost_index] = rho_[index_i];
}
void applyFarFieldBoundary(size_t ghost_index) override
{
Vecd far_field_velocity(1.0, 0.0);
Real far_field_density = 1.0;
Real far_field_pressure = fluid_.getPressure(far_field_density);
vel_[ghost_index] = far_field_velocity;
p_[ghost_index] = far_field_pressure;
rho_[ghost_index] = far_field_density;
}
protected:
Fluid &fluid_;
};
//----------------------------------------------------------------------
// Main program starts here.
//----------------------------------------------------------------------
int main(int ac, char *av[])
{
// read data from ANSYS mesh.file
ANSYSMesh ansys_mesh(ansys_mesh_file_path);
//----------------------------------------------------------------------
// Build up the environment of a SPHSystem.
//----------------------------------------------------------------------
BoundingBox system_domain_bounds(Vec2d(-DL_sponge, -DH_sponge), Vec2d(DL, DH + DH_sponge));
SPHSystem sph_system(system_domain_bounds, ansys_mesh.MinMeshEdge());
sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment();
//----------------------------------------------------------------------
// Creating body, materials and particles.
//----------------------------------------------------------------------
FluidBody water_block(sph_system, makeShared<WaterBlock>("WaterBlock"));
water_block.defineParticlesAndMaterial<BaseParticles, WeaklyCompressibleFluid>(rho0_f, c_f, mu_f);
Ghost<ReserveSizeFactor> ghost_boundary(0.5);
water_block.generateParticlesWithReserve<UnstructuredMesh>(ghost_boundary, ansys_mesh);
water_block.addBodyStateForRecording<Real>("Density");
GhostCreationFromMesh ghost_creation(water_block, ansys_mesh, ghost_boundary);
//----------------------------------------------------------------------
// Define body relation map.
//----------------------------------------------------------------------
InnerRelationInFVM water_block_inner(water_block, ansys_mesh);
//----------------------------------------------------------------------
// Define the main numerical methods used in the simulation.
// Note that there may be data dependence on the constructors of these methods.
//----------------------------------------------------------------------
/** Here we introduce the limiter in the Riemann solver and 0 means the no extra numerical dissipation.
the value is larger, the numerical dissipation larger*/
InteractionWithUpdate<fluid_dynamics::EulerianIntegration1stHalfInnerRiemann> pressure_relaxation(water_block_inner, 200.0);
InteractionWithUpdate<fluid_dynamics::EulerianIntegration2ndHalfInnerRiemann> density_relaxation(water_block_inner, 200.0);
FACBoundaryConditionSetup boundary_condition_setup(water_block_inner, ghost_creation);
ReduceDynamics<fluid_dynamics::WCAcousticTimeStepSizeInFVM> get_fluid_time_step_size(water_block, ansys_mesh.MinMeshEdge());
InteractionWithUpdate<fluid_dynamics::ViscousForceInner> viscous_force(water_block_inner);
//----------------------------------------------------------------------
// Compute the force exerted on solid body due to fluid pressure and viscosity
//----------------------------------------------------------------------
InteractionDynamics<fluid_dynamics::ViscousForceFromFluidInFVM> viscous_force_on_solid(water_block_inner, ghost_creation.each_boundary_type_contact_real_index_);
InteractionDynamics<fluid_dynamics::PressureForceFromFluidInFVM<decltype(density_relaxation)>> pressure_force_on_solid(water_block_inner, ghost_creation.each_boundary_type_contact_real_index_);
//----------------------------------------------------------------------
// Define the methods for I/O operations and observations of the simulation.
//----------------------------------------------------------------------
BodyStatesRecordingInMeshToVtp write_real_body_states(water_block, ansys_mesh);
RegressionTestDynamicTimeWarping<ReducedQuantityRecording<QuantitySummation<Vecd>>> write_total_viscous_force_on_inserted_body(water_block, "ViscousForceOnSolid");
ReducedQuantityRecording<QuantitySummation<Vecd>> write_total_pressure_force_on_inserted_body(water_block, "PressureForceOnSolid");
ReducedQuantityRecording<MaximumSpeed> write_maximum_speed(water_block);
//----------------------------------------------------------------------
// Prepare the simulation with cell linked list, configuration
// and case specified initial condition if necessary.
//----------------------------------------------------------------------
water_block_inner.updateConfiguration();
//----------------------------------------------------------------------
// Setup for time-stepping control
//----------------------------------------------------------------------
size_t number_of_iterations = 0;
int screen_output_interval = 1000;
Real end_time = 70.0;
Real output_interval = 5.0; /**< time stamps for output. */
//----------------------------------------------------------------------
// Statistics for CPU time
//----------------------------------------------------------------------
TickCount t1 = TickCount::now();
TimeInterval interval;
//----------------------------------------------------------------------
// First output before the main loop.
//----------------------------------------------------------------------
write_real_body_states.writeToFile(0);
//----------------------------------------------------------------------
// Main loop starts here.
//----------------------------------------------------------------------
while (GlobalStaticVariables::physical_time_ < end_time)
{
Real integration_time = 0.0;
while (integration_time < output_interval)
{
Real dt = get_fluid_time_step_size.exec();
boundary_condition_setup.resetBoundaryConditions();
viscous_force.exec();
pressure_relaxation.exec(dt);
boundary_condition_setup.resetBoundaryConditions();
density_relaxation.exec(dt);
integration_time += dt;
GlobalStaticVariables::physical_time_ += dt;
if (number_of_iterations % screen_output_interval == 0)
{
cout << fixed << setprecision(9) << "N=" << number_of_iterations << " Time = "
<< GlobalStaticVariables::physical_time_
<< " dt = " << dt << "\n";
viscous_force_on_solid.exec();
pressure_force_on_solid.exec();
write_total_viscous_force_on_inserted_body.writeToFile(number_of_iterations);
write_total_pressure_force_on_inserted_body.writeToFile(number_of_iterations);
write_maximum_speed.writeToFile(number_of_iterations);
}
number_of_iterations++;
}
TickCount t2 = TickCount::now();
write_real_body_states.writeToFile();
TickCount t3 = TickCount::now();
interval += t3 - t2;
}
TickCount t4 = TickCount::now();
TimeInterval tt;
tt = t4 - t1 - interval;
cout << "Total wall time for computation: " << tt.seconds() << " seconds." << endl;
write_total_viscous_force_on_inserted_body.testResult();
return 0;
}