tsolution.c - ns2dfd - 2D finite difference Navier Stokes solver for fluid dynamics
HTML git clone git://src.adamsgaard.dk/ns2dfd
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tsolution.c (7078B)
---
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4
5 /* Second order central difference (d^2 u)/(dx^2) */
6 inline double ddu_dxx(double **U, int i, int j, double dx)
7 {
8 return (U[i+1][j] - 2.0*U[i][j] + U[i-1][j])/(dx*dx);
9 }
10
11 /* Second order central difference (d^2 v)/(dx^2) */
12 inline double ddv_dxx(double **V, int i, int j, double dx)
13 {
14 return (V[i+1][j] - 2.0*V[i][j] + V[i-1][j])/(dx*dx);
15 }
16
17 /* Second order central difference (d^2 u)/(dy^2) */
18 inline double ddu_dyy(double **U, int i, int j, double dy)
19 {
20 return (U[i][j+1] - 2.0*U[i][j] + U[i][j-1])/(dy*dy);
21 }
22
23 /* Second order central difference (d^2 v)/(dy^2) */
24 inline double ddv_dyy(double **V, int i, int j, double dy)
25 {
26 return (V[i][j+1] - 2.0*V[i][j] + V[i][j-1])/(dy*dy);
27 }
28
29 /* First order upwind difference (d(uu))/(dx) */
30 double duu_dx(double **U, int i, int j, double dx, double gamma)
31 {
32 return 1.0/dx*(pow((U[i][j] + U[i+1][j])/2.0, 2.0) -
33 pow((U[i-1][j] + U[i][j])/2.0, 2.0))
34 + gamma/dx*(
35 fabs(U[i][j] + U[i+1][j])/2.0 * (U[i][j] - U[i+1][j])/2.0 -
36 fabs(U[i-1][j] + U[i][j])/2.0 * (U[i-1][j] - U[i][j])/2.0);
37 }
38
39 double dvv_dy(double **V, int i, int j, double dy, double gamma)
40 {
41 return 1.0/dy*(pow((V[i][j] + V[i][j+1])/2.0, 2.0) -
42 pow((V[i][j-1] + V[i][j])/2.0, 2.0))
43 + gamma/dy*(
44 fabs(V[i][j] + V[i][j+1])/2.0 * (V[i][j] - V[i][j+1])/2.0 -
45 fabs(V[i][j-1] + V[i][j])/2.0 * (V[i][j-1] - V[i][j])/2.0);
46 }
47
48 /* First order upwind difference (d(uv))/(dy) */
49 double duv_dy(double **U, double **V, int i, int j, double dy, double gamma)
50 {
51 return 1.0/dy*(
52 (V[i][j] + V[i+1][j])/2.0 * (U[i][j] + U[i][j+1])/2.0 -
53 (V[i][j-1] + V[i+1][j-1])/2.0 * (U[i][j-1] + U[i][j])/2.0)
54 + gamma/dy*(
55 fabs(V[i][j] + V[i+1][j])/2.0 * (U[i][j] - U[i][j+1])/2.0 -
56 fabs(V[i][j-1] + V[i+1][j-1])/2.0 * (U[i][j-1] - U[i][j])/2.0);
57 }
58
59 double duv_dx(double **U, double **V, int i, int j, double dx, double gamma)
60 {
61 return 1.0/dx*(
62 (U[i][j] + U[i][j+1])/2.0 * (V[i][j] + V[i+1][j])/2.0 -
63 (U[i-1][j] + U[i-1][j+1])/2.0 * (V[i-1][j] + V[i][j])/2.0)
64 + gamma/dx*(
65 fabs(U[i][j] + U[i][j+1])/2.0 * (V[i][j] - V[i+1][j])/2.0 -
66 fabs(U[i-1][j] + U[i-1][j+1])/2.0 * (V[i-1][j] - V[i][j])/2.0);
67 }
68
69 /* Project new velocities based on Chorin's method, saved in F and G */
70 void compute_f_g(double **F, double **G,
71 double **U, double **V, double **P, double re,
72 int nx, int ny, double dx, double dy,
73 double gx, double gy, double dt, double gamma)
74 {
75 int i, j;
76 for (i=1; i<nx; i++)
77 for (j=1; j<=ny; j++)
78 F[i][j] = U[i][j]
79 + dt*(1.0/re*(ddu_dxx(U, i, j, dx) + ddu_dyy(U, i, j, dy))
80 - duu_dx(U, i, j, dx, gamma)
81 - duv_dy(U, V, i, j, dy, gamma)
82 + gx);
83
84 for (i=1; i<=nx; i++)
85 for (j=1; j<ny; j++)
86 G[i][j] = V[i][j]
87 + dt*(1.0/re*(ddv_dxx(V, i, j, dx) + ddv_dyy(V, i, j, dy))
88 - duv_dx(U, V, i, j, dx, gamma)
89 - dvv_dy(V, i, j, dy, gamma)
90 + gy);
91
92 for (j=1; j<=ny; j++) {
93 F[0][j] = U[0][j];
94 F[nx][j] = U[nx][j];
95 }
96
97 for (i=1; i<=nx; i++) {
98 G[i][0] = V[i][0];
99 G[i][ny] = V[i][ny];
100 }
101 }
102
103 /* Find the right hand side value of the Poisson equation */
104 void compute_rhs(double **RHS, double **F, double **G,
105 double dt, double dx, double dy, int nx, int ny)
106 {
107 int i, j;
108 for (i=1; i<=nx; i++) {
109 for (j=1; j<=ny; j++) {
110 RHS[i][j] = 1.0/dt*(
111 (F[i][j] - F[i-1][j])/dx + (G[i][j] - G[i][j-1])/dy);
112 }
113 }
114 }
115
116 void e_parameters(double *e_left, double *e_right,
117 double *e_bottom, double *e_top,
118 int i, int j, int nx, int ny)
119 {
120 if (i == 1) {
121 *e_left = 0.0;
122 } else if (i > 1) {
123 *e_left = 1.0;
124 } else {
125 fprintf(stderr, "e_parameters: i value %d not understood\n", i);
126 exit(1);
127 }
128
129 if (i < nx) {
130 *e_right = 1.0;
131 } else if (i == nx) {
132 *e_right = 0.0;
133 } else {
134 fprintf(stderr, "e_parameters: i value %d not understood\n", i);
135 exit(1);
136 }
137
138 if (j == 1) {
139 *e_bottom = 0.0;
140 } else if (j > 1) {
141 *e_bottom = 1.0;
142 } else {
143 fprintf(stderr, "e_parameters: j value %d not understood\n", j);
144 exit(1);
145 }
146
147 if (j < ny) {
148 *e_top = 1.0;
149 } else if (j == ny) {
150 *e_top = 0.0;
151 } else {
152 fprintf(stderr, "e_parameters: j value %d not understood\n", j);
153 exit(1);
154 }
155 }
156
157 /* Compute a successive overrelaxation (SOR) cycle */
158 void sor_cycle(double **P, double **RHS, double omega,
159 int nx, int ny, double dx, double dy)
160 {
161 int i, j;
162 double e_left, e_right, e_top, e_bottom;
163
164 for (i=1; i<=nx; i++) {
165 for (j=1; j<=ny; j++) {
166 e_parameters(&e_left, &e_right, &e_bottom, &e_top, i, j, nx, ny);
167 P[i][j] = (1.0 - omega)*P[i][j] +
168 omega/((e_right + e_left)/(dx*dx) + (e_top + e_bottom)/(dy*dy))
169 *( (e_right*P[i+1][j] + e_left*P[i-1][j])/(dx*dx) +
170 (e_top*P[i][j+1] + e_bottom*P[i][j-1])/(dy*dy) -
171 RHS[i][j]);
172 }
173 }
174 }
175
176 void calculate_residuals(double **RES, double **P, double **RHS,
177 int nx, int ny, double dx, double dy)
178 {
179 int i, j;
180 double e_left, e_right, e_top, e_bottom;
181
182 for (i=1; i<=nx; i++) {
183 for (j=1; j<=ny; j++) {
184 e_parameters(&e_left, &e_right, &e_bottom, &e_top, i, j, nx, ny);
185 RES[i][j] =
186 (e_right*(P[i+1][j] - P[i][j]) - e_left*(P[i][j] - P[i-1][j]))
187 /(dx*dx) +
188 (e_top*(P[i][j+1] - P[i][j]) - e_bottom*(P[i][j] - P[i][j-1]))
189 /(dy*dy) - RHS[i][j];
190 if (RES[i][j] != RES[i][j]) {
191 fprintf(stderr, "calculate_residuals: "
192 "Error, residual is %f in cell i=%d, j=%d\n",
193 RES[i][j], i, j);
194 exit(1);
195 }
196 }
197 }
198 }
199
200 double max_residual_norm(double **RES, double **P, double **RHS,
201 int nx, int ny, double dx, double dy)
202 {
203 int i, j;
204 double val;
205 double max = -1.0e9;
206
207 calculate_residuals(RES, P, RHS, nx, ny, dx, dy);
208
209 for (i=1; i<=nx; i++) {
210 for (j=1; j<=ny; j++) {
211 val = fabs(RES[i][j]);
212 if (val > max)
213 max = val;
214 }
215 }
216 return max;
217 }
218
219 void correct_velocities(double **U, double **V, double **F, double **G,
220 double **P, int nx, int ny, double dt, double dx, double dy)
221 {
222 int i, j;
223 for (i=1; i<nx; i++)
224 for (j=1; j<=ny; j++)
225 U[i][j] = F[i][j] - dt/dx*(P[i+1][j] - P[i][j]);
226
227 for (i=1; i<=nx; i++)
228 for (j=1; j<ny; j++)
229 V[i][j] = G[i][j] - dt/dy*(P[i][j+1] - P[i][j]);
230 }