Calculates the stress tensor components (at the i+1/2 face of each cell) and adds explicit viscous terms to the energy and momentum equation. It is called in the during the sweep integrators. The stress tensor is given by
73 double dx_1, dy_1, dz_1;
76 double *x1r, *x2r, *x3r;
77 double dxVx, dxVy, dxVz, dyVx, dyVy, dyVz, dzVx, dzVy, dzVz;
78 static double *tau_xx, *tau_xy, *tau_xz,
79 *tau_yx, *tau_yy, *tau_yz,
80 *tau_zx, *tau_zy, *tau_zz;
81 double ***Vx, ***Vy, ***Vz;
83 double *r, *th,
r_1, dr, dr_1, s_1, tan_1, dy;
84 double dVxi,dVyi,dVzi;
85 double dVxj,dVyj,dVzj;
86 double dVxk,dVyk,dVzk;
88 static double *one_dVr, *one_dmu;
106 #if GEOMETRY != CARTESIAN
108 if (one_dVr == NULL){
111 for (i = 0; i <
NX1_TOT - 1; i++){
112 one_dVr[
i] = r[i+1]*fabs(r[i + 1]) - r[
i]*fabs(r[i]);
113 one_dVr[
i] = 2.0/one_dVr[
i];
115 for (j = 0; j <
NX2_TOT - 1; j++){
116 one_dmu[
j] = 1.0 - cos(th[j + 1]) - (1.0 - cos(th[j]))*(th[j] > 0.0 ? 1.0:-1.0);
117 one_dmu[
j] = 1.0/one_dmu[
j];
136 for (i =
beg; i <=
end; i++){
142 vi[nv] = 0.5*(V[nv][
k][
j][
i] + V[nv][
k][
j][i+1]);
143 vc[nv] = V[nv][
k][
j][
i];
148 Visc_nu(vi, x1r[i], x2[j], x3[k], &nu1, &nu2);
150 tau_xx[
i] = tau_xy[
i]= tau_xz[
i]= tau_yx[
i]=
151 tau_yy[
i] = tau_yz[
i]= tau_zx[
i]= tau_zy[
i]=
154 dVxi = dVxj = dVxk = dVyi = dVyj = dVyk = dVzi = dVzj = dVzk = 0.0;
155 dxVx = dxVy = dxVz = dyVx = dyVy = dyVz = dzVx = dzVy = dzVz = 0.0;
165 EXPAND(dxVx = dVxi*dx_1; , dxVy = dVyi*dx_1; , dxVz = dVzi*dx_1; )
167 EXPAND(dyVx = dVxj*dy_1; , dyVy = dVyj*dy_1; , dyVz = dVzj*dy_1; )
169 dzVx = dVxk*dz_1; dzVy = dVyk*dz_1; dzVz = dVzk*dz_1;
173 #if GEOMETRY == CARTESIAN
174 div_v =
D_EXPAND(dxVx, + dyVy , + dzVz);
178 tau_xx[
i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v;
179 tau_xy[
i] = nu1*(dyVx + dxVy);
180 tau_xz[
i] = nu1*(dzVx + dxVz);
181 tau_yx[
i] = tau_xy[
i];
182 tau_zx[
i] = tau_xz[
i];
186 EXPAND(ViS[i][
MX1] = 0.0; ,
193 dr_1 = 1.0/dr; dz_1 = 0.0;
194 r_1 = 1.0/grid[
IDIR].
x[
i];
198 dxVx = (Vx[
k][
j][i+1]*r[i+1] - Vx[
k][
j][
i]*fabs(r[i]))*one_dVr[i];
199 div_v =
D_EXPAND( dxVx, + dVyj*dy_1, + 0.0);
204 tau_xx[
i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v;
205 tau_xy[
i] = nu1*(dyVx + dxVy);
206 EXPAND(tau_xz[i] = nu1*(0.0);,
207 tau_xz[
i] = nu1*(0.0);,
208 tau_xz[
i] = nu1*0.5*(r[
i]+r[i+1])*dr_1*((1./r[i+1])*Vz[
k][
j][i+1]
209 - (1./r[
i])*Vz[k][j][i]);)
211 tau_yx[i] = tau_xy[i];
212 tau_zx[
i] = tau_xz[
i];
217 Visc_nu(vc, x1[i], x2[j], x3[k], &nu1, &nu2);
219 div_v =
D_EXPAND( 0.5*(Vx[k][j][i+1]-Vx[k][j][i-1])*dx_1 + Vx[k][j][i]*r_1,
220 + 0.5*(Vy[k][j + 1][i]-Vy[k][j - 1][i])*dy_1,
223 tau_zz[
i] = 2.0*nu1*r_1*Vx[
k][
j][
i]
224 + (nu2 - (2.0/3.0)*nu1)*div_v;
228 EXPAND(ViS[i][
MX1] = -tau_zz[i]*r_1; ,
235 dr = grid[
IDIR].
dx[
i]; dr_1 = 1.0/dr;
240 div_v =
D_EXPAND(r_1*vi[
VX1] + dVxi*dr_1, + r_1*dVyj*dy_1, + dVzk*dz_1);
244 tau_xx[
i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v;
245 EXPAND(tau_xy[i] = nu1*dxVy; ,
246 tau_xy[
i] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]); ,
247 tau_xy[
i] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]);)
249 tau_xz[i] = nu1*(dzVx + dxVz);
250 tau_yx[
i] = tau_xy[
i];
251 EXPAND(tau_yy[i] = 2.0*nu1*r_1*dyVy + (nu2 - (2.0/3.0)*nu1)*div_v;,
252 tau_yy[
i] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
253 + (nu2 - (2.0/3.0)*nu1)*div_v;,
254 tau_yy[
i] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
255 + (nu2 - (2.0/3.0)*nu1)*div_v;)
257 tau_zx[i] = tau_xz[i];
261 r_1 = 1.0/grid[
IDIR].
x[
i];
263 EXPAND(ViS[i][
MX1] = -0.5*(tau_yy[i - 1] + tau_yy[i])*r_1; ,
271 tan_1= 1.0/tan(th[j]); s_1 = 1.0/sin(th[j]);
276 + r_1*dVyj*dy_1 + r_1*tan_1*vi[
VX2],
277 + r_1*s_1*dVzk*dz_1);
282 tau_xx[
i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v;
283 EXPAND(tau_xy[i] = nu1*(r_1*dyVx + dxVy ); ,
284 tau_xy[
i] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]); ,
285 tau_xy[
i] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]);)
287 EXPAND(tau_xz[i] = nu1*(r_1*s_1*dzVx + dxVz) ; ,
288 tau_xz[
i] = nu1*(r_1*s_1*dzVx + dxVz ); ,
289 tau_xz[
i] = nu1*(r_1*s_1*dzVx + dxVz - r_1*vi[
VX3]);)
291 tau_yx[i] = tau_xy[i];
292 tau_yy[
i] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
293 + (nu2 - (2.0/3.0)*nu1)*div_v;
295 tau_zx[
i] = tau_xz[
i];
296 EXPAND(tau_zz[i] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[
VX1] )
297 + (nu2 - (2.0/3.0)*nu1)*div_v;,
298 tau_zz[
i] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[
VX1]
300 + (nu2 - (2.0/3.0)*nu1)*div_v;,
301 tau_zz[
i] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[
VX1]
303 + (nu2 - (2.0/3.0)*nu1)*div_v;)
308 r_1 = 1.0/grid[
IDIR].
x[i];
309 EXPAND(ViS[i][
MX1] = -0.5*( tau_yy[i - 1] + tau_yy[i]
310 + tau_zz[i - 1] + tau_zz[i])*r_1;,
318 EXPAND(ViF[i][
MX1] = tau_xx[i]; ,
319 ViF[
i][
MX2] = tau_xy[
i]; ,
320 ViF[
i][
MX3] = tau_xz[
i]; )
323 ViF[i][ENG] = EXPAND( vi[VX1]*tau_xx[i],
325 + vi[
VX3]*tau_zx[i]);
328 dcoeff[
i][
MX1] =
MAX(nu1, nu2);
340 for (j =
beg ; j <=
end; j++){
346 vi[nv] = 0.5*(V[nv][
k][j+1][
i] + V[nv][
k][
j][
i]);
347 vc[nv] = V[nv][
k][
j][
i];
352 Visc_nu(vi, x1[i], x2r[j], x3[k], &nu1, &nu2);
354 tau_xx[
j]= tau_xy[
j]= tau_xz[
j]= tau_yx[
j] = tau_yy[
j]=
355 tau_yz[
j]= tau_zx[
j]= tau_zy[
j]= tau_zz[
j] =0.0;
357 dVxi = dVxj = dVxk = dVyi = dVyj = dVyk = dVzi = dVzj = dVzk = 0.0;
358 dxVx = dxVy = dxVz = dyVx = dyVy = dyVz = dzVx = dzVy = dzVz = 0.0;
368 EXPAND(dxVx = dVxi*dx_1; , dxVy = dVyi*dx_1; , dxVz = dVzi*dx_1; )
369 EXPAND(dyVx = dVxj*dy_1; , dyVy = dVyj*dy_1; , dyVz = dVzj*dy_1; )
371 dzVx = dVxk*dz_1; dzVy = dVyk*dz_1; dzVz = dVzk*dz_1;
374 #if GEOMETRY == CARTESIAN
378 div_v =
D_EXPAND(dxVx, + dyVy, + dzVz);
382 tau_xy[
j] = nu1*(dyVx + dxVy);
383 tau_yx[
j] = tau_xy[
j];
384 tau_yy[
j] = 2.0*nu1*dyVy + (nu2 - (2.0/3.0)*nu1)*div_v;
385 tau_yz[
j] = nu1*(dyVz + dzVy);
386 tau_zy[
j] = tau_yz[
j];
390 EXPAND(ViS[j][
MX1] = 0.0; ,
397 dr = grid[
IDIR].
dx[
i]; dr_1 = 1.0/dr;
398 r_1 = 1.0/grid[
IDIR].
x[
i]; dz_1 = 0.0;
402 div_v =
D_EXPAND(r_1*vi[VX1] + dVxi*dr_1, + dVyj*dy_1, + 0.0);
406 tau_xy[
j] = nu1*(dyVx + dxVy);
407 tau_yx[
j] = tau_xy[
j];
408 tau_yy[
j] = 2.0*nu1*dyVy + (nu2 - (2.0/3.0)*nu1)*div_v;
409 EXPAND(tau_yz[j] = 0.0; ,
411 tau_yz[
j] = nu1*(dyVz);)
412 tau_zy[j] = tau_yz[j];
416 EXPAND(ViS[j][
MX1] = 0.0; ,
423 dr_1 = 1.0/dr; r_1 = 1.0/grid[
IDIR].
x[
i];
427 div_v =
D_EXPAND(r_1*vi[VX1] + dVxi*dr_1, + r_1*dVyj*dy_1, + dVzk*dz_1);
431 tau_xy[
j] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]);
433 tau_yx[
j] = tau_xy[
j];
434 tau_yy[
j] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
435 + (nu2 - (2.0/3.0)*nu1)*div_v;
437 tau_yz[
j] = nu1*(r_1*dyVz + dzVy);
438 tau_zy[
j] = tau_yz[
j];
440 r_1 = 1.0/grid[
IDIR].
x[
i];
444 EXPAND(ViS[j][
MX1] = 0.0; ,
452 tan_1= 1.0/tan(th[j]); s_1 = 1.0/sin(th[j]);
465 if (fabs(tan(th[j]))< 1.e-12) tan_1 = 0.0;
471 div_v =
D_EXPAND( 2.0*r_1*vi[VX1] + dVxi*dr_1 ,
472 + ( sin(th[j + 1])*Vy[k][j + 1][i] - fabs(sin(th[j]))*Vy[k][j][i])*r_1*one_dmu[j],
473 + r_1*s_1*dVzk*dz_1 );
479 tau_xy[
j] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]);
481 tau_yx[
j] = tau_xy[
j];
482 tau_yy[
j] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
483 + (nu2 - (2.0/3.0)*nu1)*div_v;
486 EXPAND(tau_yz[j] = nu1*r_1*dyVz; ,
487 tau_yz[
j] = nu1*r_1*dyVz; ,
488 tau_yz[
j] = nu1*(r_1*dyVz - tan_1*r_1*vi[
VX3]));
491 tau_yz[
j] = nu1*(s_1*r_1*dzVy + r_1*dyVz - tan_1*r_1*vi[
VX3]);
493 tau_zy[
j] = tau_yz[
j];
495 EXPAND(tau_zz[j] = 2.0*nu1*r_1*vi[VX1]
496 + (nu2 - (2.0/3.0)*nu1)*div_v; ,
497 tau_zz[
j] = 2.0*nu1*(r_1*vi[
VX1] + tan_1*r_1*vi[
VX2])
498 + (nu2 - (2.0/3.0)*nu1)*div_v;,
499 tau_zz[
j] = 2.0*nu1*(r_1*vi[
VX1] + tan_1*r_1*vi[
VX2])
500 + (nu2 - (2.0/3.0)*nu1)*div_v;)
503 tau_zz[j] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[VX1] + tan_1*r_1*vi[
VX2])
504 + (nu2 - (2.0/3.0)*nu1)*div_v;
509 th = grid[
JDIR].
x; tan_1= 1.0/tan(th[j]);
511 EXPAND(ViS[j][
MX1] = 0.0; ,
512 ViS[
j][
MX2] = 0.5*( tau_yx[j - 1] + tau_yx[
j])*r_1
513 - tan_1*0.5*(tau_zz[j - 1] + tau_zz[j])*
r_1; ,
520 EXPAND(ViF[j][
MX1] = tau_yx[j]; ,
521 ViF[
j][
MX2] = tau_yy[
j]; ,
522 ViF[
j][
MX3] = tau_yz[
j]; )
525 ViF[j][ENG] = EXPAND( vi[VX1]*tau_xy[j],
527 + vi[
VX3]*tau_zy[j]);
542 for (k =
beg; k <=
end; k++){
548 vi[nv] = 0.5*(V[nv][
k][
j][
i] + V[nv][k+1][
j][
i]);
549 vc[nv] = V[nv][
k][
j][
i];
554 Visc_nu(vi, x1[i], x2[j], x3r[k], &nu1, &nu2);
556 tau_xx[
k] = tau_xy[
k] = tau_xz[
k]= tau_yx[
k]= tau_yy[
k]= tau_yz[
k]=
557 tau_zx[
k] = tau_zy[
k] = tau_zz[
k]=0.0;
558 dVxi = dVxj = dVxk = dVyi = dVyj = dVyk = dVzi = dVzj = dVzk = 0.0;
566 dxVx = dVxi*dx_1; dxVy = dVyi*dx_1; dxVz = dVzi*dx_1;
567 dyVx = dVxj*dy_1; dyVy = dVyj*dy_1; dyVz = dVzj*dy_1;
568 dzVx = dVxk*dz_1; dzVy = dVyk*dz_1; dzVz = dVzk*dz_1;
570 #if GEOMETRY == CARTESIAN
574 div_v = dxVx + dyVy + dzVz;
578 tau_xz[
k] = nu1*(dzVx + dxVz);
579 tau_yz[
k] = nu1*(dyVz + dzVy);
580 tau_zx[
k] = tau_xz[
k];
581 tau_zy[
k] = tau_yz[
k];
582 tau_zz[
k] = 2.0*nu1*dzVz + (nu2 - (2.0/3.0)*nu1)*div_v;
586 EXPAND(ViS[k][
MX1] = 0.0; ,
593 dr = grid[
IDIR].
dx[
i]; dr_1 = 1.0/dr;
594 r_1 = 1.0/grid[
IDIR].
x[
i];
598 div_v = r_1*vi[
VX1] + dVxi*dr_1 + r_1*dVyj*dy_1 + dVzk*dz_1;
602 tau_xz[
k] = nu1*(dzVx + dxVz);
603 tau_yy[
k] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
604 + (nu2 - (2.0/3.0)*nu1)*div_v;
606 tau_yz[
k] = nu1*(r_1*dyVz + dzVy);
607 tau_zx[
k] = tau_xz[
k];
608 tau_zy[
k] = tau_yz[
k];
609 tau_zz[
k] = 2.0*nu1*dzVz + (nu2 - (2.0/3.0)*nu1)*div_v;
612 EXPAND(ViS[k][
MX1] = 0.0; ,
619 r_1 = 1.0/grid[
IDIR].
x[
i]; tan_1= 1.0/tan(th[j]); s_1 = 1.0/sin(th[j]);
623 div_v = 2.0*r_1*vi[
VX1] + dVxi*dr_1 + r_1*dVyj*dy_1 + r_1*tan_1*vi[
VX2]
628 tau_xz[
k] = nu1*(r_1*s_1*dzVx + dxVz - r_1*vi[
VX3]);
630 tau_yz[
k] = nu1*(s_1*r_1*dzVy + r_1*dyVz - tan_1*r_1*vi[
VX3]);
632 tau_zx[
k] = tau_xz[
k];
633 tau_zy[
k] = tau_yz[
k];
634 tau_zz[
k] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[
VX1]
636 + (nu2 - (2.0/3.0)*nu1)*div_v;
641 EXPAND(ViS[k][
MX1] = 0.0; ,
649 ViF[k][
MX1] = tau_zx[k];
650 ViF[
k][
MX2] = tau_zy[
k];
651 ViF[
k][
MX3] = tau_zz[
k];
653 #if EOS != ISOTHERMAL
654 ViF[
k][ENG] = vi[
VX1]*tau_xz[
k]
658 dcoeff[
k][
MX1] =
MAX(nu1, nu2);
int end
Global end index for the local array.
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
int g_dir
Specifies the current sweep or direction of integration.
int beg
Global start index for the local array.
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
double * r_1
Geometrical factor 1/r.
#define ARRAY_1D(nx, type)
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
void Visc_nu(double *v, double x1, double x2, double x3, double *nu1, double *nu2)
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.