1C 2C 3C Portions Copyright (C) 2001 Joseph H. Morrison 4C 5C This file is part of ISAAC. 6C 7C This program is distributed under the terms of the ISAAC Public Source 8C License. This program is distributed WITHOUT ANY WARRANTY; without 9C even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 10C PURPOSE. 11C 12C You should have received a copy of the ISAAC Public Source License 13C with this program. If you did not, you may get a copy of the license 14C at http://isaac-cfd.sourceforge.net. 15C 16C $Revision: 4.3 $ 17C $Author: jhmorr $ 18C $Date: 2001/10/29 03:25:36 $ 19C $State: Exp $ 20C $Log: komega.F,v $ 21C Revision 4.3 2001/10/29 03:25:36 jhmorr 22C Updated license information 23C 24C Revision 4.2 2001/06/08 04:56:19 jhmorr 25C Added notice 26C 27C Revision 4.1 1998/04/14 20:35:59 jhmorr 28C Alpha 4.1. 29C 30C 31 32 33 34 35 36 37 SUBROUTINE RMUTKW (IDIM, JDIM, KDIM, I, Q, RMUT) 38C 39C Routine to calculate the turbulent eddy viscosity for the k-omega 40C turbulence model. 41C 42C 43C IDIM,JDIM,KDIM : Dimensions of current block 44C I : Plane to do calculation at 45C Q : Primitive variables at cell centers 46C RMUT : Turbulent eddy viscosity at cell centers 47C of the I-plane 48C 49 include '../../header/common.h' 50C 51 DIMENSION Q (0:JDIM+2,0:KDIM+2,0:IDIM+2,NQ), 52 1 RMUT (0:JDIM+2,0:KDIM+2,0:IDIM+2) 53C 54C Calculate the turbulent eddy viscosity at cell centers from values 55C of turbulent kinetic energy (k) and specific dissipation rate (omega) 56C 57 DO 120 K = 0, KDIM + 2 58 DO 110 J = 0, JDIM + 2 59 RHO = Q(J,K,I,1) 60 TKE = Q(J,K,I,6) 61 OMEGA = Q(J,K,I,7) 62 RMUT(J,K,I) = CMU * RHO * TKE / (OMEGA + RSMALL) * 63 1 RE / FSMACH 64 110 CONTINUE 65 120 CONTINUE 66C 67C 68C Finished calculating RMUT from k and omega 69C 70 RETURN 71 END 72 73 74 75 76 77 78 SUBROUTINE SRCKW (NPTS, Q, PROPS, DQDX, DQDY, DQDZ, TAU, SRC, 79 1 NPRLIM, NPRNEG) 80C 81C Routine to calculate the source terms for the k-omega turbulence model. 82C 83C This routine is organized as follows: 84C 1. Form the source terms for the k-epsilon model 85C 86C NPTS : Number of points to calculate source terms at 87C Q : Primitive variables at cell centers 88C PROPS : Properties stored at cell centers 89C PROPS(1) = RMU molecular viscosity 90C PROPS(2) = RMUT turbulent eddy viscosity 91C PROPS(3) = YPLUS Y+ 92C DQDX,DQDY,DQDZ : Derivatives of Q at cell centers 93C TAU : Reynolds stresses 94C SRC : Source terms for the k-omega model at the cell centers 95C of the I-plane 96C NPRLIM : Counter for number of times production limited 97C NPRNEG : Number of times invoke positivity preservation for production 98C 99 include '../../header/common.h' 100C 101 DIMENSION Q (NPTS,NQ), 102 1 PROPS (NPTS,NP) 103C 104 DIMENSION DQDX (NPTS,NQ), 105 1 DQDY (NPTS,NQ), 106 2 DQDZ (NPTS,NQ), 107 3 TAU (NPTS,6), 108 4 SRC (NPTS,NF) 109C 110C 1. Form the source terms for the k-omega model 111C 112 DO 100 I = 1, NPTS 113C 114 RHO = Q(I,1) 115 TKE = Q(I,6) 116 OMEGA = Q(I,7) 117C 118C Production term: 119C PROD = TAUXX*DUDX + TAUYY*DVDY + TAUZZ*DWDZ + 120C TAUXY*(DUDY+DVDX) + TAUXZ*(DUDZ+DWDX) + TAUYZ*(DVDZ+DWDY) 121C 122 PROD = TAU(I,1) * DQDX(I,2) 123 1 + TAU(I,2) * DQDY(I,3) 124 2 + TAU(I,3) * DQDZ(I,4) 125 3 + TAU(I,4) * (DQDY(I,2) + DQDX(I,3)) 126 4 + TAU(I,5) * (DQDZ(I,2) + DQDX(I,4)) 127 5 + TAU(I,6) * (DQDZ(I,3) + DQDY(I,4)) 128C 129C Alternate Form of Production Term: PROD = mu_t * Vorticity^2 130C 131C WXY = 0.5E0 * (DQDY(I,2) - DQDX(I,3)) 132C WXZ = 0.5E0 * (DQDZ(I,2) - DQDX(I,4)) 133C WYZ = 0.5E0 * (DQDZ(I,3) - DQDY(I,4)) 134C WVORT = 4.E0 * (WXY*WXY + WXZ*WXZ + WYZ*WYZ) 135C RMUT = PROPS(I,2) 136C PROD = RMUT * WVORT * FSMACH / RE 137C 138C Transition specification 139C 140 IF (PROPS(I,4) .LT. 0.E0) PROD = 0.E0 141C 142C Positivity preservation for production 143C 144 IF (POSPRD) THEN 145 IF (PROD .LT. 0.E0) THEN 146 NPRNEG = NPRNEG + 1 147 PROD = ABS (PROD) 148 ENDIF 149 ENDIF 150C 151 PRODUL = PROD 152C 153C Limit production term for robustness: 154C 155 TSTDIS = PRDLIM * BSTRKW * RHO * TKE * OMEGA 156 IF (PROD .GT. TSTDIS) THEN 157 NPRLIM = NPRLIM + 1 158 PROD = TSTDIS 159 ENDIF 160C 161C Production for omega equation: use either limited or unlimited form 162C 163 PRODW = PRDE * PROD + PRDEM1 * PRODUL 164C 165C Calculate Source Term 166C 167 SRC(I,1) = 0.E0 168 SRC(I,2) = 0.E0 169 SRC(I,3) = 0.E0 170 SRC(I,4) = 0.E0 171 SRC(I,5) = 0.E0 172 SRC(I,6) = PROD - BSTRKW * RHO * TKE * OMEGA 173 SRC(I,7) = GKW * OMEGA / TKE * PRODW - BKW * RHO * OMEGA*OMEGA 174 100 CONTINUE 175C 176C Finished with k-omega source terms 177C 178 RETURN 179 END 180