1 SUBROUTINE VELOXY(nnod,nn,nr,nl,x,y,rex,rey,uk,vk,dfi,
2 * dt,tmin,tmax)
3 IMPLICIT double precision(A-H,O-Z)
4 include 'type'
5 real*8 x(NDIM),y(NDIM),rex(NDIM),rey(NDIM)
6 real*8 uk(NDIM),vk(NDIM),dfi(NDIM)
7 real*8 drx(NDIM),dry(NDIM),sr(NDIM),grad(NDIM)
8 real*8 dt,tmin,tmax
9 DO I=NR+2,NL-2
10 DU=DDERIV(REX(I-2),'*',0)
11 DV=DDERIV(REY(I-2),'*',0)
12 DX=DDERIV(X(I-2),'*',0)
13 DY=DDERIV(Y(I-2),'*',0)
14 Ds=DSQRT(Dx*Dx+Dy*Dy)
15 DFI(I)=DV/DS
16 UK(I)=(Du/Ds)*(Dx/Ds)+(Dv/Ds)*(Dy/Ds)
17 VK(I)=(Du/Ds)*(DY/Ds)-(Dv/Ds)*(DX/Ds)
18 drx(i)=dx
19 dry(i)=dy
20 sr(i)=Dsqrt((x(i+1)-x(i))**2+(Y(I+1)-Y(I))**2)
21 grad(i)=uk(i)*uk(i)+vk(i)*vk(I)
22 end do
23 Du=DDERIV(rex(NR-4),'-',1)
24 DV=DDERIV(REY(NR-4),'-',1)
25 Dx=DDERIV(X(NR-4),'-',1)
26 Dy=DDERIV(Y(NR-4),'-',1)
27 Ds=Dsqrt(Dx**2+Dy**2)
28 DFI(nr)=DV/DS
29 UK(NR)=(Du/Ds)*(Dx/Ds)+(Dv/Ds)*(Dy/Ds)
30 VK(NR)=(Du/Ds)*(DY/Ds)-(Dv/Ds)*(DX/Ds)
31 drx(nr)=dx
32 dry(nr)=dy
33 sr(nr)=Dsqrt((x(nr+1)-x(nr))**2+(Y(nr+1)-Y(nr))**2)
34 grad(nr)=uk(nr)*uk(nr)+vk(nr)*vk(nr)
35 Du=DDERIV(rex(NR),'+',2)
36 DV=DDERIV(REY(NR),'+',2)
37 Dx=DDERIV(X(NR),'+',2)
38 Dy=DDERIV(Y(NR),'+',2)
39 Ds=Dsqrt(Dx**2+Dy**2)
40 DFI(nr+1)=DV/DS
41 UK(NR+1)=(Du/Ds)*(Dx/Ds)+(Dv/Ds)*(Dy/Ds)
42 VK(NR+1)=(Du/Ds)*(DY/Ds)-(Dv/Ds)*(DX/Ds)
43 drx(nr+1)=dx
44 dry(nr+1)=dy
45 sr(nr+1)=Dsqrt((x(nr+2)-x(nr+1))**2+
46 * (Y(nr+2)-Y(nr+1))**2)
47 grad(nr+1)=uk(nr+1)*uk(nr+1)+vk(nr+1)*vk(nr+1)
48 Du=DDERIV(REX(NL),'+',1)
49 DV=DDERIV(REY(NL),'+',1)
50 Dx=DDERIV(X(NL),'+',1)
51 Dy=DDERIV(Y(NL),'+',1)
52 Ds=Dsqrt(Dx**2+Dy**2)
53 DFI(nl)=DV/DS
54 uk(nl)=(Du/Ds)*(Dx/Ds)+(Dv/Ds)*(Dy/Ds)
55 VK(NL)=(Du/Ds)*(DY/Ds)-(Dv/Ds)*(DX/Ds)
56 Du=DDERIV(REX(NL-4),'-',2)
57 DV=DDERIV(REY(NL-4),'-',2)
58 Dx=DDERIV(X(NL-4),'-',2)
59 Dy=DDERIV(Y(NL-4),'-',2)
60 Ds=Dsqrt(Dx**2+Dy**2)
61 DFI(nl-1)=DV/DS
62 uk(nl-1)=(Du/Ds)*(Dx/Ds)+(Dv/Ds)*(Dy/Ds)
63 VK(NL-1)=(Du/Ds)*(DY/Ds)-(Dv/Ds)*(DX/Ds)
64 drx(nl-1)=dx
65 dry(nl-1)=dy
66 sr(nl-1)=Dsqrt((x(nl-2)-x(nl-1))**2+
67 * (Y(nl-1)-Y(nl-2))**2)
68 grad(nl-1)=uk(nl-1)*uk(nl-1)+vk(nl-1)*vk(nl-1)
69 Time_step=1.d0
70 IF(Time_step.NE.0) THEN
71 R_r=0.01d0
72 taumin=tmax*R_r
73 taut=dt
74 YMAX=Y(nr)
75 Smax=SR(nr)
76 Smin=SR(nr)
77 Vmax=Dsqrt(GRAD(nr))
78 DO i = nr+1,Nl-1
79 Smax=DMAX1(Smax,SR(I))
80 Smin=DMIN1(Smin,SR(I))
81 Vmax=DMAX1(Vmax,DSQRT(GRAD(I)))
82 END DO
83 bet=Smin/Smax
84 IF(VMAX*dt.GT.R_R) THEN
85 dt=dt/1.25d0
86 END IF
87 IF(dT.GT.Tmax) dT=Tmax
88 IF(dT.LT.TMin) dT=Tmin
89 If(Dabs(dT-TauT).lt.1.d-4) Then
90 Kt=Kt+1
91 If(Kt.EQ.5) Then
92 Kt=0
93 dT=dT*1.25
94 End If
95 Else
96 Kt=0
97 End If
98 END IF
99 RETURN
100 END
101
102 Real*8 Function DDERIV(F,c,I)
103 Implicit REAL*8 (A-H,O-Z)
104 Dimension F(5)
105 Character*1 c
106 A=-1.d0/12.d0
107 IF(c.eq.'*') then
108 if(i.eq.0) DDERIV=A*(F(5)-8.d0*F(4)+8.d0*F(2)-F(1))
109 else
110 If(c.EQ.'+') Then
111 IF(I.eq.1) DDERIV=A*(3.d0*F(5)-
112 * 16.d0*F(4)+36.d0*F(3)-48.d0*F(2)+25.d0*F(1))
113 IF(I.eq.2) DDERIV=A*(3.d0*F(1)+
114 * 10.d0*F(2)-18.d0*F(3)+6.d0*F(4)-F(5))
115 Else
116 IF(I.eq.2) DDERIV=-A*(-F(1)+
117 * 6.d0*F(2)-18.d0*F(3)+10.d0*F(4)+3.d0*F(5))
118 IF(I.eq.1) DDERIV=-A*(3.d0*F(1)-
119 * 16.d0*F(2)+36.d0*F(3)-48.d0*F(4)+25.d0*F(5))
120 End if
121 End if
122 Return
123 End
Рис. 69. Листинг подпрограммы для реализации алгоритма вычисления скоростей (используется при решении нестационарной задачи)