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
NavierStokes3DParabolicFunction.c File Reference

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

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

Go to the source code of this file.

Functions

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

Detailed Description

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

Author
Debojyoti Ghosh

Definition in file NavierStokes3DParabolicFunction.c.

Function Documentation

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

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

\begin{equation} \frac {\partial} {\partial x} \left[\begin{array}{c} 0 \\ \tau_{xx} \\ \tau_{yx} \\ \tau_{zx} \\ u \tau_{xx} + v \tau_{yx} + w \tau_{zx} - q_x \end{array}\right] + \frac {\partial} {\partial y} \left[\begin{array}{c} 0 \\ \tau_{xy} \\ \tau_{yy} \\ \tau_{zy} \\ u \tau_{xy} + v \tau_{yy} + w \tau_{zy} - q_y \end{array}\right] + \frac {\partial} {\partial z} \left[\begin{array}{c} 0 \\ \tau_{xz} \\ \tau_{yz} \\ \tau_{zz} \\ u \tau_{xz} + v \tau_{yz} + w \tau_{zz} - q_z \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} - \frac{\partial w}{\partial z}\right),\\ \tau_{xy} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right),\\ \tau_{xz} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial z} + \frac{\partial w}{\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} - \frac{\partial w}{\partial z}\right),\\ \tau_{yz} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial w}{\partial y} + \frac{\partial v}{\partial z}\right),\\ \tau_{zx} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial z} + \frac{\partial w}{\partial x}\right),\\ \tau_{zy} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}\right),\\ \tau_{zz} &= \frac{2}{3}\left(\frac{\mu}{Re}\right)\left(-\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y} + 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}, \\ q_z &= -\frac{mu}{\left(\gamma-1\right)Re Pr}\frac{\partial T}{\partial z} \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. NavierStokes3DInitialize() 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 50 of file NavierStokes3DParabolicFunction.c.

57 {
58  HyPar *solver = (HyPar*) s;
59  MPIVariables *mpi = (MPIVariables*) m;
60  NavierStokes3D *physics = (NavierStokes3D*) solver->physics;
61  int i,j,k,v;
62  _DECLARE_IERR_;
63 
64  int ghosts = solver->ghosts;
65  int imax = solver->dim_local[0];
66  int jmax = solver->dim_local[1];
67  int kmax = solver->dim_local[2];
68  int *dim = solver->dim_local;
69  int size = solver->npoints_local_wghosts;
70 
71  _ArraySetValue_(par,size*_MODEL_NVARS_,0.0);
72  if (physics->Re <= 0) return(0); /* inviscid flow */
73  solver->count_par++;
74 
75  static double two_third = 2.0/3.0;
76  double inv_gamma_m1 = 1.0 / (physics->gamma-1.0);
77  double inv_Re = 1.0 / physics->Re;
78  double inv_Pr = 1.0 / physics->Pr;
79 
80 #if defined(CPU_STAT)
81  clock_t startEvent, stopEvent;
82 #endif
83 
84  double *Q; /* primitive variables */
85  Q = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
86 
87 #if defined(CPU_STAT)
88  startEvent = clock();
89 #endif
90 
91  for (i=-ghosts; i<(imax+ghosts); i++) {
92  for (j=-ghosts; j<(jmax+ghosts); j++) {
93  for (k=-ghosts; k<(kmax+ghosts); k++) {
94  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
95  double energy,pressure;
96  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
97  _NavierStokes3DGetFlowVar_( (u+p),
98  _NavierStokes3D_stride_,
99  Q[p+0],
100  Q[p+1],
101  Q[p+2],
102  Q[p+3],
103  energy,
104  pressure,
105  physics->gamma);
106  Q[p+4] = physics->gamma*pressure/Q[p+0]; /* temperature */
107  }
108  }
109  }
110 
111 #if defined(CPU_STAT)
112  stopEvent = clock();
113  printf("%-50s CPU time (secs) = %.6f\n",
114  "NaverStokes3DParabolicFunction_Q",(double)(stopEvent-startEvent)/CLOCKS_PER_SEC);
115 #endif
116 
117  double *QDerivX = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
118  double *QDerivY = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
119  double *QDerivZ = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
120 
121  IERR solver->FirstDerivativePar(QDerivX,Q,_XDIR_,1,solver,mpi); CHECKERR(ierr);
122  IERR solver->FirstDerivativePar(QDerivY,Q,_YDIR_,1,solver,mpi); CHECKERR(ierr);
123  IERR solver->FirstDerivativePar(QDerivZ,Q,_ZDIR_,1,solver,mpi); CHECKERR(ierr);
124 
125  IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,
126  solver->ghosts,mpi,QDerivX); CHECKERR(ierr);
127  IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,
128  solver->ghosts,mpi,QDerivY); CHECKERR(ierr);
129  IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,
130  solver->ghosts,mpi,QDerivY); CHECKERR(ierr);
131 
132  for (i=-ghosts; i<(imax+ghosts); i++) {
133  for (j=-ghosts; j<(jmax+ghosts); j++) {
134  for (k=-ghosts; k<(kmax+ghosts); k++) {
135  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
136  double dxinv, dyinv, dzinv;
137  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
138  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
139  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
140  _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv);
141  _ArrayScale1D_((QDerivX+p),dxinv,_MODEL_NVARS_);
142  _ArrayScale1D_((QDerivY+p),dyinv,_MODEL_NVARS_);
143  _ArrayScale1D_((QDerivZ+p),dzinv,_MODEL_NVARS_);
144  }
145  }
146  }
147 
148  double *FViscous = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
149  double *FDeriv = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
150 
151  /* Along X */
152  for (i=-ghosts; i<(imax+ghosts); i++) {
153  for (j=-ghosts; j<(jmax+ghosts); j++) {
154  for (k=-ghosts; k<(kmax+ghosts); k++) {
155  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
156  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
157 
158  double uvel, vvel, wvel, T, Tx,
159  ux, uy, uz, vx, vy, wx, wz;
160  uvel = (Q+p)[1];
161  vvel = (Q+p)[2];
162  wvel = (Q+p)[3];
163  T = (Q+p)[4];
164  Tx = (QDerivX+p)[4];
165  ux = (QDerivX+p)[1];
166  vx = (QDerivX+p)[2];
167  wx = (QDerivX+p)[3];
168  uy = (QDerivY+p)[1];
169  vy = (QDerivY+p)[2];
170  uz = (QDerivZ+p)[1];
171  wz = (QDerivZ+p)[3];
172 
173  /* calculate viscosity coeff based on Sutherland's law */
174  double mu = raiseto(T, 0.76);
175  /*
176  if (p/_MODEL_NVARS_ == 49841) {
177  printf("[CPU] 49841 mu = %.20e T = %.20e log(0.76) = %.20e T*log(0.76) = %.20e\n",
178  mu, T, log(0.76), T*log(0.76));
179  }
180  */
181 
182  double tau_xx, tau_xy, tau_xz, qx;
183  tau_xx = two_third * (mu*inv_Re) * (2*ux - vy - wz);
184  tau_xy = (mu*inv_Re) * (uy + vx);
185  tau_xz = (mu*inv_Re) * (uz + wx);
186  qx = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Tx;
187 
188  (FViscous+p)[0] = 0.0;
189  (FViscous+p)[1] = tau_xx;
190  (FViscous+p)[2] = tau_xy;
191  (FViscous+p)[3] = tau_xz;
192  (FViscous+p)[4] = uvel*tau_xx + vvel*tau_xy + wvel*tau_xz + qx;
193  }
194  }
195  }
196 
197  IERR solver->FirstDerivativePar(FDeriv,FViscous,_XDIR_,-1,solver,mpi); CHECKERR(ierr);
198  for (i=0; i<imax; i++) {
199  for (j=0; j<jmax; j++) {
200  for (k=0; k<kmax; k++) {
201  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
202  double dxinv;
203  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
204  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
205  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dxinv * (FDeriv+p)[v] );
206  }
207  }
208  }
209 
210  /* Along Y */
211  for (i=-ghosts; i<(imax+ghosts); i++) {
212  for (j=-ghosts; j<(jmax+ghosts); j++) {
213  for (k=-ghosts; k<(kmax+ghosts); k++) {
214  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
215  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
216 
217  double uvel, vvel, wvel, T, Ty,
218  ux, uy, vx, vy, vz, wy, wz;
219  uvel = (Q+p)[1];
220  vvel = (Q+p)[2];
221  wvel = (Q+p)[3];
222  T = (Q+p)[4];
223  Ty = (QDerivY+p)[4];
224  ux = (QDerivX+p)[1];
225  vx = (QDerivX+p)[2];
226  uy = (QDerivY+p)[1];
227  vy = (QDerivY+p)[2];
228  wy = (QDerivY+p)[3];
229  vz = (QDerivZ+p)[2];
230  wz = (QDerivZ+p)[3];
231 
232  /* calculate viscosity coeff based on Sutherland's law */
233  double mu = raiseto(T, 0.76);
234 
235  double tau_yx, tau_yy, tau_yz, qy;
236  tau_yx = (mu*inv_Re) * (uy + vx);
237  tau_yy = two_third * (mu*inv_Re) * (-ux + 2*vy - wz);
238  tau_yz = (mu*inv_Re) * (vz + wy);
239  qy = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Ty;
240 
241  (FViscous+p)[0] = 0.0;
242  (FViscous+p)[1] = tau_yx;
243  (FViscous+p)[2] = tau_yy;
244  (FViscous+p)[3] = tau_yz;
245  (FViscous+p)[4] = uvel*tau_yx + vvel*tau_yy + wvel*tau_yz + qy;
246  }
247  }
248  }
249 
250  IERR solver->FirstDerivativePar(FDeriv,FViscous,_YDIR_,-1,solver,mpi); CHECKERR(ierr);
251  for (i=0; i<imax; i++) {
252  for (j=0; j<jmax; j++) {
253  for (k=0; k<kmax; k++) {
254  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
255  double dyinv;
256  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
257  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
258  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dyinv * (FDeriv+p)[v] );
259  }
260  }
261  }
262 
263  /* Along Z */
264  for (i=-ghosts; i<(imax+ghosts); i++) {
265  for (j=-ghosts; j<(jmax+ghosts); j++) {
266  for (k=-ghosts; k<(kmax+ghosts); k++) {
267  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
268  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
269 
270  double uvel, vvel, wvel, T, Tz,
271  ux, uz, vy, vz, wx, wy, wz;
272  uvel = (Q+p)[1];
273  vvel = (Q+p)[2];
274  wvel = (Q+p)[3];
275  T = (Q+p)[4];
276  Tz = (QDerivZ+p)[4];
277  ux = (QDerivX+p)[1];
278  wx = (QDerivX+p)[3];
279  vy = (QDerivY+p)[2];
280  wy = (QDerivY+p)[3];
281  uz = (QDerivZ+p)[1];
282  vz = (QDerivZ+p)[2];
283  wz = (QDerivZ+p)[3];
284 
285  /* calculate viscosity coeff based on Sutherland's law */
286  double mu = raiseto(T,0.76);
287 
288  double tau_zx, tau_zy, tau_zz, qz;
289  tau_zx = (mu*inv_Re) * (uz + wx);
290  tau_zy = (mu*inv_Re) * (vz + wy);
291  tau_zz = two_third * (mu*inv_Re) * (-ux - vy + 2*wz);
292  qz = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Tz;
293 
294  (FViscous+p)[0] = 0.0;
295  (FViscous+p)[1] = tau_zx;
296  (FViscous+p)[2] = tau_zy;
297  (FViscous+p)[3] = tau_zz;
298  (FViscous+p)[4] = uvel*tau_zx + vvel*tau_zy + wvel*tau_zz + qz;
299  }
300  }
301  }
302 
303  IERR solver->FirstDerivativePar(FDeriv,FViscous,_ZDIR_,-1,solver,mpi); CHECKERR(ierr);
304  for (i=0; i<imax; i++) {
305  for (j=0; j<jmax; j++) {
306  for (k=0; k<kmax; k++) {
307  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
308  double dzinv;
309  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
310  _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv);
311  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dzinv * (FDeriv+p)[v] );
312  }
313  }
314  }
315 
316  free(Q);
317  free(QDerivX);
318  free(QDerivY);
319  free(QDerivZ);
320  free(FViscous);
321  free(FDeriv);
322 
323  if (solver->flag_ib) _ArrayBlockMultiply_(par,solver->iblank,size,_MODEL_NVARS_);
324  return(0);
325 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_NavierStokes3DGetFlowVar_
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
Definition: navierstokes3d.h:98
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
raiseto
#define raiseto(x, a)
Definition: math_ops.h:37
HyPar::iblank
double * iblank
Definition: hypar.h:436
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::flag_ib
int flag_ib
Definition: hypar.h:441
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_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
_ZDIR_
#define _ZDIR_
Definition: navierstokes3d.h:76
HyPar::count_par
int count_par
Definition: hypar.h:418
_ArrayBlockMultiply_
#define _ArrayBlockMultiply_(x, a, n, bs)
Definition: arrayfunctions.h:388
HyPar::FirstDerivativePar
int(* FirstDerivativePar)(double *, double *, int, int, void *, void *)
Definition: hypar.h:243
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
HyPar::dxinv
double * dxinv
Definition: hypar.h:110
NavierStokes3D
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
Definition: navierstokes3d.h:482
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
_NavierStokes3D_stride_
static const int _NavierStokes3D_stride_
Definition: navierstokes3d.h:559