tFinished velocity projection - ns2dfd - 2D finite difference Navier Stokes solver for fluid dynamics
(HTM) git clone git://src.adamsgaard.dk/ns2dfd
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit 108a9b876d69391731de42da1ab178384198f345
(DIR) parent ccc17b2737056c0104d110fa098e7f330be5c919
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 3 Mar 2014 17:16:17 +0100
Finished velocity projection
Diffstat:
A src/solution.c | 89 +++++++++++++++++++++++++++++++
A src/solution.h | 9 +++++++++
2 files changed, 98 insertions(+), 0 deletions(-)
---
(DIR) diff --git a/src/solution.c b/src/solution.c
t@@ -0,0 +1,89 @@
+#include <math.h>
+
+/* Second order central difference (d^2 u)/(dx^2) */
+inline double ddu_dxx(double **U, int i, int j, double dx)
+{
+ return (U[i+1][j] - 2.0*U[i][j] + U[i-1][j])/(dx*dx);
+}
+
+/* Second order central difference (d^2 v)/(dx^2) */
+inline double ddv_dxx(double **V, int i, int j, double dx)
+{
+ return (V[i+1][j] - 2.0*V[i][j] + V[i-1][j])/(dx*dx);
+}
+
+/* Second order central difference (d^2 u)/(dy^2) */
+inline double ddu_dyy(double **U, int i, int j, double dy)
+{
+ return (U[i][j+1] - 2.0*U[i][j] + U[i][j-1])/(dy*dy);
+}
+
+/* Second order central difference (d^2 v)/(dy^2) */
+inline double ddv_dyy(double **V, int i, int j, double dy)
+{
+ return (V[i][j+1] - 2.0*V[i][j] + V[i][j-1])/(dy*dy);
+}
+
+/* First order upwind difference (d(uu))/(dx) */
+double duu_dx(double **U, int i, int j, double dx, double gamma)
+{
+ return 1.0/dx*(pow((U[i][j] + U[i+1][j])/2.0, 2.0) -
+ pow((U[i-1][j] + U[i][j])/2.0, 2.0))
+ + gamma/dx*(
+ fabs(U[i][j] + U[i+1][j])/2.0 * (U[i][j] - U[i+1][j])/2.0 -
+ fabs(U[i-1][j] + U[i][j])/2.0 * (U[i-1][j] - U[i][j])/2.0);
+}
+
+double dvv_dy(double **V, int i, int j, double dy, double gamma)
+{
+ return 1.0/dy*(pow((V[i][j] + V[i][j+1])/2.0, 2.0) -
+ pow((V[i][j-1] + V[i][j])/2.0, 2.0))
+ + gamma/dy*(
+ fabs(V[i][j] + V[i][j+1])/2.0 * (V[i][j] - V[i][j+1])/2.0 -
+ fabs(V[i][j-1] + V[i][j])/2.0 * (V[i][j-1] - V[i][j])/2.0);
+}
+
+/* First order upwind difference (d(uv))/(dy) */
+double duv_dy(double **U, double **V, int i, int j, double dy, double gamma)
+{
+ return 1.0/dy*(
+ (V[i][j] + V[i+1][j])/2.0 * (U[i][j] + U[i][j+1])/2.0 -
+ (V[i][j-1] + V[i+1][j-1])/2.0 * (U[i][j-1] + U[i][j])/2.0)
+ + gamma/dy*(
+ fabs(V[i][j] + V[i+1][j])/2.0 * (U[i][j] - U[i][j+1])/2.0 -
+ fabs(V[i][j-1] + V[i+1][j-1])/2.0 * (U[i][j-1] - U[i][j])/2.0);
+}
+
+double duv_dx(double **U, double **V, int i, int j, double dx, double gamma)
+{
+ return 1.0/dx*(
+ (U[i][j] + U[i][j+1])/2.0 * (V[i][j] + V[i+1][j])/2.0 -
+ (U[i-1][j] + U[i-1][j+1])/2.0 * (V[i-1][j] + V[i][j])/2.0)
+ + gamma/dx*(
+ fabs(U[i][j] + U[i][j+1])/2.0 * (V[i][j] - V[i+1][j])/2.0 -
+ fabs(U[i-1][j] + U[i-1][j+1])/2.0 * (V[i-1][j] - V[i][j])/2.0);
+}
+
+/* Project new velocities based on Chorin's method, saved in F and G */
+void compute_f_g(double **F, double **G,
+ double **U, double **V, double **P, double re,
+ int nx, int ny, double dx, double dy,
+ double gx, double gy, double dt, double gamma)
+{
+ int i, j;
+ for (i=1; i<=nx; i++) {
+ for (j=1; j<=ny; j++) {
+ F[i][j] = U[i][j]
+ + dt*(1.0/re*(ddu_dxx(U, i, j, dx) + ddu_dyy(U, i, j, dy))
+ - duu_dx(U, i, j, dx, gamma)
+ - duv_dy(U, V, i, j, dy, gamma)
+ + gx);
+
+ G[i][j] = V[i][j]
+ + dt*(1.0/re*(ddv_dxx(V, i, j, dx) + ddv_dyy(V, i, j, dy))
+ - duv_dx(U, V, i, j, dx, gamma)
+ - dvv_dy(V, i, j, dy, gamma)
+ + gy);
+ }
+ }
+}
(DIR) diff --git a/src/solution.h b/src/solution.h
t@@ -0,0 +1,9 @@
+#ifndef SOLUTION_H_
+#define SOLUTION_H_
+
+void compute_f_g(double **F, double **G,
+ double **U, double **V, double **P, double re,
+ int nx, int ny, double dx, double dy,
+ double gx, double gy, double dt, double gamma);
+
+#endif