...

APPENDIX A

by user

on
Category: Documents
1

views

Report

Comments

Transcript

APPENDIX A
117
APPENDIX A
Visual basic code for Crowley's cone test: 03 scheme
********************
CLS
CLEAR
DIM u(35, 35), v(35, 35), e(35, 35), d(35, 35), dx(35, 35), dy(35, 35), xp(35,
&35), yp(35, 35)
h = 1
11 = 1 35 35 ta = 1 * 288 ww = 7.2722 s = (2 * 3.141593) /
HA
VA =
(ta * ww) '****************************************************
OPEN "e:/matlab/vogx.m" FOR OUTPUT AS #2
'****************************************************
, - Velocity field
FOR i = 1 TO 35 = (i - 18) * h FOR j = 1 TO 35 v(i, j) = -ww * x NEXT j NEXT i x
FOR j = 1 TO
Y = (j - 18)
FOR i = 1 TO
uti, j) = ww
NEXT i NEXT j 35 * 11 35 * Y
'****************************************************
, - Initial cone distribution
FOR i = 1 TO 35 (i - 18) * h FOR j = 1 TO 35 Y = (j - 18) * 11 IF (x + 8) A 2 + Y A 2 <= 16 THEN eli, j) = -(25) * SQR((x + 8) A 2 + Y
ELSE eli, j) = 0 END IF NEXT j NEXT i x =
A
2) + 100 '****************************************************
, - Calculation of departure points
FOR i = 3 TO 33 x = (i - 18) * h FOR j = 3 TO 33 118
y = (j - 18) * 11
ux
(u(i + 1, j) - uti - 1, j) ) / (2 * h)
uy = (u (i, j + 1) - uti, j - 1) ) / (2 * 11)
vx = (v (i + 1, j) - v(i - 1, j) ) / (2 * h)
vy = (v (i, j + 1)
vii, j - 1) ) / (2 * 11)
uxx = (u(i + 2, j) - 2 * uti, j) + uti - 2, j) ) / (4 * h A 2) 2) uyy = {u (i, j + 2) - 2 * uti, j) + uti, j - 2) ) / (4 * 11
uxy = (u(i + 1, j + 1) - uti - 1, j + 1) - uti + 1, j - 1) & + uti - 1, j - 1) ) / (4 * h * 11) vxx = (v{i + 2, j) - 2 * vii, j) + v(i - 2, j ) ) / (4 * h "- 2) vyy = (v (i, j + 2) - 2 * vIi, j) + vIi, j - 2) ) / (4 * 11 " 2) vxy = (v (i + 1, j + 1) - v(i + 1, j - 1) - v(i - 1, j + 1) & + v(i - 1, j - 1) ) / &(4 * h * 11) A
dx ( i , j) = x - 5 * U (i , j) + « 5 " 2) / 2)
& - « 5 " 3 ) / 6 ) * « u (i, j) " 2) * uxx +
& vIi, j) * uxy + uti, j) * vx * uy + vii,
& + vIi, j) * &vy * uy) dy(i, j) = Y - 5 * vIi, j) + «5 A 2) / 2)
& « 5
3) / 6) * « u (i, j) " 2) * vxx +
& vIi, j) * vxy + uti, j) * vx * vy + vii,
& + & vIi, j) * &«vy) "2})
A
* (u (i , j) * ux + v ( i , j) * u y) u (i , j) * ( ux
2) + 2 * u (i , j) * j) * uy * ux + (v(i, j) " 2) * (uyy) A
* (u(i, j) * vx + vii, j) * vy)
u (i, j) * ux * vx + 2 * u (i, j) *
j) * uy * vx + (v(i, j) A 2) * (vyy)
FOR 1 = 1 TO 35
=
=
18)
(1 - 18)
IF dx (i, j) >
IF dy(i, j) >
NEXT 1
xt
yt
(1 -
* h
* 11
xt THEN xp(i, j) = 1 + 1
yt THEN yp(i, j)
1 + 1
NEXT j
NEXT i
'************************************
,-- Bicubic interpolation
FOR t
1 TO ta
PRINT t
FOR i
3 TO 33
FOR j = 3 TO 33
d(i,
j)
0
=0
1p = 0
50m = 0
w
50ml =
0
.-- Bilinear interpolation may be required
IF xp(i, j)
34 OR yp(i, j) = 34 THEN Ip
1
IF xp (i, j) = 3 OR yp(i, j) = 3 THEN Ip = 1
= 1p + xp(i, j)
(k - 18) * h
FOR 1 = Ip + yp(i, j)
y1 = (1 - 18) * 11
FOR k
xk
=
50mx =
1
-
2 TO xp(i,
j)
+ 1 - Ip
-
2 TO yp(i,
j)
+ 1 - 1p
119
FOR m = Ip + xp(i, j) - 2 TO xp(i, j) + 1 - Ip
xm = (m - 18) * h
IF m <> k THEN sornx = ((dx(i, j) - xm) I (xk - xm)) * sornx
NEXT m
somy = 1
FOR m = Ip + yp(i, j) - 2 TO yp(i, j) + 1 - Ip
yrn = (m - 18) * 11
IF m <> 1 THEN somy = ((dy(i, j) - yrn) I (yl NEXT m
w = sornx * somy
som = som + w * elk, 1)
soml = soml + w
NEXT 1
NEXT k
d(i, j)
= som
NEXT j
NEXT i
FOR i = 1 TO 35
FOR j = 1 TO 35
eli, j) = d(i, j)
NEXT j
NEXT i
NEXT
t
'*********************************
,-- Writing to Matlab
PRINT #2, "x= [" ; FOR i = 1 TO HA - 1 x = (i - 18) * h PRINT #2, x; ","; NEXT i i = HA x = (i - 18) * h PRINT #2, x; "];" PRINT #2, " " PRINT #2, "y=["; FOR j = 1 TO VA - 1 Y = (j - 18) * 11 PRINT #2, y; ","; NEXT j j = VA Y = (j - 18) * 11 PRINT #2, y; "];" PRINT #2, " " PRINT #2, "[X,Y]=rneshgrid(x,y);" PRINT #2, " "
PRINT #2, nz=[";
FOR j = 1 TO VA - 1
FOR i = 1 TO HA
yrn))
* somy
120
PRINT #2, c(i, j);
PRINT #2, tt tf., NEXT i PRINT #2, NEXT j
"
j
ff
= VA
FOR i = 1
PRINT #2,
PRINT #2,
NEXT i
i = HA
PRINT #2,
PRINT #2,
PRINT #2,
PRINT #2,
CLOSE 2
TO HA - 1
c(i, j)i
It
II;
c(i, j)i
"J ; "
"
tt
"surf(X,Y,Z)"
'****************************************************
END
'****************************************************
'****************************************************
121
APPENDIX B
Fortran code for Smolarkiewicz's deformational flow: D3 scheme
******************** program smolar integer A, L, CI, IER, t, Ip, xp(lOl,lOl), yp(lOl,lOl),ta
real tt,aa
real pi, kk, tsi(lOl,lOl), c(lOl,lOl), ALEV(13),KEELY(lO), TR(6)
real u(lOl,lOl), v(lOl,lOl), s, dx(101,lOl), dy(101,101)
real zt, d(101,lOl), som, som1, xk, yl, xm, ym, sornx, somy, w
real f(101,101),b(101,101)
character*(4) LABEL
A=8 L=lOO pi=3.141592654 kk=4*pi/L s=0.7 c-"'*
Iterations correspond to the times in the paper of Staniforth et al. (1987)
c
c
c
c
c
ta=19
ta=38
ta=57
ta=75
ta=377
ta=3768
c-"'*
Calculate grid paint values of.. ........... . do i
x =
do j
Y
= 1,
(i-1)
= 1,
(j-1)
101 101 c-"'* Streamlines
tsi(i,j) = A*(sin(kk*x))*(cos(kk*y))
c-"'* Velocity field
u(i,j) = A*(kk)*(sin(kk*x))*(sin(kk*y)) v(i,j) = A*(kk)*(cos(kk*x))*(cos(kk*y)) enddo enddo c-"'* Calculate the grid point values of initial cone distribution
do i = 1, 101 x = (i-1) do j = 1, 101 Y = (j-1) c(i,j) = 0.0 if ((x-50)**2 + (y-50)**2 .le. 225.0) c(i,j) = & ( 1.0/15.0)*(((x-50)**2 + (y-50)**2)**(O.5)) +1 b(i,j)=c(i,j) enddo enddo c*********************
c-"'*
Determination of departure paints
do i=2, 100 .
122
x=i-l do j=2,100 y=j-l (u(i+l,j)-u(i-l,j)}/2.
ux
(u(i,j+l)-u(i,j-l»/2.
uy
vx
(v(i+l,j)-v(i-l,j»/2.
vy = (v(i,j+l)-v(i,j-l»/2.
uxx = (u(i+2,j)-2.*u(i,j)+u(i-2,j»/4.
uyy
(u(i,j+2)-2.*u(i,j)+u(i,j-2»/4.
uxy = (u(i+l,j+l)-u(i-l,j+l)-u(i+l,j-l)+u(i-l,j-l»/4.
vxx = (v(i+2,j}-2.*v(i,j)+v(i-2,j»/4.
vyy = (v(i,j+2)-2.*v(i,j)+v(i,j-2»/4.
vxy = (v(i+l,j+l)-v(i+l,j-l)-v(i-l,j+l)+v(i-l,j-l»/4.
dX(i,j)=x-s*u(i,j)+«s**2)/2.)*(u(i,j)*ux+
&v (i, j) *uy) - ( (5 * * 3) /6. ) * ( (u (i, j ) * * 2) *uxx+u (i, j )
&*(ux**2.}+2.*u(i,j)*v(i,j)*uxy+u(i,j)*vx*uy+
&v(i,jl*uy*ux+(v(i,j)**2)*(uyyl+v(i,j)*vy*uy)
dy(i,j)=y-s*v(i,j)+«s**2)/2.)*(u(i,j)*vx+
&v(i,j)*vy)-«s**3)/6.)*«u(i,j)**2)*vxx+u(i,j)
&*ux*vx+2.*u(i,j)*v(i,j)*vxy+u(i,j)*vx*vy+
&v(i,j)*uy*vx+(v(i,j)**2)*(vyy)+v(i,j)*«vy)**2»
do 1=1,101 xt= (1-1) yt= (1-1) if (dx(i,j) .gt.xt) xp(i,j)=1+1 if (dy(i,j).gt.yt) yp(i,j)=1+1 enddo enddo enddo c*******************************************************
c*******************************************************
c-* Bicubic interpolation
do t=l,ta
do i=2,100 do j=2,100 d(i,j)=O
w=0
Ip=O som=O soml=O if
if
if
if
(xp(i,j)
(yp(i,j)
(xp(i,j)
(yp(i,j)
.eq.l0l) Ip=1 .eq.l01) Ip=l .eq.2) Ip=1 .eq.2) Ip=1 do k=lp+xp(i,j)-2,xp(i,j)+1-lp xk=k-l do 1=lp+yp(i,j)-2,yp(i,j)+1-1p yl=1-1 somx=1 do m=lp+xp(i,j)-2,xp(i,j)+1-1p 123
xm=m-l if(m.ne.k) sornx=«dx(i,j)-xm)/(xk-xm»*sornx enddo somy=l do m=lp+yp(i,j)-2,yp(i,j)+1-lp ym=m-l if (m.ne.l) somy=«dy(i,j)-ym)/(yl-ym»*somy enddo w=sornx* 5 omy som=som+w*c(k,l) soml=soml+w enddo enddo d(i,j)=som
enddo enddo do i=l,lOl do j=l,lOl c(i,j)=d(i,j) enddo enddo enddo
c********************************************************
c**** Conservation properties
conl=O con2=0 con3=0 cona=O conb=O conc=O do i=l,lOl do j=l,lOl conl=conl+c(i,j) con2=con2+(c(i,j»**2 con3=con3+abs(c(i,j» cona=cona+b(i,j) conb=conb+(b(i,j»**2 conc=conc+abs(b(i,jll enddo enddo conl=conl/cona con2=con2/conb con3=con3/conc print*, conl,con2,con3 c*********************************************************
c*********************************************************
c**** Writing out for MATLAB
20
open(1,file='/frae/WEERLIG/MAT/matd36.m',form='formatted'
&,recl=lOOOOO)
write(1,20) c
format (101f20.10)
124
c-'-* Writing statistics
open(2,file='/frae/WEERLIG/STA/stad36') write(2,30) con1 write(2,30) con2 write(2,30) con3 30
format (f20.10)
c*********************************************************
c**** Interactive display using PGPLOT
c**** Open plot
IF (PGOPEN('?')
.LT. 1) STOP
c**** Set contours of streamlines to be plotted
do i = 1,13 ALEV(I) = -6 + (i-1) enddo c**** Set contours of cone to be plotted
KEELY(l)
KEELY (2)
KEELY(3)
KEELY(4)
KEELY(S)
=
=
=
0.02 0.1 0.2 0.3 0.4 KEELY ( 6 ) = O. 5 KEELY(7)
KEELY(S)
KEELY(9)
KEELY (10)
= 0.6 0.7 O.S = O. 9 c**** Set axis
CALL PGENV(0.,100.,0.,100.,l,1) TR( 1) = -1. 0 TR(2) = 1. 0 TR(3) = 0.0 TR (4)
-1. 0 TR(5) = 0.0 TR(6) = 1.0 c**** Label axis, title
CALL PGLAB('x',
'y', 'Smolarkiewicz's Test
D3 scheme
c-'-* Set colour and plot contours of streamlines
CALL
CALL
CALL
CALL
PGSCI(3) PGCONT(tsi,101,101,1,101,1,101,ALEV,-7,TR) PGSCI(4) PGCONT(tsi,101,101,1,101,1,101,ALEV(8),-6,TR) c**** Set colour and plot contours of cone
CALL PGSCI(l) CALL PGCONT(b,101,101,1,101,1,101,KEELY,10,TR) c-.-* Set colour and plot contours of solution
CALL PGSCI(2) CALL PGCONT(c,lOl,lOl,l,lOl,l,lOl,KEELY,lO,TR) c**** CLOSE PLOT
CALL PGCLOS
end
t=T')
Fly UP