HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
NavierStokes2DParabolicFunction.c File Reference

Compute the viscous terms for the 2D Navier Stokes equations. More...

#include <stdlib.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <physicalmodels/navierstokes2d.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions

int NavierStokes2DParabolicFunction (double *par, double *u, void *s, void *m, double t)
 

Detailed Description

Compute the viscous terms for the 2D Navier Stokes equations.

Author
Debojyoti Ghosh

Definition in file NavierStokes2DParabolicFunction.c.

Function Documentation

int NavierStokes2DParabolicFunction ( double *  par,
double *  u,
void *  s,
void *  m,
double  t 
)

Compute the viscous terms in the 2D Navier Stokes equations: this function computes the following:

\begin{equation} \frac {\partial} {\partial x} \left[\begin{array}{c} 0 \\ \tau_{xx} \\ \tau_{yx} \\ u \tau_{xx} + v \tau_{yx} - q_x \end{array}\right] + \frac {\partial} {\partial y} \left[\begin{array}{c} 0 \\ \tau_{xy} \\ \tau_{yy} \\ u \tau_{xy} + v \tau_{yy} - q_y \end{array}\right] \end{equation}

where

\begin{align} \tau_{xx} &= \frac{2}{3}\left(\frac{\mu}{Re}\right)\left(2\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y}\right),\\ \tau_{xy} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right),\\ \tau_{yx} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right),\\ \tau_{yy} &= \frac{2}{3}\left(\frac{\mu}{Re}\right)\left(-\frac{\partial u}{\partial x} +2\frac{\partial v}{\partial y}\right),\\ q_x &= -\frac{mu}{\left(\gamma-1\right)Re Pr}\frac{\partial T}{\partial x}, \\ q_y &= -\frac{mu}{\left(\gamma-1\right)Re Pr}\frac{\partial T}{\partial y} \end{align}

and the temperature is \(T = \gamma p/\rho\). \(Re\) and \(Pr\) are the Reynolds and Prandtl numbers, respectively. Note that this function computes the entire parabolic term, and thus bypasses HyPar's parabolic function calculation interfaces. NavierStokes2DInitialize() assigns this function to HyPar::ParabolicFunction.

Reference:

  • Tannehill, Anderson and Pletcher, Computational Fluid Mechanics and Heat Transfer, Chapter 5, Section 5.1.7 (However, the non-dimensional velocity and the Reynolds number is based on speed of sound, instead of the freestream velocity).
Parameters
parArray to hold the computed viscous terms
uSolution vector array
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent simulation time

Definition at line 38 of file NavierStokes2DParabolicFunction.c.

45 {
46  HyPar *solver = (HyPar*) s;
47  MPIVariables *mpi = (MPIVariables*) m;
48  NavierStokes2D *physics = (NavierStokes2D*) solver->physics;
49  int i,j,v;
50  _DECLARE_IERR_;
51 
52  int ghosts = solver->ghosts;
53  int imax = solver->dim_local[0];
54  int jmax = solver->dim_local[1];
55  int *dim = solver->dim_local;
56  int nvars = solver->nvars;
57  int ndims = solver->ndims;
58  int size = (imax+2*ghosts)*(jmax+2*ghosts)*nvars;
59 
60  _ArraySetValue_(par,size,0.0);
61  if (physics->Re <= 0) return(0); /* inviscid flow */
62  solver->count_par++;
63 
64  static double two_third = 2.0/3.0;
65  double inv_gamma_m1 = 1.0 / (physics->gamma-1.0);
66  double inv_Re = 1.0 / physics->Re;
67  double inv_Pr = 1.0 / physics->Pr;
68 
69  double *Q; /* primitive variables */
70  Q = (double*) calloc (size,sizeof(double));
71  for (i=-ghosts; i<(imax+ghosts); i++) {
72  for (j=-ghosts; j<(jmax+ghosts); j++) {
73  int p,index[2]; index[0]=i; index[1]=j;
74  double energy,pressure;
75  _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
76  _NavierStokes2DGetFlowVar_( (u+p),Q[p+0],Q[p+1],Q[p+2],energy,
77  pressure,physics->gamma);
78  Q[p+3] = physics->gamma*pressure/Q[p+0]; /* temperature */
79  }
80  }
81 
82  double *QDerivX = (double*) calloc (size,sizeof(double));
83  double *QDerivY = (double*) calloc (size,sizeof(double));
84 
85  IERR solver->FirstDerivativePar(QDerivX,Q,_XDIR_,1,solver,mpi); CHECKERR(ierr);
86  IERR solver->FirstDerivativePar(QDerivY,Q,_YDIR_,1,solver,mpi); CHECKERR(ierr);
87 
88  IERR MPIExchangeBoundariesnD(solver->ndims,solver->nvars,dim,
89  solver->ghosts,mpi,QDerivX); CHECKERR(ierr);
90  IERR MPIExchangeBoundariesnD(solver->ndims,solver->nvars,dim,
91  solver->ghosts,mpi,QDerivY); CHECKERR(ierr);
92 
93  for (i=-ghosts; i<(imax+ghosts); i++) {
94  for (j=-ghosts; j<(jmax+ghosts); j++) {
95  int p,index[2]; index[0]=i; index[1]=j;
96  double dxinv, dyinv;
97  _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
98  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
99  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
100  _ArrayScale1D_((QDerivX+p),dxinv,nvars);
101  _ArrayScale1D_((QDerivY+p),dyinv,nvars);
102  }
103  }
104 
105  double *FViscous = (double*) calloc (size,sizeof(double));
106  double *FDeriv = (double*) calloc (size,sizeof(double));
107 
108  /* Along X */
109  for (i=-ghosts; i<(imax+ghosts); i++) {
110  for (j=-ghosts; j<(jmax+ghosts); j++) {
111  int p,index[2]; index[0]=i; index[1]=j;;
112  _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
113 
114  double uvel, vvel, T, Tx, ux, uy, vx, vy;
115  uvel = (Q+p)[1];
116  vvel = (Q+p)[2];
117  T = (Q+p)[3];
118  Tx = (QDerivX+p)[3];
119  ux = (QDerivX+p)[1];
120  vx = (QDerivX+p)[2];
121  uy = (QDerivY+p)[1];
122  vy = (QDerivY+p)[2];
123 
124  /* calculate viscosity coeff based on Sutherland's law */
125  double mu = raiseto(T, 0.76);
126 
127  double tau_xx, tau_xy, qx;
128  tau_xx = two_third * (mu*inv_Re) * (2*ux - vy);
129  tau_xy = (mu*inv_Re) * (uy + vx);
130  qx = ( (mu*inv_Re) * inv_gamma_m1 * inv_Pr ) * Tx;
131 
132  (FViscous+p)[0] = 0.0;
133  (FViscous+p)[1] = tau_xx;
134  (FViscous+p)[2] = tau_xy;
135  (FViscous+p)[3] = uvel*tau_xx + vvel*tau_xy + qx;
136  }
137  }
138  IERR solver->FirstDerivativePar(FDeriv,FViscous,_XDIR_,-1,solver,mpi); CHECKERR(ierr);
139  for (i=0; i<imax; i++) {
140  for (j=0; j<jmax; j++) {
141  int p,index[2]; index[0]=i; index[1]=j;
142  double dxinv;
143  _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
144  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
145  for (v=0; v<nvars; v++) (par+p)[v] += (dxinv * (FDeriv+p)[v] );
146  }
147  }
148 
149  /* Along Y */
150  for (i=-ghosts; i<(imax+ghosts); i++) {
151  for (j=-ghosts; j<(jmax+ghosts); j++) {
152  int p,index[2]; index[0]=i; index[1]=j;
153  _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
154 
155  double uvel, vvel, T, Ty, ux, uy, vx, vy;
156  uvel = (Q+p)[1];
157  vvel = (Q+p)[2];
158  T = (Q+p)[3];
159  Ty = (QDerivY+p)[3];
160  ux = (QDerivX+p)[1];
161  vx = (QDerivX+p)[2];
162  uy = (QDerivY+p)[1];
163  vy = (QDerivY+p)[2];
164 
165  /* calculate viscosity coeff based on Sutherland's law */
166  double mu = raiseto(T, 0.76);
167 
168  double tau_yx, tau_yy, qy;
169  tau_yx = (mu*inv_Re) * (uy + vx);
170  tau_yy = two_third * (mu*inv_Re) * (-ux + 2*vy);
171  qy = ( (mu*inv_Re) * inv_gamma_m1 * inv_Pr ) * Ty;
172 
173  (FViscous+p)[0] = 0.0;
174  (FViscous+p)[1] = tau_yx;
175  (FViscous+p)[2] = tau_yy;
176  (FViscous+p)[3] = uvel*tau_yx + vvel*tau_yy + qy;
177  }
178  }
179  IERR solver->FirstDerivativePar(FDeriv,FViscous,_YDIR_,-1,solver,mpi); CHECKERR(ierr);
180  for (i=0; i<imax; i++) {
181  for (j=0; j<jmax; j++) {
182  int p,index[2]; index[0]=i; index[1]=j;
183  double dyinv;
184  _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
185  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
186  for (v=0; v<nvars; v++) (par+p)[v] += (dyinv * (FDeriv+p)[v] );
187  }
188  }
189 
190  free(Q);
191  free(QDerivX);
192  free(QDerivY);
193  free(FViscous);
194  free(FDeriv);
195 
196  return(0);
197 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
raiseto
#define raiseto(x, a)
Definition: math_ops.h:37
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
_ArrayScale1D_
#define _ArrayScale1D_(x, a, size)
Definition: arrayfunctions.h:190
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
HyPar::count_par
int count_par
Definition: hypar.h:418
HyPar::FirstDerivativePar
int(* FirstDerivativePar)(double *, double *, int, int, void *, void *)
Definition: hypar.h:243
HyPar::nvars
int nvars
Definition: hypar.h:29
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
_NavierStokes2DGetFlowVar_
#define _NavierStokes2DGetFlowVar_(u, rho, vx, vy, e, P, gamma)
Definition: navierstokes2d.h:81
HyPar::dxinv
double * dxinv
Definition: hypar.h:110
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
HyPar::ndims
int ndims
Definition: hypar.h:26
IERR
#define IERR
Definition: basic.h:16
MPIExchangeBoundariesnD
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Definition: MPIExchangeBoundariesnD.c:42
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17