1CC*************************************************************************** 2CC Copyright (c) 1994-1999 ACUSIM Software, Inc. 3CC All rights reserved. 4CC This source code is confidential and may not be disclosed. 5CC*************************************************************************** 6 7CC=========================================================================== 8CC 9CC "elm3Keps.f": 3D k-epsiron (LB model) turbulence element formation. 10CC 11CC Original: Farzin Shakib (Nov 98) 12CC K-epsilon modification by Je-Hoon, Kim (Jan 98) 13CC Wlicox's book, pp.140, Lam-Bremhorst model 14CC 15CC Includes k-epsilons stagnation abnormally correction (Feb 99) 16CC ( Durbin, P. A., "On the k-3 stagnation point 17CC abnormally", pp.89-90, 1996, Int. J. of Heat and Fluid flow). 18CC 19 20 21c... Note 22 23c.. 1. K-epsilon model seems more sensitive to grid resolution than one equation models. 24c It seems good (fine enough) grid seldom get blow-up ('nan' error) while coarse grids 25c sometimes blow up. It is strongly recommended to have very fine grid near the wall 26c at least y+ (1st grid off the wall) less than or = 0.5. 27 28c.. 2. when it blew up while converging to the solution, alternative srcJac's were 29c.. used to prevent blow up. It eliminated blow up problem, but it does not seem 30c.. to be a good Jacobian to solve the problem from the very beginning (just good 31c.. as a restarter). 32c.. Further numerical test will be done to find a better Jacobian 33c.. which is good for either case, I have not finished it yet. 34c.. I have coded all possible terms to play with and am studying the behavior now. 35 36c.. 3. Durbin stated that we could limit the T2nd (second time scale) to control 37c.. the abnormal growth of k , but it does not seem to help much. 38 39c.. 4. About good IC's for K-epsilon model, 40c.. K_(O) = k+ max(=4.5~6) * (u_tau)^2 41c.. eps_(O) = eps+ max(=0.1~0.25) * (u_tau)^4 / kinematic_viscosity 42c.. velocity(0)= u_tau *( 1/0.418*Log(y+ max)+5.5 ) or mass velocity if it can be estimated. 43c.. mu_t(0) = (500~1000) * mu 44 45 46 47c.. (example) for periodic channel problem 48 49c.. when u_tau=0.707, k+ max = 6, k IC = 3 50c.. eps+ max = 0.1, eps IC = 1250 51c.. vel IC ~ 20 52 53c.. when u_tau=1, k+ max = 6, k IC = 6 54c.. eps+ max = 0.1, eps IC = 5000 55c.. vel IC ~ 30 56c.. u_tau=1's IC's works for lower dpdx cases. 57 58 59c.. I used possible maximum for k and epsilon. 60c.. From one of review papers on k-epsilon models, we know approximate k+ max 61c.. and epslon+ max for attached flows. 62c.. Shear velocity (utau) could also be estimated for particular probelm which 63c.. actually set k and epsilon along with another flow properties. 64c.. Too small epsilon or too small mu_t might let k grow unchecked therefore 65c.. make the flow field unstable. 66 67c.. 5. Restart works better from higher Reynolds number to lower Reynolds number flows but 68c.. not the other way (Before Durbin's modification, and I haven't check after Durbin's 69c.. modification). 70 71c.. 6. In terms of terminology and naming convention, I pretty much followed DR. Shakib's 72c.. way. If there are any confusion, pls let me know. 73 74 75c.. 7. Slight modifications are done to Jacobians (basically sign control) to accelerate 76c.. convergence and to prevent a blow-up before the solutions converge 77 78 79c.. 8. Periodic channel problem tells : 80c.. The lower bound of y+ of the 1st node is 0.6 and 12 nodes within y+ = 20 (near wall region) 81c.. (STAR-CD recommend y+~1 and 20 nodes within y+=20). 82c.. Our 1st y+ requirement is sligtly more than what STAR-CD recommends. 83 84 85c.. 9. Fin tube problem (pls neglect the sample problem. it will not converge) tells : 86c.. By inspecting y+ resolution, when the solver blowed up, there were 87c.. about 6~7 points which had y+ (1st node) more than 0.6 near the stagnation points. 88c.. The GR's also did not go down below GR of 1.e-4. 89 90c.. After fixing y+ resolution near the stagnation points, it converged well (below 1.e-4). 91c.. 92c.. (example) IC selection for fin tube 93 94c.. To find out candidiate Initial Condition for k and epsilon, I assumed that 95c.. the wall shear stress could reach upto 10 times the average wall shear stress 96c.. at the stagnation point (I think it will vary of course, depending on the problem). 97c.. 98c.. u_tau_average can be estimated from experimaental data. 99c.. u_tau_max = SQRT(10) * u_tau_average 100c.. k+ max ~ 6, epsilon max ~0.25 101c.. 102CC=========================================================================== 103 104 105C============================================================================ 106C 107C "fElm3KepsCoef": 108C 109C============================================================================ 110 subroutine elm3keps (kay, epsilon, 111 & dwi, gradV, 112 & srcRat1,src1, srcJac) 113c 114c.... Data declaration 115c 116 include "common.h" 117 real*8 kay(npro), epsilon(npro), 118 & dwi(npro) 119 real*8 gradV(npro,3,3) 120 121 real*8 rmu(npro), 122 & rho(npro) 123 real*8 srcRat(npro,2), 124 & src(npro,2), srcJac(npro,4), 125 & srcRat1(npro), src1(npro) 126c 127 integer advdiff 128 integer e 129 real*8 fct1 130 real*8 tmp1, tmp2 , tmp0 131c 132 real*8 k, kInv, kq, eps, y, epsInv, epsP2, 133 & ss, mut, mut_k, mut_eps, rat, jac 134 135 136 real*8 Ce1, Ce2, epsP2Inv, CmuKE, sigmaKE, nu, 137 & kk, kqInv, nuInv, Rey, Ret, Rey_k, 138 & Ret_k, Ret_eps, ff1, f2, kkInv, 139 & RetInv, fmukeInv, fmukeP2Inv 140 141 real*8 fmuke, fmuke_k, fmuke_eps 142 143c... Addings due to Durbin's correction 144 145 real*8 T1st, T2nd, T2ndInv, two3rdCmuInv, tri8thq, 146 & ssq, ssqInv, two3rdq, pk1, pk2, 147 & pk, pk_k, pk_eps, Tscale,Tinv, 148 & Tscale_k, Tscale_eps, ff1_k, ff1_eps, f2_k,f2_eps, 149 & ff1Inv, f2Inv 150 151 152 153c 154c include "elmf.h" 155c 156c.... Compute src and its jacobians 157c 158c.... set Lam-Bremhorst' k-epsilon Model constants 159 rho(:)=datmat(1,1,1) 160 rmu(:)=datmat(1,2,1) 161 Ce1 = 1.44 162 Ce2 = 1.92 163 CmuKE = 0.09 164 sigmaKE = 1.3 165 166 two3rdCmuInv = 2./3./CmuKE 167 tri8thq = SQRT(3./8.) 168 two3rdq = SQRT(2./3.) 169c 170 advdiff = 0 171 if(advdiff.eq.0)then ! not advection-diffusion 172 do e = 1, npro 173 174 nuInv = rho(e)/rmu(e) 175 k = abs(kay(e)) 176 if (k.lt.1.e-32) k=0 177 eps = abs(epsilon(e)) 178 y = dwi(e) 179 180c--------------patch 181c if(k.gt.0.45) k=0.45 182c if(eps.gt.2158) eps=2158 183c-------------------------------- 184 185c 186 kInv = 0 187 kq = sqrt(k) 188 kqInv = 0 189 kkInv = 0 190 191c.... limiting k. Instead of saying k.ne.0 192 193 if ( k .gt.1.e-32 ) then 194 kInv = 1. / k 195 kqInv = 1./sqrt(k) 196 kkInv = kInv*kInv 197 endif 198 199 kk = k * k 200 201 epsP2 = eps * eps 202 epsInv = 0 203 epsP2Inv = 0 204 205 206c..... limiting epsilon. Instead of saying eps.ne.0 207 208 if ( eps .gt.1.e-32) then 209 epsInv = 1. / eps 210 epsP2Inv = epsInv*epsInv 211 endif 212 213c 214 ss = gradV(e,1,1) ** 2 215 & + gradV(e,2,2) ** 2 216 & + gradV(e,3,3) ** 2 217 & + 0.5 218 & * ( (gradV(e,2,3) + gradV(e,3,2)) ** 2 219 & + (gradV(e,3,1) + gradV(e,1,3)) ** 2 220 & + (gradV(e,1,2) + gradV(e,2,1)) ** 2 ) 221c 222 ssq = sqrt(ss) 223 ssqInv = 0 224 225 if(ss.ne.0) ssqInv = 1./sqrt(ss) 226 227 Rey = kq * y * nuInv 228 Ret = kk * epsInv * nuInv 229 RetInv = 0 230 231c... limitng Ret so that it does not get 'nan' error 232 233 if(Ret.lt.1.d100.AND.Ret.gt.zero) RetInv=1./Ret 234 235 Rey_k = 0.5 * y * kqInv * nuInv 236 Ret_k = 2. * k * epsInv * nuInv 237 Ret_eps = -kk * epsP2Inv * nuInv 238 239 tmp1 = exp(-0.0165*Rey) 240 241 fmuKE = (1. -tmp1) ** 2 * (1.+20.5*RetInv) ! fmu of Lam-Bremhorst 242 243 fmuKEInv = 0.0 244 fmuKEP2Inv = 0.0 245c.... limiting fmuKE. fmuke max ~ 1, and it could get very small near the wall. 246 247 if(fmuKe.gt.1e-32) then 248 fmuKEInv = 1./fmuke 249 fmuKEP2Inv = fmuKEInv*fmuKEInv 250 endif 251 252 fmuKE_k = 0.033*(1.-tmp1)*(1.+20.5*RetInv)*Rey_k*tmp1 253 & -20.5 *(1.-tmp1)**2 * Ret_k*RetInv**2 254 255 fmuKE_eps= -20.5*(1.-tmp1)**2* Ret_eps*RetInv**2 256 257 ff1 = 1. + ( 0.05*fmuKEInv) ** 3 ! f1 as in Lam-Bremhorst 258 f2 = 1. - exp(- Ret ** 2 ) ! f2 as in Lam-Bremhorst 259 260 ff1Inv=zero 261 f2Inv=zero 262 if(ff1.gt.1.0e-32) ff1Inv=1./ff1 263 if(f2.gt.1.0e-32) f2Inv =1./f2 264 265 266 ff1_k = -0.000375*fmuKE_k *fmuKEInv**4 267 ff1_eps = -0.000375*fmuKE_eps*fmuKEInv**4 268 269 f2_k = 2.* Ret * Ret_k * exp(-Ret**2) 270 f2_eps = 2.* Ret * Ret_eps * exp(-Ret**2) 271 272 T1st = k * epsInv ! 1st time scale 273 T2nd = fmuKEInv*two3rdCmuInv*tri8thq*ssqInv ! 2nd time scale 274 275c... Depending on the choice of T (to limit k growth near stagnation) 276 277 if (T1st.lt.T2nd) then 278 279 Tscale = T1st 280 Tinv = 0 281 282c if(Tscale.ne.0) Tinv=1./Tscale 283 284c... Limiting time scale s.t. Tinv**4 not go over 1.e160 285 286 if(Tscale.gt.1.e-32) Tinv=1./Tscale 287cccccTHIS IS WHAT WAS HERE WHEN I GOT IT FROM JE HOON 288ccccc Tscale_k = fmuKE_k*k*epsInv+fmuKE*epsInv 289ccccc Tscale_eps = fmuKE_eps*k*epsInv -fmuKE*k*epsP2Inv 290 Tscale_k= epsInv ! Acusolve's choice 291 Tscale_eps= -k*epsP2Inv ! AcuSolve's choice 292 293 else 294 295 Tscale = T2nd 296 Tinv = 0 297 298 if(Tscale.gt.1.e-32) Tinv=1./Tscale 299cccccTHIS IS WHAT WAS HERE WHEN I GOT IT FROM JE HOON 300ccccc Tscale_k = 0 301ccccc Tscale_eps =0 302 Tscale_k= two3rdCmuInv*tri8thq*ssqInv* 303 & (-fmuKEP2Inv*fmuKE_k) ! acusolve's choice 304 Tscale_eps= two3rdCmuInv*tri8thq*ssqInv* 305 & (-fmuKEP2Inv*fmuKE_eps) ! acusolve's choice 306 307 endif 308 309c... After Limiting all values, i feel it's unnecessary to limit jacobians 310c... since they are already limited. 311c... Two other routines which defines these quantities are the same 312 313 314 315 316c... recall that tmp1= exp(-0.0165*Rey) 317 318 319 320 mut = rho(e)*CmuKE*fmuKE* k *Tscale 321 322 mut_k = rho(e)*CmuKE*( 323 & fmuKE_k*k*Tscale 324 & + fmuKE *Tscale 325 & + fmuKE *k*Tscale_k 326 & ) 327 328 mut_eps = rho(e) * CmuKE *k*( 329 & fmuKE_eps*Tscale+ 330 & fmuKE *Tscale_eps 331 & ) 332 333 tmp0 = 2 * ss 334 335 336 pk1 = mut*tmp0 337 pk2 = rho(e)*two3rdq*k*sqrt(ss) 338 339 340 if (pk1.lt.pk2) then 341 pk = pk1 342 pk_k = mut_k * tmp0 343 pk_eps = mut_eps * tmp0 344 else 345 pk = pk2 346 pk_k = kInv*pk2 347 pk_eps =0 348 endif 349 350 351 src(e,1) = pk - rho(e) * eps 352 src(e,2) = ff1 * Ce1 * pk * Tinv 353 & -rho(e)* Ce2 * f2 * eps * Tinv 354 355 srcRat(e,1) = -kInv*( pk - rho(e) * eps ) 356 srcRat(e,2) = -epsInv*( 357 & ff1 * Ce1 * pk * Tinv 358 & -rho(e)* Ce2 * f2 * eps * Tinv 359 & ) 360 361 srcJac(e,1) = -(pk_k) ! jacobian for k PDE alone 362 srcJac(e,3) = -(pk_eps-rho(e)) ! d(Fsrck)/d(epsilon) 363 srcJac(e,2) = -( 364 & ff1_k * Ce1 * pk * Tinv 365 & +ff1 * Ce1 * pk_k * Tinv 366 & -ff1 * Ce1 * pk * Tscale_k * Tinv**2 367 & -rho(e) * Ce2 * f2_k * eps * Tinv 368 & +rho(e) * Ce2 * f2 * eps * Tscale_k * Tinv**2 ! do not touch 369 & ) ! d(Fsrce)/d(k) 370 srcJac(e,4) = -( 371 & ff1_eps * Ce1 * pk * Tinv 372 & +ff1 * Ce1 * pk_eps * Tinv 373 & -ff1 * Ce1 * pk * Tscale_eps * Tinv**2 374 & -rho(e) * Ce2 * f2_eps * eps * Tinv 375 & -rho(e) * Ce2 * f2 * Tinv 376 & +rho(e) * Ce2 * f2* eps * Tscale_eps * Tinv**2 377 & ) ! jacobian for epsilon PDE alone 378 379 380 381 382 enddo 383c 384c.... Ensure positivity of srcJac 385c 386 do e = 1, npro 387 388 389 srcJac(e,1) = max( srcJac(e,1), srcRat(e,1), 0.d0 ) 390 srcJac(e,4) = max( srcJac(e,4), srcRat(e,2), 0.d0 ) 391 392c write(777,*) e, srcJac(e,1), srcJac(e,4) 393 394c 395 tmp1 = min( srcJac(e,1) * srcJac(e,4), 396 & (srcJac(e,1)-srcRat(e,1)) * 397 & (srcJac(e,4)-srcRat(e,2)) ) 398 tmp2 = srcJac(e,2) * srcJac(e,3) 399 400 401 if ( tmp2 .gt. tmp1 ) then 402 tmp2 = sqrt(tmp1/tmp2) 403 srcJac(e,2) = tmp2 * srcJac(e,2) 404 srcJac(e,3) = tmp2 * srcJac(e,3) 405 endif 406C srcJac(e,2) = 0 407C srcJac(e,3) = 0 408 enddo 409 if(isclr.eq.1)then ! kay 410 srcrat1 = srcrat(:,1) 411 src1 = src(:,1) 412 else if (isclr.eq.2) then ! epsilon 413 srcrat1 = srcrat(:,2) 414 src1 = src(:,2) 415 endif 416 417 else ! Advection-diffusion 418c Advection-diffusion case 419 srcrat1 = zero 420 src1 = zero 421 srcjac = zero 422 endif 423c 424c.... Compute viscosity 425c 426c do e = 1, npro 427c viscTot(e,1) = rmu(e) + xmuT(e) 428c 429c viscTot(e,2) = rmu(e) + xmuT(e) 430c & / (SigmaKE) 431c enddo 432c 433c.... Compute PDE residual 434c 435c do e = 1, nElems 436c pdeRes(e,1) = dens 437c 1 * ( masFct * tkeTd(e) 438c 2 + velK(e,1) * gradK(e,1) 439c 3 + velK(e,2) * gradK(e,2) 440c 4 + velK(e,3) * gradK(e,3) ) 441c 5 - src(e,1) 442c 6 - diffK(e) 443c pdeRes(e,2) = dens 444c 1 * ( masFct * tepsTd(e) 445c 2 + veleps(e,1) * gradeps(e,1) 446c 3 + veleps(e,2) * gradeps(e,2) 447c 4 + veleps(e,3) * gradeps(e,3) ) 448c 5 - src(e,2) 449c 6 - diffeps(e) 450c enddo 451c 452c.... end of fElm3KepsCoef() 453c 454 return 455 end 456