1 subroutine e3source(xx, src) 2c----------------------------------------------------------------------- 3c 4c this routine computes the body force term. 5c 6c currently this computes a swirl body with the axis alligned with 7c the z-coordinate 8c 9c----------------------------------------------------------------------- 10 11 include "common.h" 12 13 real*8 xx(npro,nsd), src(npro,nsd) 14 15 real*8 nu 16 17 real*8 r, Stheta, dpdz, rP5 18 if(datmat(2,5,1).eq. 0) then ! calc swirl 19c 20c This is the body force which will drive a swirl in a pipe flow 21c 22 bigR = 0.5d0 23 dpdz = datmat(1,5,1) 24 do iel = 1, npro 25 26 r = sqrt( xx(iel,1)**2 + xx(iel,2)**2) 27 rP5 = (r/bigR)**5 28 29 Stheta = dpdz * sin(0.5*pi*rP5) 30 31 src(iel,1) = -xx(iel,2)/r * Stheta 32 src(iel,2) = xx(iel,1)/r * Stheta 33 src(iel,3) = dpdz 34 enddo 35 else ! contrived test problem 36 37c$$$ nu=datmat(1,2,1) 38c$$$ omag=datmat(3,5,1) ! frame rotation rate 39c$$$ e1 = one/sqrt(3.0d0) ! axis of rotation 40c$$$ e2 = e1 41c$$$ e3 = e1 42c$$$ om1=omag*e1 43c$$$ om2=omag*e2 44c$$$ om3=omag*e3 45c$$$ w=0 46c$$$ 47c$$$ do iel = 1, npro 48c$$$ 49c$$$ x = xx(iel,1) 50c$$$ y = xx(iel,2) 51c$$$ z = xx(iel,3) 52c$$$ 53c$$$c 54c$$$c The following are MAPLE generated forcing functions for 55c$$$c a contrived flow in a rotating reference frame. 56c$$$c 57c$$$ 58c$$$ t1 = x**2 59c$$$ t2 = t1**2 60c$$$ t5 = dexp(14.D0*x) 61c$$$ t6 = t2*x*t5 62c$$$ t7 = y**2 63c$$$ t8 = t7**2 64c$$$ t9 = t8*t7 65c$$$ t12 = t2*t5 66c$$$ t15 = nu*t1 67c$$$ t17 = dexp(7.D0*x) 68c$$$ t18 = t7*y 69c$$$ t19 = t17*t18 70c$$$ t22 = t1*x 71c$$$ t23 = nu*t22 72c$$$ t24 = t17*y 73c$$$ t29 = t17*t7 74c$$$ t34 = nu*x 75c$$$ t43 = nu*t2 76c$$$ t46 = -40.D0*t6*t9-6.D0*t12*t7+92.D0*t15*t19+132.D0*t23*t24+80.D0* 77c$$$ #t6*t18-252.D0*t23*t29+168.D0*t23*t19+96.D0*t34*t29-64.D0*t34*t19+2 78c$$$ #2.D0*t15*t24-138.D0*t15*t29+294.D0*t43*t29 79c$$$ t50 = t2*t22*t5 80c$$$ t57 = nu*t17 81c$$$ t60 = t8*y 82c$$$ t65 = t2**2 83c$$$ t66 = t65*t5 84c$$$ t73 = t22*t5 85c$$$ t77 = t2*t1*t5 86c$$$ t80 = -196.D0*t43*t19-96.D0*t50*t9-122.D0*t43*t24-32.D0*t34*t24-4. 87c$$$ #D0*t57*y+36.D0*t12*t60+192.D0*t50*t18+98.D0*t66*t8+24.D0*t12*t18+2 88c$$$ #8.D0*t66*t9-24.D0*t73*t60-336.D0*t77*t60 89c$$$ t106 = -336.D0*t50*t8-12.D0*t12*t9+8.D0*t73*t9-140.D0*t6*t8+28.D0* 90c$$$ #t73*t8+12.D0*t15*t17+12.D0*t43*t17+288.D0*t50*t60+112.D0*t77*t9-22 91c$$$ #4.D0*t77*t18-56.D0*t66*t18+120.D0*t6*t60 92c$$$ t131 = -8.D0*t57*t18-20.D0*t6*t7+56.D0*t77*t7-42.D0*t12*t8-24.D0*t 93c$$$ #23*t17-48.D0*t50*t7+4.D0*t73*t7-16.D0*t73*t18+392.D0*t77*t8+14.D0* 94c$$$ #t66*t7-84.D0*t66*t60+12.D0*t57*t7 95c$$$ fx = t46+t80+t106+t131 96c$$$ 97c$$$ 98c$$$ t1 = x**2 99c$$$ t2 = t1**2 100c$$$ t3 = nu*t2 101c$$$ t5 = dexp(7.D0*x) 102c$$$ t6 = t5*y 103c$$$ t9 = y**2 104c$$$ t10 = t9**2 105c$$$ t11 = nu*t10 106c$$$ t18 = t1*x 107c$$$ t22 = nu*x 108c$$$ t25 = nu*t18 109c$$$ t32 = dexp(14.D0*x) 110c$$$ t33 = t2*t32 111c$$$ t39 = t2*t1*t32 112c$$$ t40 = t10*t9 113c$$$ t43 = t1*t32 114c$$$ t46 = t9*y 115c$$$ t47 = t10*t46 116c$$$ t50 = -84.D0*t3*t6+343.D0*t11*t5*t2+66.D0*t11*t5*x-98.D0*t11*t5*t1 117c$$$ #8-24.D0*t22*t6+120.D0*t25*t6-287.D0*t11*t5*t1-140.D0*t33*t10+14.D0 118c$$$ #*t3*t5-56.D0*t39*t40-28.D0*t43*t40+8.D0*t43*t47 119c$$$ t52 = t2*x*t32 120c$$$ t55 = t18*t32 121c$$$ t62 = t10*y 122c$$$ t69 = nu*t5 123c$$$ t80 = 120.D0*t52*t10-32.D0*t55*t47+56.D0*t33*t47+30.D0*t11*t5-216. 124c$$$ #D0*t52*t62-144.D0*t55*t62+72.D0*t39*t62-60.D0*t69*t46+30.D0*t69*t9 125c$$$ #-20.D0*t25*t5-20.D0*t43*t10+36.D0*t43*t62 126c$$$ t86 = nu*t1 127c$$$ t89 = t5*t9 128c$$$ t92 = t5*t46 129c$$$ t109 = 112.D0*t55*t40+80.D0*t55*t10-12.D0*t86*t6-275.D0*t86*t89+57 130c$$$ #4.D0*t86*t92-218.D0*t25*t89+196.D0*t25*t92+90.D0*t22*t89-132.D0*t2 131c$$$ #2*t92+427.D0*t3*t89+8.D0*t39*t46+4.D0*t43*t46 132c$$$ t134 = 16.D0*t39*t47+28.D0*t33*t46-48.D0*t52*t47+252.D0*t33*t62+16 133c$$$ #8.D0*t52*t40-24.D0*t52*t46-686.D0*t3*t92+2.D0*t86*t5-40.D0*t39*t10 134c$$$ #-196.D0*t33*t40-16.D0*t55*t46+4.D0*t22*t5 135c$$$ fy = t50+t80+t109+t134 136c$$$ 137c$$$ 138c$$$ t3 = om3*x 139c$$$ t5 = dexp(7.D0*x) 140c$$$ t7 = x**2 141c$$$ t12 = y**2 142c$$$ t15 = (-1.D0+y)**2 143c$$$ fxa = 2.D0*om2*w+2.D0*t3*t5*(2.D0+x-10.D0*t7+7.D0*t7*x)*t12*t15+om 144c$$$ #2*(om1*y-om2*x)-om3*(t3-om1*z) 145c$$$ 146c$$$ t1 = x**2 147c$$$ t4 = (-1.D0+x)**2 148c$$$ t7 = dexp(7.D0*x) 149c$$$ t10 = y**2 150c$$$ fya = 4.D0*om3*t1*t4*t7*y*(1.D0-3.D0*y+2.D0*t10)-2.D0*om1*w+om3*(o 151c$$$ #m2*z-om3*y)-om1*(om1*y-om2*x) 152c$$$ 153c$$$ 154c$$$ t3 = dexp(7.D0*x) 155c$$$ t5 = x**2 156c$$$ t10 = y**2 157c$$$ t13 = (-1.D0+y)**2 158c$$$ t19 = (-1.D0+x)**2 159c$$$ fza = -2.D0*om1*x*t3*(2.D0+x-10.D0*t5+7.D0*t5*x)*t10*t13-4.D0*om2* 160c$$$ #t5*t19*t3*y*(1.D0-3.D0*y+2.D0*t10)+om1*(om3*x-om1*z)-om2*(om2*z-om 161c$$$ #3*y) 162c$$$ 163c$$$ src(iel,1) = fx + fxa 164c$$$ src(iel,2) = fy + fya 165c$$$ src(iel,3) = fza 166c$$$ enddo 167!Analytic lid case 168 169 do iel = 1, npro 170 x = xx(iel,1) 171 y = xx(iel,2) 172 z = xx(iel,3) 173 Re = 1.0/datmat(1,2,1) 174c 175c The following are MAPLE generated forcing functions for 176c a Lid Driven cavity flow with an analytic solution 177c 178 t2 = x**2 179 t3 = t2**2 180 t4 = t3*x 181 t7 = t2*x 182 t8 = 8*t7 183 t13 = y**2 184 t20 = t13**2 185 t21 = t20-t13 186 t26 = t3**2 187 t30 = t3*t2 188 t37 = t13*y 189 t54 = -8/Re*(24.E0/5.E0*t4-12*t3+t8+2*(4*t7-6*t2+2*x)*(12 190 & *t13-2)+(24*x-12)*t21)-64*(t26/2-2*t3*t7+3*t30-2*t4+t3 191 & /2)*(-24*t20*y+8*t37-4*y)+64*t21*(4*t37-2*y)*(-4*t30+12 192 & *t4-14*t3+t8-2*t2) 193 194 src(iel,1) = 0.0 195 src(iel,2) = t54 196 src(iel,3) = 0.0 197 enddo 198 199 endif 200 return 201 end 202 203 204!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 205!****************************************************************** 206!------------------------------------------------------------------- 207 208 209 210 subroutine e3sourceSclr ( Sclr, Sdot, gradS, 211 & dwl, shape_funct, shg, 212 & yl, dxidx, rmu, 213 & u1, u2, u3, xl, 214 & srcR, srcL, uMod, 215 & srcRat) 216 217 218c----------------------------------------------------------------------- 219c 220c calculate the coefficients for the Spalart-Allmaras 221c turbulence model and its jacobian 222c 223c input: 224c Sclr : The scalar that is being transported in this pde. 225c gradS : spatial gradient of Sclr 226c rmu : kinematic viscosity 227c 228c output: 229c rmu : diffusion for eddy viscosity equation 230c srcR : residual terms for 231c (srcR * turb) 232c srcL : jacobian of srcR 233c 234c----------------------------------------------------------------------- 235 use turbSA 236 use turbKE 237 include "common.h" 238c coming in 239 real*8 Sclr(npro), Sdot(npro), 240 & gradS(npro,nsd), dwl(npro,nenl), 241 & shape_funct(npro,nshl), shg(npro,nshl,nsd), 242 & yl(npro,nenl,ndof), dxidx(npro,nsd,nsd), 243 & rmu(npro), u1(npro), 244 & u2(npro), u3(npro), 245 & xl(npro,nenl,nsd) 246c going out 247 real*8 srcR(npro), srcL(npro), 248 & uMod(npro,nsd) 249c used locally 250 251 real*8 gradv(npro,nsd,nsd), absVort(npro), 252 & dwall(npro) 253 real*8 chi, chiP3, fv1, fv2, st, r, 254 & g, fw, s, viscInv, k2d2Inv, 255 & dP2Inv, sixth, tmp, tmp1, p, dp, 256 & fv1p, fv2p, stp, gp, fwp, rp, 257 & chiP2, mytmp(npro), sign_levelset(npro), 258 & sclr_ls(npro) 259c Kay-Epsilon 260 real*8 Ry(npro), Rt(npro), RtP2(npro), f2(npro), f1(npro), 261 & kay(npro), epsilon(npro), fmu(npro), fmui(npro), 262 & srcjac(npro,4), srcRat(npro) 263 integer e 264c---------------------------------------------------------------------- 265c 266c -- -- -- -- 267c | cb2 | 1 | | 268c phi, + |u - -----phi, | phi, - -----|(nu+phi)phi, | 269c t | i sigma i| i sigma| i| 270c -- -- -- --, 271c i 272c 273c ~ phi^2 274c - cb1*S*phi + cw1*fw*----- 275c d^2 276c 277c---------------------------------------------------------------------- 278c 279 if(iRANS.eq.-1) then ! spalart almaras model 280c 281c.... compute the global gradient of u 282c 283 gradV = zero 284 do n = 1, nshl 285c 286c du_i/dx_j 287c 288c i j indices match array where V is the velocity (u in our notes) 289 gradV(:,1,1) = gradV(:,1,1) + shg(:,n,1) * yl(:,n,2) 290 gradV(:,2,1) = gradV(:,2,1) + shg(:,n,1) * yl(:,n,3) 291 gradV(:,3,1) = gradV(:,3,1) + shg(:,n,1) * yl(:,n,4) 292c 293 gradV(:,1,2) = gradV(:,1,2) + shg(:,n,2) * yl(:,n,2) 294 gradV(:,2,2) = gradV(:,2,2) + shg(:,n,2) * yl(:,n,3) 295 gradV(:,3,2) = gradV(:,3,2) + shg(:,n,2) * yl(:,n,4) 296c 297 gradV(:,1,3) = gradV(:,1,3) + shg(:,n,3) * yl(:,n,2) 298 gradV(:,2,3) = gradV(:,2,3) + shg(:,n,3) * yl(:,n,3) 299 gradV(:,3,3) = gradV(:,3,3) + shg(:,n,3) * yl(:,n,4) 300c a j u a i 301c from our notes where we had N_{a,j} = dN_a/dx_j note that i is off by one because p was first in yl vector 302c 303 enddo 304c 305c.... magnitude of vorticity 306c 307 absVort = sqrt( (gradV(:,2,3) - gradV(:,3,2)) ** 2 308 & + (gradV(:,3,1) - gradV(:,1,3)) ** 2 309 & + (gradV(:,1,2) - gradV(:,2,1)) ** 2 ) 310 dwall = zero 311 do n = 1, nenl 312 dwall = dwall + shape_funct(:,n) * dwl(:,n) 313 enddo 314 315 sixth = 1.0/6.0 316c 317c.... compute source and its jacobian 318c 319 do e = 1, npro 320 SclrNN= max(Sclr(e),zero) 321 322 if(iles.lt.0) then ! for DES 323 dx=maxval(xl(e,:,1))-minval(xl(e,:,1)) 324 dy=maxval(xl(e,:,2))-minval(xl(e,:,2)) 325 dz=maxval(xl(e,:,3))-minval(xl(e,:,3)) 326 dmax=max(dx,max(dy,dz)) 327 dmax=0.65*dmax 328c 329c Next 5 lines are DDES 330c 331 qfac = sqrt(gradV(e,1,1)**2+gradV(e,1,2)**2+gradV(e,1,3)**2 332 & +gradV(e,2,1)**2+gradV(e,2,2)**2+gradV(e,2,3)**2 333 & +gradV(e,3,1)**2+gradV(e,3,2)**2+gradV(e,3,3)**2) 334 rd=SclrNN*saKappaP2Inv/(dwall(e)**2*qfact) 335 fd=one-tanh((8.0000000000000000d0*rd)**3) 336 dwall(e)=dwall(e)-fd*max(zero,dwall(e)-dmax) 337c DES97 dwall(e)=min(dmax,dwall(e)) 338 endif 339 340 if (dwall(e) .ne. 0) dP2Inv = 1.0 / dwall(e)**2 341 342 k2d2Inv = dP2Inv * saKappaP2Inv 343 344 viscInv = 1.0/datmat(1,2,1) 345 chi = SclrNN * viscInv !skipping the kiy kie for now 346 chiP2 = chi * chi 347 chiP3 = chi * chiP2 348 349 ep=zero 350 if(Sclr(e).gt.zero) ep=one 351 352 tmp = 1.0 / ( chiP3 + saCv1P3 ) 353 fv1 = chiP3 * tmp 354 fv1p = 3.0 * viscInv * saCv1P3 * chiP2 * tmp**2 355 356 tmp = 1.0 / (1.0 + chi * fv1) 357 fv2 = 1.0 - chi * tmp 358 fv2p = (chiP2 * fv1p - viscInv) * tmp**2 359 360 s = absVort(e) !prd(e,2) 361 tmp = SclrNN * k2d2Inv !eOkdP2 362 st = s + fv2 * tmp 363 stp = ep*k2d2Inv * fv2 + tmp*fv2p 364 365 r = zero 366 rp = zero 367 if (st .gt. epsM ) then 368 r = tmp / st 369 r = min( max(r, -8.0d0), 8.0d0) 370 rp = k2d2Inv / st**2 * (ep*st - SclrNN*stp) 371 endif 372 rP5 = r * (r * r) ** 2 373 tmp = 1.0 + saCw2 * (rP5 - 1.0) 374 g = r * tmp 375 gp = rp * (tmp + 5.0 * saCw2 * rP5) 376 377 gP6 = (g * g * g) ** 2 378 tmp = 1.0 / (gP6 + saCw3P6) 379 tmp1 = ( (1.0 + saCw3P6) * tmp ) ** sixth 380 fw = g * tmp1 381 fwp = gp * tmp1 * saCw3P6 * tmp 382 if(abs(r).gt.7.0d0) fwp=zero 383 384 tmp1 = saCw1 * dP2Inv*SclrNN 385 prat = ep*saCb1*st !prodRat 386 srcRat(e)= prat - ep*fw * tmp1 387 388 tmp = zero 389 if(prat .lt. zero) tmp=prat 390 if(stp .lt. zero) tmp=tmp+SaCb1*stp*SclrNN 391 if(fw .gt. zero) tmp=tmp-two*tmp1*fw 392 if(fwp .gt. zero) tmp=tmp-tmp1*SclrNN*fwp 393 394 srcL(e) = -1.0*tmp 395 srcR(e) = srcRat(e)*SclrNN 396 enddo 397c 398c.... One source term has the form (beta_i)*(phi,_i). Since 399c the convective term has (u_i)*(phi,_i), it is useful to treat 400c beta_i as a "correction" to the velocity. In calculating the 401c stabilization terms, the new "modified" velocity (u_i-beta_i) is 402c then used in place of the pure velocity for stabilization terms, 403c and the source term sneaks into the RHS and LHS. Here, the term 404c is given by beta_i=(cb_2)*(phi,_i)/(sigma) 405c 406 407 tmp = saCb2 * saSigmaInv 408 uMod(:,1) = u1(:) - tmp * gradS(:,1) 409 uMod(:,2) = u2(:) - tmp * gradS(:,2) 410 uMod(:,3) = u3(:) - tmp * gradS(:,3) 411 412 elseif(iRANS.eq.-2) then ! K-epsilon 413c.... compute the global gradient of u 414c 415 gradV = zero 416 do n = 1, nshl 417c 418c du_i/dx_j 419c 420c i j indices match array where V is the velocity (u in our notes) 421 gradV(:,1,1) = gradV(:,1,1) + shg(:,n,1) * yl(:,n,2) 422 gradV(:,2,1) = gradV(:,2,1) + shg(:,n,1) * yl(:,n,3) 423 gradV(:,3,1) = gradV(:,3,1) + shg(:,n,1) * yl(:,n,4) 424c 425 gradV(:,1,2) = gradV(:,1,2) + shg(:,n,2) * yl(:,n,2) 426 gradV(:,2,2) = gradV(:,2,2) + shg(:,n,2) * yl(:,n,3) 427 gradV(:,3,2) = gradV(:,3,2) + shg(:,n,2) * yl(:,n,4) 428c 429 gradV(:,1,3) = gradV(:,1,3) + shg(:,n,3) * yl(:,n,2) 430 gradV(:,2,3) = gradV(:,2,3) + shg(:,n,3) * yl(:,n,3) 431 gradV(:,3,3) = gradV(:,3,3) + shg(:,n,3) * yl(:,n,4) 432c a j u a i 433c from our notes where we had N_{a,j} = dN_a/dx_j note that i is off by one because p was first in yl vector 434c 435 enddo 436 437 dwall = zero 438 do n = 1, nenl 439 dwall = dwall + shape_funct(:,n) * dwl(:,n) 440 enddo 441 442 kay(:)=zero 443 epsilon(:)=zero 444 do ii=1,npro 445 do jj = 1, nshl 446 kay(ii) = kay(ii) + shape_funct(ii,jj) * yl(ii,jj,6) 447 epsilon(ii) = epsilon(ii) 448 & + shape_funct(ii,jj) * yl(ii,jj,7) 449 enddo 450 enddo 451c no source term in the LHS yet 452 srcL(:)=zero 453 454 call elm3keps (kay, epsilon, 455 & dwall, gradV, 456 & srcRat, srcR, srcJac) 457 if(isclr.eq.1) srcL = srcJac(:,1) 458 if(isclr.eq.2) srcL = srcJac(:,4) 459 iadvdiff=0 ! scalar advection-diffusion flag 460 if(iadvdiff.eq.1)then 461 srcL(:)=zero 462 srcR(:)=zero 463 srcRat(:)=zero 464 endif 465c 466c.... No source terms with the form (beta_i)*(phi,_i) for K or E 467c 468 uMod(:,1) = u1(:) - zero 469 uMod(:,2) = u2(:) - zero 470 uMod(:,3) = u3(:) - zero 471 472 473 elseif (iLSet.ne.0) then 474 if (isclr.eq.1) then 475 476 srcR = zero 477 srcL = zero 478 479 elseif (isclr.eq.2) then !we are redistancing level-sets 480 481CAD If Sclr(:,1).gt.zero, result of sign_term function 1 482CAD If Sclr(:,1).eq.zero, result of sign_term function 0 483CAD If Sclr(:,1).lt.zero, result of sign_term function -1 484 485 sclr_ls = zero !zero out temp variable 486 487 do ii=1,npro 488 489 do jj = 1, nshl ! first find the value of levelset at point ii 490 491 sclr_ls(ii) = sclr_ls(ii) 492 & + shape_funct(ii,jj) * yl(ii,jj,6) 493 494 enddo 495 496 if (sclr_ls(ii) .gt. epsilon_ls)then 497 498 sign_levelset(ii) = one 499 500 elseif (abs(sclr_ls(ii)) .le. epsilon_ls)then 501 502 sign_levelset(ii) = sclr_ls(ii)/epsilon_ls 503 & + sin(pi*sclr_ls(ii)/epsilon_ls)/pi 504 505 elseif (sclr_ls(ii) .lt. epsilon_ls) then 506 507 sign_levelset(ii) = - one 508 509 endif 510 511 srcR(ii) = sign_levelset(ii) 512 513 enddo 514 515 srcL = zero 516 517cad The redistancing equation can be written in the following form 518cad 519cad d_{,t} + sign(phi)*( d_{,i}/|d_{,i}| )* d_{,i} = sign(phi) 520cad 521cad This is rewritten in the form 522cad 523cad d_{,t} + u * d_{,i} = sign(phi) 524cad 525 526CAD For the redistancing equation the "velocity" term is calculated as 527CAD follows 528 529 530 531 mytmp = sign_levelset / sqrt( gradS(:,1)*gradS(:,1) 532 & + gradS(:,2)*gradS(:,2) 533 & + gradS(:,3)*gradS(:,3)) 534 535 uMod(:,1) = mytmp * gradS(:,1) ! These are for the LHS 536 uMod(:,2) = mytmp * gradS(:,2) 537 uMod(:,3) = mytmp * gradS(:,3) 538 539 u1 = umod(:,1) ! These are for the RHS 540 u2 = umod(:,2) 541 u3 = umod(:,3) 542 543 544 endif ! close of scalar 2 of level set 545 546 else ! NOT turbulence and NOT level set so this is a simple 547 ! scalar. If your scalar equation has a source term 548 ! then add your own like the above but leave an unforced case 549 ! as an option like you see here 550 551 srcR = zero 552 srcL = zero 553 endif 554 555 return 556 end 557 558 559 560