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