Skip to content

Commit 73cce2e

Browse files
author
nianhuawong
committed
平头钝体滑移壁面，计算完成，计算会出负，所以加上SolutionFix.
1 parent be3ed66 commit 73cce2e

File tree

3 files changed

+88
-88
lines changed

3 files changed

+88
-88
lines changed

‎src/Compute_Boundary.cpp‎

Lines changed: 70 additions & 75 deletions
Original file line numberDiff line numberDiff line change
@@ -46,31 +46,7 @@ void Compute_Boundary_Blunt_Body()
4646
qField[i][j][IP] = qField[i - 1][j][IP];
4747
}
4848
}
49-
50-
//下边界, slip wall
51-
for (int i = ist; i <= ied; i++)
52-
{
53-
for (int j = jst; j >= 0; j--)
54-
{
55-
qField[i][j][IR] = qField[i][j + 1][IR];
56-
//qField[i][j][IU] = qField[i][j + 1][IU];
57-
//qField[i][j][IV] = qField[i][j + 1][IV];
58-
//qField[i][j][IU] = 0.0;
59-
//qField[i][j][IV] = 0.0;
60-
qField[i][j][IP] = qField[i][j + 1][IP];
61-
}
62-
63-
//以下是滑移物面
64-
qField[i][jst ][IU] = 0.0;
65-
qField[i][jst - 1][IU] = -qField[i][jst + 1][IU];
66-
qField[i][jst - 2][IU] = -qField[i][jst + 1][IU];
67-
68-
qField[i][jst ][IV] = 0.0;
69-
qField[i][jst - 1][IV] = -qField[i][jst + 1][IV];
70-
qField[i][jst - 2][IV] = -qField[i][jst + 1][IV];
71-
}
72-
73-
//左边界：超声速入口 ！！！左边界条件设置放在下边界后面，确保左下角点设置为入口边界
49+
//左边界：超声速入口
7450
for (int i = 0; i <= ist; i++)
7551
{
7652
for (int j = jst; j <= jed; j++)
@@ -82,69 +58,88 @@ void Compute_Boundary_Blunt_Body()
8258
}
8359
}
8460

85-
//左物面, slip wall
86-
for (int j = jst + Jw1; j <= jst + Jw2; j++)
61+
//下边界, slip wall
62+
for (int i = ist; i <= ied; i++)
8763
{
88-
for (int i = ist + Iw; i <= ist + Iw + num_ghost_point; i++)
89-
{
90-
qField[i][j][IR] = qField[i - 1][j][IR];
91-
//qField[i][j][IU] = 0.0;
92-
//qField[i][j][IV] = 0.0;
93-
qField[i][j][IP] = qField[i - 1][j][IP];
94-
}
64+
qField[i][jst][IR] = qField[i][jst + 1][IR]; //weno格式不加这一行在下边界不会反射激波？
65+
//qField[i][jst][IU] = qField[i][jst + 1][IU];
66+
qField[i][jst][IV] = 0.0;
67+
qField[i][jst][IP] = qField[i][jst + 1][IP]; //weno格式不加这一行在下边界不会反射激波？
9568

96-
qField[ist + Iw ][j][IU] = 0.0;
97-
qField[ist + Iw + 1][j][IU] = -qField[ist + Iw - 1][j][IU];
98-
qField[ist + Iw + 2][j][IU] = -qField[ist + Iw - 1][j][IU];
69+
qField[i][jst - 1][IR] = qField[i][jst + 1][IR];
70+
qField[i][jst - 1][IU] = qField[i][jst + 1][IU];
71+
qField[i][jst - 1][IV] = -qField[i][jst + 1][IV];
72+
qField[i][jst - 1][IP] = qField[i][jst + 1][IP];
9973

100-
qField[ist + Iw ][j][IV] = 0.0;
101-
qField[ist + Iw + 1][j][IV] = -qField[ist + Iw - 1][j][IV];
102-
qField[ist + Iw + 2][j][IV] = -qField[ist + Iw - 1][j][IV];
74+
qField[i][jst - 2][IR] = qField[i][jst + 2][IR];
75+
qField[i][jst - 2][IU] = qField[i][jst + 2][IU];
76+
qField[i][jst - 2][IV] = -qField[i][jst + 2][IV];
77+
qField[i][jst - 2][IP] = qField[i][jst + 2][IP];
10378
}
10479

10580
//上物面, slip wall
81+
int JJ = jst + Jw2;
10682
for (int i = ist + Iw; i <= ied; i++)
10783
{
108-
for (int j = jst + Jw2; j >= Jw2; j--)
109-
{
110-
qField[i][j][IR] = qField[i][j + 1][IR];
111-
//qField[i][j][IU] = 0.0;
112-
//qField[i][j][IV] = 0.0;
113-
qField[i][j][IP] = qField[i][j + 1][IP];
114-
}
115-
116-
qField[i][jst + Jw2 ][IU] = 0.0;
117-
qField[i][jst + Jw2 - 1][IU] = -qField[i][jst + Jw2 + 1][IU];
118-
qField[i][jst + Jw2 - 2][IU] = -qField[i][jst + Jw2 + 1][IU];
119-
120-
qField[i][jst + Jw2 ][IV] = 0.0;
121-
qField[i][jst + Jw2 - 1][IV] = -qField[i][jst + Jw2 + 1][IV];
122-
qField[i][jst + Jw2 - 2][IV] = -qField[i][jst + Jw2 + 1][IV];
84+
//qField[i][JJ][IR] = qField[i][JJ + 1][IR];
85+
//qField[i][JJ][IU] = qField[i][JJ + 1][IU];
86+
qField[i][JJ][IV] = 0.0;
87+
//qField[i][JJ][IP] = qField[i][JJ + 1][IP];
88+
89+
qField[i][JJ - 1][IR] = qField[i][JJ + 1][IR];
90+
qField[i][JJ - 1][IU] = qField[i][JJ + 1][IU];
91+
qField[i][JJ - 1][IV] = -qField[i][JJ + 1][IV];
92+
qField[i][JJ - 1][IP] = qField[i][JJ + 1][IP];
93+
94+
qField[i][JJ - 2][IR] = qField[i][JJ + 2][IR];
95+
qField[i][JJ - 2][IU] = qField[i][JJ + 2][IU];
96+
qField[i][JJ - 2][IV] = -qField[i][JJ + 2][IV];
97+
qField[i][JJ - 2][IP] = qField[i][JJ + 2][IP];
12398
}
12499

125100
//下物面, slip wall
101+
JJ = jst + Jw1;
126102
for (int i = ist + Iw; i <= ied; i++)
127103
{
128-
for (int j = jst + Jw1; j <= jst + Jw1 + num_ghost_point; j++)
129-
{
130-
qField[i][j][IR] = qField[i][j - 1][IR];
131-
//qField[i][j][IU] = 0.0;
132-
//qField[i][j][IV] = 0.0;
133-
qField[i][j][IP] = qField[i][j - 1][IP];
134-
}
135-
136-
qField[i][jst + Jw1 ][IU] = 0.0;
137-
qField[i][jst + Jw1 + 1][IU] = -qField[i][jst + Jw1 - 1][IU];
138-
qField[i][jst + Jw1 + 2][IU] = -qField[i][jst + Jw1 - 1][IU];
104+
//qField[i][JJ][IR] = qField[i][JJ - 1][IR];
105+
//qField[i][JJ][IU] = qField[i][JJ - 1][IU];
106+
qField[i][JJ][IV] = 0.0;
107+
//qField[i][JJ][IP] = qField[i][JJ - 1][IP];
108+
109+
qField[i][JJ + 1][IR] = qField[i][JJ - 1][IR];
110+
qField[i][JJ + 1][IU] = qField[i][JJ - 1][IU];
111+
qField[i][JJ + 1][IV] = -qField[i][JJ - 1][IV];
112+
qField[i][JJ + 1][IP] = qField[i][JJ - 1][IP];
113+
114+
qField[i][JJ + 2][IR] = qField[i][JJ - 2][IR];
115+
qField[i][JJ + 2][IU] = qField[i][JJ - 2][IU];
116+
qField[i][JJ + 2][IV] = -qField[i][JJ - 2][IV];
117+
qField[i][JJ + 2][IP] = qField[i][JJ - 2][IP];
118+
}
139119

140-
qField[i][jst + Jw1 ][IV] = 0.0;
141-
qField[i][jst + Jw1 + 1][IV] = -qField[i][jst + Jw1 - 1][IV];
142-
qField[i][jst + Jw1 + 2][IV] = -qField[i][jst + Jw1 - 1][IV];
120+
//左物面, slip wall
121+
int II = ist + Iw;
122+
for (int j = jst + Jw1; j <= jst + Jw2; j++)
123+
{
124+
//qField[II][j][IR] = qField[II - 1][j][IR];
125+
qField[II][j][IU] = 0.0;
126+
//qField[II][j][IV] = 0.0; //qField[II - 1][j][IV];
127+
//qField[II][j][IP] = qField[II - 1][j][IP];
128+
129+
qField[II + 1][j][IR] = qField[II - 1][j][IR];
130+
qField[II + 1][j][IU] = -qField[II - 1][j][IU];
131+
qField[II + 1][j][IV] = qField[II - 1][j][IV];
132+
qField[II + 1][j][IP] = qField[II - 1][j][IP];
133+
134+
qField[II + 2][j][IR] = qField[II - 2][j][IR];
135+
qField[II + 2][j][IU] = -qField[II - 2][j][IU];
136+
qField[II + 2][j][IV] = qField[II - 2][j][IV];
137+
qField[II + 2][j][IP] = qField[II - 2][j][IP];
143138
}
144139

145140
//角点处理，左上，设置参考角点标号
146-
int II = ist + Iw;
147-
int JJ = jst + Jw2;
141+
II = ist + Iw;
142+
JJ = jst + Jw2;
148143
//1点，对角方向
149144
qField[II + 2][JJ - 2][IR] = qField[II][JJ][IR];
150145
qField[II + 2][JJ - 2][IU] = -qField[II][JJ][IU];
@@ -154,12 +149,12 @@ void Compute_Boundary_Blunt_Body()
154149
//2点, x方向取值
155150
qField[II + 1][JJ - 2][IR] = qField[II - 1][JJ - 2][IR];
156151
qField[II + 1][JJ - 2][IU] = -qField[II - 1][JJ - 2][IU];
157-
qField[II + 1][JJ - 2][IV] = -qField[II - 1][JJ - 2][IV];
152+
qField[II + 1][JJ - 2][IV] = qField[II - 1][JJ - 2][IV];
158153
qField[II + 1][JJ - 2][IP] = qField[II - 1][JJ - 2][IP];
159154

160155
//3点, y方向取值
161156
qField[II + 2][JJ - 1][IR] = qField[II + 2][JJ + 1][IR];
162-
qField[II + 2][JJ - 1][IU] = -qField[II + 2][JJ + 1][IU];
157+
qField[II + 2][JJ - 1][IU] = qField[II + 2][JJ + 1][IU];
163158
qField[II + 2][JJ - 1][IV] = -qField[II + 2][JJ + 1][IV];
164159
qField[II + 2][JJ - 1][IP] = qField[II + 2][JJ + 1][IP];
165160

@@ -183,12 +178,12 @@ void Compute_Boundary_Blunt_Body()
183178
//2点, x方向取值
184179
qField[II + 1][JJ + 2][IR] = qField[II - 1][JJ + 2][IR];
185180
qField[II + 1][JJ + 2][IU] = -qField[II - 1][JJ + 2][IU];
186-
qField[II + 1][JJ + 2][IV] = -qField[II - 1][JJ + 2][IV];
181+
qField[II + 1][JJ + 2][IV] = qField[II - 1][JJ + 2][IV];
187182
qField[II + 1][JJ + 2][IP] = qField[II - 1][JJ + 2][IP];
188183

189184
//3点, y方向取值
190185
qField[II + 2][JJ + 1][IR] = qField[II + 2][JJ - 1][IR];
191-
qField[II + 2][JJ + 1][IU] = -qField[II + 2][JJ - 1][IU];
186+
qField[II + 2][JJ + 1][IU] = qField[II + 2][JJ - 1][IU];
192187
qField[II + 2][JJ + 1][IV] = -qField[II + 2][JJ - 1][IV];
193188
qField[II + 2][JJ + 1][IP] = qField[II + 2][JJ - 1][IP];
194189

‎src/Global.cpp‎

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -62,22 +62,22 @@ void Input_Parameters()
6262
}
6363
else if (global_case_id == 2)
6464
{
65-
num_grid_point_x = 1 * 200 + 1;
66-
num_grid_point_y = 1 * 100 + 1;
65+
num_grid_point_x = 4 * 100 + 1;
66+
num_grid_point_y = 4 * 50 + 1;
6767

68-
max_num_of_steps = 2000;
68+
max_num_of_steps = 5000;
6969
residual_output_steps = 2; //残差输出间隔步数
70-
flow_save_steps = 50; //流场输出间隔步数
70+
flow_save_steps = 250; //流场输出间隔步数
7171
converge_criterion = 1e-8; //残差收敛标准
7272
tec_file_name = "./results/flow.plt";
7373

74-
cfl_num = 0.3;
75-
max_simu_time = 2.0;
74+
cfl_num = 0.4;
75+
max_simu_time = 10;
7676

7777
method_of_half_q = 1; //1-MUSCL, 2-WENO(不插值), 3-WCNS
7878
muscl_k = 1.0/3; //0.0-二阶迎风偏置， 1/3-二阶迎风偏置
7979
method_of_limiter = 1; //0-nolim, 1-vanleer, 2-minmod, 3-superbee, 4-1st
80-
method_of_flux = 1; //1-Roe, 2-Steger, 3-VanLeer, 4-WENO, 5-WCNS
80+
method_of_flux = 2; //1-Roe, 2-Steger, 3-VanLeer, 4-WENO, 5-WCNS
8181
entropy_fix_coeff = 0.1; //Roe格式熵修正系数epsilon
8282
}
8383
}
@@ -92,11 +92,11 @@ void Init_Global_Param()
9292
time_step = 0.0; //时间步长要根据最大特征值确定，这里只是初始化
9393
solve_direction = 'x';
9494

95-
num_of_RK_stages = 3;
96-
RK_Coeff = { {1.0, 0.0, 1.0},{3.0 / 4.0, 1.0 / 4.0, 1.0 / 4.0},{1.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0} };
95+
//num_of_RK_stages = 3;
96+
//RK_Coeff = { {1.0, 0.0, 1.0},{3.0 / 4.0, 1.0 / 4.0, 1.0 / 4.0},{1.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0} };
9797

98-
//num_of_RK_stages = 2;
99-
//RK_Coeff = { {1.0, 0.0, 1.0},{0.5, 0.5, 0.5} };
98+
num_of_RK_stages = 2;
99+
RK_Coeff = { {1.0, 0.0, 1.0},{0.5, 0.5, 0.5} };
100100

101101
//设置运行参数
102102
Input_Parameters();

‎src/Time_Integral.cpp‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
#include "QlQr_Solver.h"
77
#include "Time_Integral.h"
88
#include <cmath>
9+
#include <iostream>
910

1011
void Time_Integration()
1112
{
@@ -81,6 +82,7 @@ void Update_Flowfield_X(int iStage)
8182
#pragma omp for
8283
#endif
8384
for (int i = 1; i <= ied + 1; i++)
85+
//for (int i = ist; i <= ied; i++)
8486
{
8587
for (int j = jst; j <= jed; j++)
8688
{
@@ -104,7 +106,7 @@ void Update_Flowfield_X(int iStage)
104106
if (qPrimitive1[IR] < minimumDensityLimit || qPrimitive1[IP] < minimumPressureLimit)
105107
{
106108
int kkk = 1;
107-
//SolutionFix(qPrimitive1, i, j);
109+
SolutionFix(qPrimitive1, i, j);
108110
}
109111
//RK公式里左端项，q1、q2、q3，即下一stage的q值，还要继续用该值计算rhs(q1)、rhs(q2)
110112
qField_N1[i][j] = qPrimitive1;
@@ -150,6 +152,7 @@ void Update_Flowfield_Y(int iStage)
150152
for (int i = ist; i <= ied; i++)
151153
{
152154
for (int j = 1; j <= jed + 1; j++)
155+
//for (int j = jst; j <= jed; j++)
153156
{
154157
if (marker[i][j] == 0) continue;
155158

@@ -171,7 +174,7 @@ void Update_Flowfield_Y(int iStage)
171174
if (qPrimitive1[IR] < minimumDensityLimit || qPrimitive1[IP] < minimumPressureLimit)
172175
{
173176
int kkk = 1;
174-
//SolutionFix(qPrimitive1, i, j);
177+
SolutionFix(qPrimitive1, i, j);
175178
}
176179
//RK公式里左端项，q1、q2、q3，即下一stage的q值，还要继续用该值计算rhs(q1)、rhs(q2)
177180
qField_N1[i][j] = qPrimitive1;
@@ -190,6 +193,8 @@ void Set_Field()
190193

191194
void SolutionFix(VDouble& primitiveVector, int iCell, int jCell)
192195
{
196+
//cout << "压力出负！i = " << iCell << " j = " << jCell << endl;
197+
193198
for (int iEquation = 0; iEquation < num_of_prim_vars; ++iEquation)
194199
{
195200
primitiveVector[iEquation] = 0.0;

0 commit comments

Comments
 (0)