初学者,编了一个用simple算法解couette程序,为什么收敛不了?
浏览:1005 回答: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"); } }