第1页 / 共28页
第2页 / 共28页
第3页 / 共28页
第4页 / 共28页
第5页 / 共28页
第6页 / 共28页
第7页 / 共28页
第8页 / 共28页
immersed boundary method.docx
2
4 5
6 7
3
0
0 1
0 +1 -1
0 0
1. // This is an example 2D immersed boundary lattice Boltzmann method code.
2. // It uses the D2Q9 lattice with Guo's forcing term.
3. // Rigid bottom and top walls are parallel to the x-axis (channel).
4. // The flow is periodic along the x-axis.
5. // One initially cylindrical particle is positioned in the flow.
6. // This particle can be rigid/deformable and stationary/moving.
7. //
8. // Last update 21-Aug-2011 by Timm Krüger.
9. // This code may be changed and distributed freely.
10. //
11. // The lattice velocities are defined according to the following scheme:
12. // index:
8
13. // ----------------------------------
14. // x:
0 +1 -1 +1 -1
15. // y:
0 +1 -1 +1 -1 -1 +1
16. //
17. // 8 3 5
18. //
\|/
19. // 2-0-1
20. //
/|\
21. // 6 4 7
22.
23. /// *********************
24. /// PREPROCESSOR COMMANDS
25. /// *********************
26.
27. #include // vector containers
28. #include // mathematical library
29. #include // for the use of 'cout'
30. #include // file streams
31. #include // string streams
32. #include // standard library
33. #define SQ(x) ((x) * (x)) // square function; replaces SQ(x) by ((x) * (x))
^y
|
x
--->
in the code
34.
35. using namespace std; // permanently use the standard namespace
36.
37. /// *********************
38. /// SIMULATION PARAMETERS
39. /// *********************
40.
41. // These are the relevant simulation parameters.
42. // They can be changed by the user.
43. // If a bottom or top wall shall move in negative x-direction, a negative ve
locity has to be specified.
44. // Moving walls and gravity can be switched on simultaneously.
45.
46. /// Simulation types
47. // Exactly one of the following options has to be defined.
48. // RIGID_CYLINDER
49. // - for simulation of, e.g., Karman vortex street
50. // - the cylinder is kept in space, it does not move
51. // DEFORMABLE_CYLINDER
52. // - for simulation of a deformable cylinder
53. // - the cylinder moves along with the flow
54. // DEFORMABLE_RBC
55. // - for simulation of a deformable red blood cell
56. // - the cell moves along with the flow
57.
58. //#define RIGID_CYLINDER
59. //#define DEFORMABLE_CYLINDER
60. #define DEFORMABLE_RBC
61.
62. /// Fluid/lattice properties
63.
64. #ifdef RIGID_CYLINDER
65.
const int Nx = 300; // number of lattice nodes along the x-axis (periodic)
66.
const int Ny = 62; // number of lattice nodes along the y-axis (including
two wall nodes)
67.
68.
const double tau = 0.53; // relaxation time
const int t_num = 100000; // number of time steps (running from 1 to t_num)
69.
const int t_disk = 200; // disk write time step (data will be written to t
he disk every t_disk step)
70.
71.
const int t_info = 5000; // info time step (screen message will be printed
every t_info step)
const double gravity = 0.000005; // force density due to gravity (in posit
ive x-direction)
72.
const double wall_vel_bottom = 0; // velocity of the bottom wall (in posit
ive x-direction)
73.
const double wall_vel_top = 0; // velocity of the top wall (in positive x-
direction)
74. #else
75.
const int Nx = 30; // number of lattice nodes along the x-axis (periodic)
76.
const int Ny = 32; // number of lattice nodes along the y-axis (including
two wall nodes)
77.
78.
const double tau = 1; // relaxation time
const int t_num = 50000; // number of time steps (running from 1 to t_num)
79.
const int t_disk = 200; // disk write time step (data will be written to t
he disk every t_disk step)
80.
81.
const int t_info = 1000; // info time step (screen message will be printed
every t_info step)
const double gravity = 0; // force density due to gravity (in positive x-d
irection)
82.
const double wall_vel_bottom = -0.02; // velocity of the bottom wall (in p
ositive x-direction)
83.
const double wall_vel_top = 0.02; // velocity of the top wall (in positive
x-direction)
const int particle_num_nodes = 36; // number of surface nodes
const double particle_radius = 8; // radius
const double particle_stiffness = 0.3; // stiffness modulus
const double particle_center_x = 15; // center position (x-component)
const double particle_center_y = 29; // center position (y-component)
const double particle_center_y = 15; // center position (y-component)
84. #endif
85.
86. /// Particle properties
87.
88. #ifdef RIGID_CYLINDER
89.
90.
91.
92.
93.
94. #else
95.
96.
97.
98.
99.
100.
101. #endif
102.
103. /// *****************
104. /// DECLARE VARIABLES
105. /// *****************
106.
107. // The following code should not be modified when it is first used.
108.
109. const double omega = 1. / tau; // relaxation frequency (inverse of relaxati
const int particle_num_nodes = 32; // number of surface nodes
const double particle_radius = 6; // radius
const double particle_stiffness = 0.1; // stiffness modulus
const double particle_bending = 0.001; // stiffness modulus
const double particle_center_x = 15; // center position (x-component)
on time)
110. double ***pop, ***pop_old; // LBM populations (old and new)
111. double **density; // fluid density
112. double **velocity_x; // fluid velocity (x-component)
113. double **velocity_y; // fluid velocity (y-component)
114. double **force_x; // fluid force (x-component)
115. double **force_y; // fluid force (y-component)
116. double force_latt[9]; // lattice force term entering the lattice Boltzmann
equation
117. double pop_eq[9]; // equilibrium populations
118. const double weight[9] = {4./9., 1./9., 1./9., 1./9., 1./9., 1./36., 1./36.,
1./36., 1./36.}; // lattice weights
119.
120. /// ******************
121. /// PARTICLE STRUCTURE
122. /// ******************
123.
124. // The following code handles the object immersed in the flow.
125. // In the present implementation, only a single object can be put into the
flow.
126.
127. /// Structure for surface nodes
128. // Each node has a current x- and y-position and a reference x- and y-posit
ion.
/// Constructor
node_struct() {
x = 0;
y = 0;
x_ref = 0;
y_ref = 0;
vel_x = 0;
vel_y = 0;
force_x = 0;
force_y = 0;
}
129.
130. struct node_struct {
131.
132.
133.
134.
135.
136.
137.
138.
139.
140.
141.
142.
143.
144.
145.
146.
147.
148.
149.
150.
151.
/// Elements
double x; // current x-position
double y; // current y-position
double x_ref; // reference x-position
double y_ref; // reference y-position
double vel_x; // node velocity (x-component)
/// Constructor
particle_struct() {
double vel_y; // node velocity (y-component)
double force_x; // node force (x-component)
double force_y; // node force (y-component)
152.
153.
154.
155. };
156.
157. /// Structure for object (either cylinder or red blood cell)
158.
159. struct particle_struct {
160.
161.
162.
163.
164.
165.
166.
167.
168.
169.
170.
171.
172.
173.
174.
175.
num_nodes = particle_num_nodes;
radius = particle_radius;
stiffness = particle_stiffness;
center.x = particle_center_x;
center.y = particle_center_y;
center.x_ref = particle_center_x;
center.y_ref = particle_center_y;
node = new node_struct[num_nodes];
// The initial shape of the object is set in the following.
// For a cylinder (rigid or deformable), the nodes define a circle.
// For a red blood cell, the y-position has to be changed in order to d
escribe a red blood cell.
176.
// Initially, the current node positions and reference node positions a
re identical.
// During the simulation, only the current positions are updated,
// the reference node positions are fixed.
for(int n = 0; n < num_nodes; ++n) {
#if defined RIGID_CYLINDER || defined DEFORMABLE_CYLINDER
node[n].x = center.x + radius * sin(2. * M_PI * (double) n / num_no
177.
178.
179.
180.
181.
182.
des);
183.
node[n].x_ref = center.x + radius * sin(2. * M_PI * (double) n / nu
m_nodes);
184.
node[n].y = center.y + radius * cos(2. * M_PI * (double) n / num_no
des);
185.
node[n].y_ref = center.y + radius * cos(2. * M_PI * (double) n / nu
m_nodes);
186.
187.
188.
#endif
#ifdef DEFORMABLE_RBC
191.
192.
193.
194.
195.
196.
198.
199.
200.
node[n].y = center.y + sqrt(1 - SQ((center.x - node[n].x) / radiu
s)) * (0.207 + 2.00 * SQ((center.x - node[n].x) / radius) - 1.12 * SQ(SQ((ce
nter.x - node[n].x) / radius))) * radius / 2;
197.
node[n].y_ref = center.y + sqrt(1 - SQ((center.x - node[n].x) / r
adius)) * (0.207 + 2.00 * SQ((center.x - node[n].x) / radius) - 1.12 * SQ(SQ
((center.x - node[n].x) / radius))) * radius / 2;
}
else {
189.
node[n].x = center.x + radius * sin(2. * M_PI * (double) n / num_no
des);
190.
node[n].x_ref = center.x + radius * sin(2. * M_PI * (double) n / nu
m_nodes);
node[n].y = radius * cos(2. * M_PI * (double) n / num_nodes);
// Parametrization of the red blood cell shape in 2D
if(node[n].y > 0) {
node[n].y = center.y - sqrt(1 - SQ((center.x - node[n].x) / radiu
s)) * (0.207 + 2.00 * SQ((center.x - node[n].x) / radius) - 1.12 * SQ(SQ((ce
nter.x - node[n].x) / radius))) * radius / 2;
201.
node[n].y_ref = center.y - sqrt(1 - SQ((center.x - node[n].x) / r
adius)) * (0.207 + 2.00 * SQ((center.x - node[n].x) / radius) - 1.12 * SQ(SQ
((center.x - node[n].x) / radius))) * radius / 2;
}
#endif
}
}
/// Elements
int num_nodes; // number of surface nodes
double radius; // object radius
double stiffness; // stiffness modulus
node_struct center; // center node
node_struct *node; // list of nodes
202.
203.
204.
205.
206.
207.
208.
209.
210.
211.
212.
213.
214. };
215.
216. /// *****************
217. /// DECLARE FUNCTIONS
218. /// *****************
219.
220. // The following functions are used in the simulation code.
221.
222. void initialize(); // allocate memory and initialize variables
223. void LBM(); // perform LBM operations
224. void momenta(); // compute fluid density and velocity from the populations
225. void equilibrium(double, double, double); // compute the equilibrium popula
tions from the fluid density and velocity
226. void compute_particle_forces(particle_struct); // compute the forces acting
on the object nodes
227. void spread(particle_struct); // spread node forces to fluid lattice
228. void interpolate(particle_struct); // interpolate node velocities from flui
d velocity
229. void update_particle_position(particle_struct); // update object center pos
ition
230. void write_fluid_vtk(int); // write the fluid state to the disk as VTK file
231. void write_particle_vtk(int, particle_struct); // write the particle state
to the disk as VTK file
232. void write_data(int, particle_struct); // write data to the disk (drag/lift,
center position)
233.
234. /// *************
235. /// MAIN FUNCTION
236. /// *************
237.
238. // This is the main function, containing the simulation initialization and
the simulation loop.
initialize(); // allocate memory and initialize variables
particle_struct particle; // create immersed object
/// Compute derived quantities
/// ************
/// PREPARATIONS
/// ************
239.
240. int main() {
241.
242.
243.
244.
245.
246.
247.
248.
249.
250.
251.
252.
253.
const double D = Ny - 2; // inner channel diameter
const double nu = (tau - 0.5) / 3; // lattice viscosity
const double umax = gravity / (2 * nu) * SQ(0.5 * D); // expected maximum
velocity for Poiseuille flow without immersed object
254.
const double Re = D * umax / nu; // Reynolds number for Poiseuille flow w
ithout immersed object
255.
256.
257.
258.
259.
260.
261.
262.
263.
264.
265.
266.
267.
268.
269.
270.
271.
272.
273.
274.
275.
276.
277.
278.
279.
280.
/// Report derived parameters
cout << "simulation parameters" << endl;
cout << "=====================" << endl;
cout << "D = " << D << endl;
cout << "nu = " << nu << endl;
cout << "umax = " << umax << endl;
cout << "Re = " << Re << endl;
cout << endl;
/// ***************
/// SIMULATION LOOP
/// ***************
// Overview of simulation algorithm:
// 1) compute the node forces based on the object's deformation
// 2) spread the node forces to the fluid lattice
// 3) update the fluid state via LBM
// 4) interpolate the fluid velocity to the object nodes
// 5) update node positions and object's center
// 6) if desired, write data to disk and report status
cout << "starting simulation" << endl;
for(int t = 1; t <= t_num; ++t) { // run over all times between 1 and t_n
um
281.
282.
283.
mesh
284.
285.
286.
287.
288.
289.
290.
291.
292.
293.
294.
295.
296.
297.
compute_particle_forces(particle); // compute particle forces
spread(particle); // spread forces from the Lagrangian to the Eulerian
LBM(); // perform collision, propagation, and bounce-back
interpolate(particle); // interpolate velocity
update_particle_position(particle); // update particle position
/// Write fluid and particle to VTK files
// The data is only written each t_info time step.
if(t % t_disk == 0) {
write_fluid_vtk(t);
write_particle_vtk(t, particle);
write_data(t, particle);
}
/// Report end of time step
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
