初学者,编了一个用simple算法解couette程序,为什么收敛不了?
浏览:1224 回答:1
#include<stdio.h>
#include<stdlib.h>
#define UE 1
#define rou 0.002377
#define miu 3.737e-7
#define dx 0.025
#define dy 0.001
#define dt 0.001
#define sor 0.1
int i, j, k, steps=0, steps_sum;
double p[22][12] = { 0 }, u[23][12] = { 0 }, v[24][13] = { 0 };
double mass_u[23][12] = { 0 }, mass_v[24][13] = { 0 };
void show()
{
printf("\n************************************%d************************************\n",steps);
printf("\nu\n");
for (j = 11; j >=1; j--)
{
for (i = 1; i <= 22; i++)
printf("%.3lf ", u[i][j]);
printf("\n");
}
printf("\nv\n");
for (j = 12; j >= 1; j--)
{
for (i = 1; i <= 23; i++)
printf("%.3lf ", v[i][j]);
printf("\n");
}
printf("\np\n");
for (j = 11; j >= 1; j--)
{
for (i = 1; i <= 21; i++)
printf("%.3lf ", v[i][j]);
printf("\n");
}
/* printf("\nmass_u\n");
for (j = 11; j >= 1; j--)
{
for (i = 1; i <= 22; i++)
printf("%.3lf ", mass_u[i][j]);
printf("\n");
}
printf("\nmass_v\n");
for (j = 11; j >= 1; j--)
{
for (i = 1; i <= 22; i++)
printf("%.3lf ", mass_v[i][j]);
printf("\n");
}*/
}
void initialize()
{
for (i = 1; i <= 22; i++)
u[i][11] = UE;
v[15][5] = 0; //disturbancr
for (i = 1; i <= 22; i++)
mass_u[i][11] = rou*UE;
mass_v[15][5] = rou*v[15][5];
}
void compute_mass()
{
double A, B, v_1, v_2, u_1, u_2;
for (i = 2; i <= 21; i++)
for (j = 2; j <= 10; j++)
{
v_1 = 0.5*(v[i][j + 1] + v[i + 1][j + 1]);
v_2 = 0.5*(v[i][j] + v[i + 1][j]);
A = -((rou*u[i + 1][j] * u[i + 1][j] - rou*u[i - 1][j] * u[i - 1][j]) / (2 * dx) + (rou*u[i][j + 1] * v_1 - rou*u[i][j - 1] * v_2) / (2 * dy))
+ miu*((u[i + 1][j] - 2 * u[i][j] + u[i - 1][j]) / (dx*dx) + (u[i][j + 1] - 2 * u[i][j] + u[i][j - 1]) / (dy*dy));
mass_u[i][j] = rou*u[i][j] + A*dt - (dt / dx)*(p[i][j] - p[i - 1][j]);
}
for (i = 2; i <= 22; i++)
for (j = 2; j <= 11; j++)
{
u_1 = 0.5*(u[i][j - 1] + u[i][j]);
u_2 = 0.5*(u[i - 1][j - 1] + u[i - 1][j]);
B = -((rou*v[i + 1][j] * u_1 - rou*v[i - 1][j] * u_2) / (2 * dx) + (rou*v[i][j + 1] * v[i][j + 1] - rou*v[i][j - 1] * v[i][j - 1]) / (2 * dy))
+ miu*((v[i + 1][j] - 2 * v[i][j] + v[i - 1][j]) / (dx*dx) + (v[i][j + 1] - 2 * v[i][j] + v[i][j - 1]) / (dy*dy));
mass_v[i][j] = rou*v[i][j] + B*dt - (dt / dy)*(p[i - 1][j] - p[i - 1][j - 1]);
}
}
void pressure_correct()
{
double a, b, c;
double d[22][12] = { 0 };
double _p[22][12] = { 0 };//pressure correct
double _p_;//mid
a = 2 * (dt / (dx*dx) + dt / (dy*dy));
b = -dt / (dx*dx);
c = -dt / (dy*dy);
for (i = 2; i <= 20; i++)
for (j = 2; j <= 10; j++)
d[i][j] = (mass_u[i + 1][j] - mass_u[i][j]) / dx + (mass_v[i + 1][i + 1] - mass_v[i + 1][j]) / dy;
for (k = 1; k <= 100; k++)
{
for (i = 2; i <= 20; i++)
for (j = 2; j <= 10; j++)
{
_p_ = -(b*_p[i + 1][j] + b*_p[i - 1][j] + c*_p[i][j + 1] + c*_p[i][j - 1] + d[i][j]) / a;
_p[i][j] = _p[i][j] + 0.1*(_p_ - _p[i][j]);
}
}
for (i = 2; i <= 20; i++)
for (j = 2; j <= 20; j++)
p[i][j] = p[i][j]+(sor*_p[i][j]);
}
void compute_velocity()
{
for (i = 2; i <= 21; i++)
for (j = 2; j <= 10; j++)
u[i][j] = mass_u[i][j] / rou;
for (i = 2; i <= 22; i++)
for (j = 2; j <= 11; j++)
v[i][j] = mass_v[i][j] / rou;
for (j = 1; j <= 11; j++)
{
u[1][j] = u[2][j]; u[22][j] = u[21][j];
}
for (j = 1; j <= 12; j++)
v[23][j] = v[22][j];
}
void main()
{
double store[6][12] = { 0 };
int m = 1;
initialize();
show();
steps_sum = 20;
for (steps = 1; steps <= steps_sum; steps++)
{
compute_mass();
compute_velocity();
pressure_correct(); show();
if (steps == 4 || steps == 20 || steps == 50 || steps == 150 || steps == 300)
{
show();
for (j = 1; j <= 11; j++)
store[m][j] = u[15][j];
m += 1;
}
}
printf("\n");
for (j = 11; j >= 1; j--)
{
for (m = 1; m <= 5; m++)
printf("%.9lf ", store[m][j]);
printf("\n");
}
}



















