CG  Version 25
lineSolveBL.h
Go to the documentation of this file.
1 ! This file ins included by insLineSolveNew.bf
2 
3 c Define the turbulent eddy viscosity and its derivatives for the BL model
4 #beginMacro defineBLDerivatives(dim,gt)
5  nuT = u(i1,i2,i3,nc)
6  #If #gt == "rectangular"
7  nuTx(0)=ux2(nc)
8  nuTx(1)=uy2(nc)
9  #If #dim == "3"
10  nuTx(2)=uz2(nc)
11  #End
12  #Else
13  #If #dim == "2"
14  nuTx(0)=ux2c(nc)
15  nuTx(1)=uy2c(nc)
16  #Else
17  nuTx(0)=ux3c(nc)
18  nuTx(1)=uy3c(nc)
19  nuTx(2)=uz3c(nc)
20  #End
21  #End
22 #endMacro
23 
24 
25 c Define the turbulent eddy viscosity and its derivatives
26 #beginMacro defineValuesBL(dim,gt)
27  ! chi3=0.
28  defineBLDerivatives(dim,gt)
29 #endMacro
30 
31 
32 c **************************************************************
33 c Macro to compute Baldwin-Lomax Turbulent viscosity
34 c **************************************************************
35 #beginMacro computeBLNuT()
36 
37  maxvt=0
38  indexRange(0,0)=n1a
39  indexRange(1,0)=n1b
40  indexRange(0,1)=n2a
41  indexRange(1,1)=n2b
42  indexRange(0,2)=n3a
43  indexRange(1,2)=n3b
44  ! assign loop variables to correspond to the boundary
45 
46  do axis=0,nd-1
47  do side=0,1
48 c write(*,*) "SIDE, AXIS, BC ",side,axis,boundaryCondition(side,axis)
49  if( boundaryCondition(side,axis).eq.noSlipWall )then
50  is1=0
51  is2=0
52  is3=0
53  if( axis.eq.0 )then
54  is1=1-2*side
55  n1a=indexRange(side,axis) !-is1 ! boundary is 1 pt outside
56  n1b=n1a
57  else if( axis.eq.1 )then
58  is2=1-2*side
59  n2a=indexRange(side,axis) !-is2
60  n2b=n2a
61  else
62  is3=1-2*side
63  n3a=indexRange(side,axis) !-is3
64  n3b=n3a
65  end if
66 
67  io(1)=0
68  io(2)=0
69  io(3)=0
70  io(axis+1)=1-2*side
71 
72  ibb=indexRange(0,axis)
73  ibe=indexRange(1,axis)-1
74 c write(*,*) ibb,ibe
75 
76  do ii3=n3a,n3b
77  do ii2=n2a,n2b
78  do ii1=n1a,n1b
79 
80  if ( ii3.ge.ktrip .and. ii2.ge.jtrip .and. ii1.ge.itrip ) then
81  i1 = ii1
82  i2 = ii2
83  i3 = ii3
84 
85  if ( nd.eq.2 ) then
86  if ( axis.eq.0 ) then
87  ditrip = ii2-jtrip
88  else
89  ditrip = ii1-itrip
90  endif
91  else
92  if ( axis.eq.0 ) then
93  ditrip = min((ii3-ktrip),(ii2-jtrip))
94  else if ( axis.eq.1 ) then
95  ditrip = min((ii1-itrip),(ii3-ktrip))
96  else
97  ditrip = min((ii1-itrip),(ii2-jtrip))
98  endif
99  endif
100 
101  ctrans = (1-exp(-ditrip/3.))**2
102 c ctrans=1
103 c write(*,*) i1,i2,i3,ctrans
104  norm(1) = 0
105  norm(2) = 0
106  norm(3) = 0
107 
108  norm(1) = rxi(axis,0)
109  norm(2) = rxi(axis,1)
110  if ( nd.eq.3 )norm(3) = rxi(axis,2)
111 
112  nmag=sqrt(norm(1)*norm(1)+norm(2)*norm(2)+norm(3)*norm(3))
113 
114  norm(1) = norm(1)/nmag
115  norm(2) = norm(2)/nmag
116  norm(3) = norm(3)/nmag
117 
118  ftan(1) = 0
119  ftan(2) = 0
120  ftan(3) = 0
121 
122  if ( nd.eq.2 ) then
123 
124  if ( gridType.eq.rectangular ) then
125 
126  ftan(1) = 2*norm(1)*ux2(uc) + norm(2)*(ux2(vc)+uy2(uc))
127  ftan(2) = norm(1)*(uy2(uc)+ux2(vc)) + 2*norm(2)*uy2(vc)
128 
129  else
130 
131  ftan(1) = 2*norm(1)*ux2c(uc) + norm(2)*(ux2c(vc)+uy2c(uc))
132  ftan(2) = norm(1)*(uy2c(uc)+ux2c(vc)) + 2*norm(2)*uy2c(vc)
133 
134  end if
135 
136  else
137 
138  if ( gridType.eq.rectangular ) then
139 
140  ftan(1)=2*norm(1)*ux2(uc)+norm(2)*(ux2(vc)+uy2(uc)) + norm(3)*(ux2(wc)+uz2(uc))
141 
142  ftan(2)=norm(1)*(ux2(vc)+uy2(uc)) + 2*norm(2)*uy2(vc) + norm(3)*(uy2(wc)+uz2(vc))
143 
144  ftan(3)=norm(1)*(ux2(wc)+uz2(uc)) + norm(2)*(uy2(wc)+uz2(vc)) + 2*norm(3)*uz2(wc)
145 
146  else
147 
148  ftan(1)=2*norm(1)*ux3c(uc)+ norm(2)*(ux3c(vc)+uy3c(uc)) + norm(3)*(ux3c(wc)+uz3c(uc))
149 
150  ftan(2)=norm(1)*(ux3c(vc)+uy3c(uc)) + 2*norm(2)*uy3c(vc) + norm(3)*(uy3c(wc)+uz3c(vc))
151 
152  ftan(3)=norm(1)*(ux3c(wc)+uz3c(uc)) + norm(2)*(uy3c(wc)+uz3c(vc)) + 2*norm(3)*uz3c(wc)
153 
154  end if
155 
156  end if
157 
158  fdotn = ftan(1)*norm(1)+ftan(2)*norm(2)+ftan(3)*norm(3)
159 
160 
161  ftan(1) = ftan(1) - norm(1)*fdotn
162  ftan(2) = ftan(2) - norm(2)*fdotn
163  ftan(3) = ftan(3) - norm(3)*fdotn
164 
165  tauw=nu*sqrt(ftan(1)*ftan(1)+ftan(2)*ftan(2)+ftan(3)*ftan(3))
166 
167 c yplus = y*yscale
168  yscale = sqrt(tauw)/nu ! assuming density=1 here...
169 
170  ymax=0
171  lmixmax=0
172  lmix2max=0
173 
174  maxumag=0
175  ulmax=0
176 
177  do i=ibb,ibe
178 
179  i1 = ii1 + io(1)*i
180  i2 = ii2 + io(2)*i
181  i3 = ii3 + io(3)*i
182  u(i1,i2,i3,nc) = 0
183 
184  if (gridType.eq.rectangular) then
185  if (nd.eq.2) then
186  vort = abs(ux2(vc)-uy2(uc))
187  else
188  vort = sqrt( (uy2(wc)-uz2(vc))*(uy2(wc)-uz2(vc)) - (ux2(wc)-uz2(uc))*(ux2(wc)-uz2(uc)) + (ux2(vc)-uy2(uc))*(ux2(vc)-uy2(uc)) )
189  end if
190  else
191  if (nd.eq.2) then
192  vort = abs(ux2c(vc)-uy2c(uc))
193  else
194  vort = sqrt( (uy3c(wc)-uz3c(vc))*(uy3c(wc)-uz3c(vc))- (ux3c(wc)-uz3c(uc))*(ux3c(wc)-uz3c(uc))+ (ux3c(vc)-uy3c(uc))*(ux3c(vc)-uy3c(uc)))
195  end if
196  end if
197 
198  yplus = dw(i1,i2,i3)*yscale
199  lmixw = vort* kbl*kbl*dw(i1,i2,i3)*dw(i1,i2,i3)*(1.-exp(-yplus/a0p))**2
200 
201 c write(*,*) "yplus, vort ",yplus, vort
202 c write(*,*) "dw, yscale, yplus, lmixw is ",dw(i1,i2,i3)," ",yscale," ",yplus," " ,lmixw
203  magu = u(i1,i2,i3,uc)*u(i2,i2,i3,uc) + u(i1,i2,i3,vc)*u(i1,i2,i3,vc)
204 
205  if ( nd.eq.3 ) magu = magu + u(i1,i2,i3,wc)*u(i1,i2,i3,wc)
206 
207  magumax = max(magu,maxumag)
208 
209  if ( (vort*kbl*dw(i1,i2,i3)*(1.-exp(-yplus/a0p))).gt.lmixmax ) then
210  ymax = dw(i1,i2,i3)
211  ulmax = magu
212  lmixmax = vort*kbl*dw(i1,i2,i3)*(1.-exp(-yplus/a0p))
213  lmix2max = lmixw
214 c write(*,*) "--",i,ymax,lmixmax,lmix2max
215  end if
216 
217  u(i1,i2,i3,nc) = lmixw
218 
219  end do ! i=ibb,ibe
220 
221 c now that we know lmixmax, ulmax and maxumag we can compute the eddy viscosity
222 
223  magumax = sqrt(magumax)
224  ulmax = sqrt(ulmax)
225 
226 c write(*,*) "ymax is ",ymax," lmix2max ",lmix2max
227  iswitch=0
228  do i=ibb,ibe
229 
230  i1 = ii1 + io(1)*i
231  i2 = ii2 + io(2)*i
232  i3 = ii3 + io(3)*i
233 
234  vto = alpha*ccp*min(ymax*lmixmax/kbl, cwk*ymax*(maxumag-ulmax)*(maxumag-ulmax)*kbl/lmixmax) / (1+5.5*(dw(i1,i2,i3)*ckleb/ymax)**6)
235 c vto = alpha*ccp*ymax*lmixmax/kbl/(1+5.5*(dw(i1,i2,i3)*ckleb/ymax)**6)
236 c write(*,*) ymax,dw(i1,i2,i3)
237 c write(*,*) (1+5.5*(dw(i1,i2,i3)*ckleb/ymax)**6)
238 c write(*,*) "i,j,k, yplus, vti, vto ",i1,i2,i3,dw(i1,i2,i3)*yscale,u(i1,i2,i3,nc), vto
239 
240 c write(*,*) yscale*dw(i1,i2,i3),u(i1,i2,i3,nc),vto,iswitch
241  if ( (iswitch.eq.0 .and. vto.lt.u(i1,i2,i3,nc)).or. iswitch.gt.0 ) then
242 c write(*,*) "switched at ",i, u(i1,i2,i3,nc), vto
243  u(i1,i2,i3,nc) = vto
244 
245  if ( iswitch.eq.0 ) iswitch = i
246  endif
247 
248  u(i1,i2,i3,nc) = ctrans*u(i1,i2,i3,nc)
249  maxvt = max(maxvt,u(i1,i2,i3,nc))
250 
251  end do ! i=ibb,ibe
252 
253  ! smooth the eddy viscosity a bit near the switch from inner to outter solutions
254  do i=max(ibb+1,iswitch-5),min(iswitch+5,ibe-2)
255 
256  i1 = ii1 + io(1)*i
257  i2 = ii2 + io(2)*i
258  i3 = ii3 + io(3)*i
259 
260 c yes, the relaxation coeff. is 1. I'm just setting it equal to the neighbors now
261 c yes, the i+1 node uses the updated version of the i node's value
262  u(i1,i2,i3,nc) = .5*(u(i1+io(1),i2+io(2),i3+io(3),nc)+u(i1-io(1),i2-io(2),i3-io(3),nc))
263 
264 c also, it seems the region for this smoothing should increase as the boundary
265 c layer increases in order to improve convergence. +- 5 was chosen through trial and
266 c error but could be made a function of iswitch or ymax for instance.
267  enddo
268 
269  else
270  do i=ibb,ibe
271  i1 = ii1 + io(1)*i
272  i2 = ii2 + io(2)*i
273  i3 = ii3 + io(3)*i
274 
275  u(i1,i2,i3,nc) = 0
276  end do
277  end if
278 
279  end do ! i3=i3a,i3b
280  end do ! i2=i2a,2b
281  end do ! i1=i1a,i1b
282 
283  ! reset values
284  if( axis.eq.0 )then
285  n1a=indexRange(0,axis)
286  n1b=indexRange(1,axis)
287  else if( axis.eq.1 )then
288  n2a=indexRange(0,axis)
289  n2b=indexRange(1,axis)
290  else
291  n3a=indexRange(0,axis)
292  n3b=indexRange(1,axis)
293  end if
294 
295  end if !bc
296 
297  end do ! do side
298  end do ! do axis
299 
300 c write(*,*) "maxvt is ",maxvt
301 #endMacro
Generated on Fri Jan 4 2013 10:19:23 for CG by   doxygen 1.8.3