📄 vbsimple.bas
字号:
Call DIFLOW
AIM(i + 1, j) = ACOF + Max(0#, FLOW)
AIP(i, j) = AIM(i + 1, j) - FLOW
If j <> M2 Then
FL = XCVI(i) * f(i, j + 1, 2) * (FY(j + 1) * RHO(i, j + 1) + FYM(j + 1) * RHO(i, j))
FLM = XCVIP(i - 1) * f(i - 1, j + 1, 2) * (FY(j + 1) * RHO(i - 1, j + 1) + FYM(j + 1) * _
RHO(i - 1, j))
GM = GAM(i, j) * GAM(i, j + 1) / (YCV(j) * GAM(i, j + 1) + YCV(j + 1) * GAM(i, j) + _
1E-30) * XCVI(i)
GMM = GAM(i - 1, j) * GAM(i - 1, j + 1) / (YCV(j) * GAM(i - 1, j + 1) + YCV(j + 1) * _
GAM(i - 1, j) + 1E-30) * XCVIP(i - 1)
DIFF = RMN(j + 1) * 2# * (GM + GMM)
Else
FL = XCVI(i) * f(i, M1, 2) * RHO(i, M1)
FLM = XCVIP(i - 1) * f(i - 1, M1, 2) * RHO(i - 1, M1)
DIFF = R(M1) * (XCVI(i) * GAM(i, M1) + XCVIP(i - 1) * GAM(i - 1, M1)) / YDIF(M1)
End If
FLOW = RMN(j + 1) * (FL + FLM)
Call DIFLOW
AJM(i, j + 1) = ACOF + Max(0#, FLOW)
AJP(i, j) = AJM(i, j + 1) - FLOW
VOL = YCVR(j) * XCVS(i)
APT = (RHO(i, j) * XCVI(i) + RHO(i - 1, j) * XCVIP(i - 1)) _
/ (XCVS(i) * DT)
AP(i, j) = AP(i, j) - APT
CON(i, j) = CON(i, j) + APT * f(i, j, 1)
AP(i, j) = (-AP(i, j) * VOL + AIP(i, j) + AIM(i, j) + AJP(i, j) + AJM(i, j)) _
/ RELAX(NF)
CON(i, j) = CON(i, j) * VOL + REL * AP(i, j) * f(i, j, 1)
DU(i, j) = VOL / (XDIF(i) * SX(j))
CON(i, j) = CON(i, j) + DU(i, j) * (P(i - 1, j) - P(i, j))
DU(i, j) = DU(i, j) / AP(i, j)
Next i
Next j
Call Solve
End If
'---COEFFICIENTS FOR THE V EQUATION----
NF = 2
If LSolve(NF) Then
IST = 2
JST = 3
Call GAMSOR
REL = 1# - RELAX(NF)
For i = 2 To L2
AREA = R(1) * XCV(i)
FLOW = AREA * f(i, 2, 2) * RHO(i, 1)
DIFF = AREA * GAM(i, 1) / YCV(2)
Call DIFLOW
AJM(i, 3) = ACOF + Max(0#, FLOW)
Next i
For j = 3 To M2
FL = ARXJ(j) * f(2, j, 1) * RHO(1, j)
FLM = ARXJP(j - 1) * f(2, j - 1, 1) * RHO(1, j - 1)
FLOW = FL + FLM
DIFF = (ARXJ(j) * GAM(1, j) + ARXJP(j - 1) * GAM(1, j - 1)) / (XDIF(2) * SXMN(j))
Call DIFLOW
AIM(2, j) = ACOF + Max(0#, FLOW)
For i = 2 To L2
If i <> L2 Then
FL = ARXJ(j) * f(i + 1, j, 1) * (FX(i + 1) * RHO(i + 1, j) + FXM(i + 1) * RHO(i, j))
FLM = ARXJP(j - 1) * f(i + 1, j - 1, 1) * (FX(i + 1) * RHO(i + 1, j - 1) + FXM(i + 1) * _
RHO(i, j - 1))
GM = GAM(i, j) * GAM(i + 1, j) / (XCV(i) * GAM(i + 1, j) + XCV(i + 1) * GAM(i, j) + _
1E-30) * ARXJ(j)
GMM = GAM(i, j - 1) * GAM(i + 1, j - 1) / (XCV(i) * GAM(i + 1, j - 1) + XCV(i + 1) * _
GAM(i, j - 1) + 1E-30) * ARXJP(j - 1)
DIFF = 2# * (GM + GMM) / SXMN(j)
Else
FL = ARXJ(j) * f(L1, j, 1) * RHO(L1, j)
FLM = ARXJP(j - 1) * f(L1, j - 1, 1) * RHO(L1, j - 1)
DIFF = (ARXJ(j) * GAM(L1, j) + ARXJP(j - 1) * GAM(L1, j - 1)) / (XDIF(L1) * SXMN(j))
End If
FLOW = FL + FLM
Call DIFLOW
AIM(i + 1, j) = ACOF + Max(0#, FLOW)
AIP(i, j) = AIM(i + 1, j) - FLOW
If j <> M2 Then
AREA = R(j) * XCV(i)
FL = f(i, j, 2) * (FY(j) * RHO(i, j) + FYM(j) * RHO(i, j - 1)) * RMN(j)
FLP = f(i, j + 1, 2) * (FY(j + 1) * RHO(i, j + 1) + FYM(j + 1) * RHO(i, j)) * RMN(j + 1)
FLOW = (FV(j) * FL + FVP(j) * FLP) * XCV(i)
DIFF = AREA * GAM(i, j) / YCV(j)
Else
AREA = R(M1) * XCV(i)
FLOW = AREA * f(i, M1, 2) * RHO(i, M1)
DIFF = AREA * GAM(i, M1) / YCV(M2)
End If
Call DIFLOW
AJM(i, j + 1) = ACOF + Max(0#, FLOW)
AJP(i, j) = AJM(i, j + 1) - FLOW
VOL = YCVRS(j) * XCV(i)
SXT = SX(j)
If j = M2 Then SXT = SX(M1)
SXB = SX(j - 1)
If j = 3 Then SXB = SX(1)
APT = (ARXJ(j) * RHO(i, j) * 0.5 * (SXT + SXMN(j)) + ARXJP(j - 1) * RHO(i, j - 1) * _
0.5 * (SXB + SXMN(j))) / (YCVRS(j) * DT)
AP(i, j) = AP(i, j) - APT
CON(i, j) = CON(i, j) + APT * f(i, j, 2)
AP(i, j) = (-AP(i, j) * VOL + AIP(i, j) + AIM(i, j) + AJP(i, j) + AJM(i, j)) _
/ RELAX(NF)
CON(i, j) = CON(i, j) * VOL + REL * AP(i, j) * f(i, j, 2)
DV(i, j) = VOL / YDIF(j)
CON(i, j) = CON(i, j) + DV(i, j) * (P(i, j - 1) - P(i, j))
DV(i, j) = DV(i, j) / AP(i, j)
Next i
Next j
Call Solve
End If
'---COEFIICIENTS FOR THE PRESSURE CORRECTION EQUATION----
NF = 3
If LSolve(NF) Then
IST = 2
JST = 2
Call GAMSOR
SMAX = 0#
SSUM = 0#
For j = 2 To M2
For i = 2 To L2
VOL = YCVR(j) * XCV(i)
CON(i, j) = CON(i, j) * VOL
Next i
Next j
For i = 2 To L2
ARHO = R(1) * XCV(i) * RHO(i, 1)
CON(i, 2) = CON(i, 2) + ARHO * f(i, 2, 2)
AJM(i, 2) = 0#
Next i
For j = 2 To M2
ARHO = ARX(j) * RHO(1, j)
CON(2, j) = CON(2, j) + ARHO * f(2, j, 1)
AIM(2, j) = 0#
For i = 2 To L2
If i <> L2 Then
ARHO = ARX(j) * (FX(i + 1) * RHO(i + 1, j) + FXM(i + 1) * RHO(i, j))
FLOW = ARHO * f(i + 1, j, 1)
CON(i, j) = CON(i, j) - FLOW
CON(i + 1, j) = CON(i + 1, j) + FLOW
AIP(i, j) = ARHO * DU(i + 1, j)
AIM(i + 1, j) = AIP(i, j)
Else
ARHO = ARX(j) * RHO(L1, j)
CON(i, j) = CON(i, j) - ARHO * f(L1, j, 1)
AIP(i, j) = 0#
End If
If j <> M2 Then
ARHO = RMN(j + 1) * XCV(i) * (FY(j + 1) * RHO(i, j + 1) + FYM(j + 1) * RHO(i, j))
FLOW = ARHO * f(i, j + 1, 2)
CON(i, j) = CON(i, j) - FLOW
CON(i, j + 1) = CON(i, j + 1) + FLOW
AJP(i, j) = ARHO * DV(i, j + 1)
AJM(i, j + 1) = AJP(i, j)
Else
ARHO = RMN(M1) * XCV(i) * RHO(i, M1)
CON(i, j) = CON(i, j) - ARHO * f(i, M1, 2)
AJP(i, j) = 0#
End If
AP(i, j) = AIP(i, j) + AIM(i, j) + AJP(i, j) + AJM(i, j)
f(i, j, 3) = 0#
SMAX = Max(SMAX, Abs(CON(i, j)))
SSUM = SSUM + CON(i, j)
Next i
Next j
Call Solve
'---COMEE HERE TO CORRECT THE PRESSURE AND VELOCITIES
For j = 2 To M2
For i = 2 To L2
P(i, j) = P(i, j) + f(i, j, 3) * RELAX(NP)
If i <> 2 Then f(i, j, 1) = f(i, j, 1) + DU(i, j) * (f(i - 1, j, 3) - f(i, j, 3))
If j <> 2 Then f(i, j, 2) = f(i, j, 2) + DV(i, j) * (f(i, j - 1, 3) - f(i, j, 3))
Next i
Next j
End If
'---COEFFICIENTS FOR OTHER EQUATIONS----
IST = 2
JST = 2
For N = 4 To NFMAX
NF = N
If LSolve(NF) Then
Call GAMSOR
REL = 1# - RELAX(NF)
For i = 2 To L2
AREA = R(1) * XCV(i): FLOW = AREA * f(i, 2, 2) * RHO(i, 1)
DIFF = AREA * GAM(i, 1) / YDIF(2)
Call DIFLOW
AJM(i, 2) = ACOF + Max(0#, FLOW)
Next i
For j = 2 To M2
FLOW = ARX(j) * f(2, j, 1) * RHO(1, j): DIFF = ARX(j) * GAM(1, j) / (XDIF(2) * SX(j))
Call DIFLOW
AIM(2, j) = ACOF + Max(0#, FLOW)
For i = 2 To L2
If i <> L2 Then
FLOW = ARX(j) * f(i + 1, j, 1) * (FX(i + 1) * RHO(i + 1, j) + FXM(i + 1) * RHO(i, j))
DIFF = ARX(j) * 2# * GAM(i, j) * GAM(i + 1, j) / ((XCV(i) * GAM(i + 1, j) + _
XCV(i + 1) * GAM(i, j) + 1E-30) * SX(j))
Else
FLOW = ARX(j) * f(L1, j, 1) * RHO(L1, j): DIFF = ARX(j) * GAM(L1, j) / (XDIF(L1) * SX(j))
End If
Call DIFLOW
AIM(i + 1, j) = ACOF + Max(0#, FLOW)
AIP(i, j) = AIM(i + 1, j) - FLOW: AREA = RMN(j + 1) * XCV(i)
If j <> M2 Then
FLOW = AREA * f(i, j + 1, 2) * (FY(j + 1) * RHO(i, j + 1) + FYM(j + 1) * RHO(i, j))
DIFF = AREA * 2# * GAM(i, j) * GAM(i, j + 1) / (YCV(j) * GAM(i, j + 1) + _
YCV(j + 1) * GAM(i, j) + 1E-30)
Else
FLOW = AREA * f(i, M1, 2) * RHO(i, M1): DIFF = AREA * GAM(i, M1) / YDIF(M1)
End If
Call DIFLOW
AJM(i, j + 1) = ACOF + Max(0#, FLOW): AJP(i, j) = AJM(i, j + 1) - FLOW
VOL = YCVR(j) * XCV(i): APT = RHO(i, j) / DT
AP(i, j) = AP(i, j) - APT
CON(i, j) = CON(i, j) + APT * f(i, j, NF)
AP(i, j) = (-AP(i, j) * VOL + AIP(i, j) + AIM(i, j) + AJP(i, j) + AJM(i, j)) _
/ RELAX(NF)
CON(i, j) = CON(i, j) * VOL + REL * AP(i, j) * f(i, j, NF)
Next i
Next j
Call Solve
End If
Next N
TIME = TIME + DT
ITER = ITER + 1
If ITER >= LAST Then LStop = True
End Sub
Sub Supply()
' 10 FORMAT(1X,26(1H*),3X,A10,3X,26(1H*))
' 20 FORMAT(1X,4H I =,I6,6I9)
' 30 FORMAT(1X,1HJ)
' 40 FORMAT(1X,I2,3X,1P7E9.2)
' 50 FORMAT(1X,1H )
' 51 FORMAT(1X,' I =',2X,7(I4,5X))
' 52 FORMAT(1X,' X =',1P7E9.2)
' 53 FORMAT(1X,'TH =',1P7E9.2)
' 54 FORMAT(1X,'J =',2X,7(I4,5X))
' 55 FORMAT(1X,'Y =',1P7E9.2)
End Sub
'************************************************************************
Sub UGRID()
XU(2) = 0#
DX = XL / CDbl(L1 - 2)
For i = 3 To L1
XU(i) = XU(i - 1) + DX
Next i
YV(2) = 0#
DY = YL / CDbl(M1 - 2)
For j = 3 To M1
YV(j) = YV(j - 1) + DY
Next
End Sub
'************************************************************************
Sub ToPrint()
If LPrint(3) Then
'---CALCULATE THE STREAM FUNTION----------------------------------------
f(2, 2, 3) = 0#
For i = 2 To L1
If i <> 2 Then f(i, 2, 3) = f(i - 1, 2, 3) - RHO(i - 1, 1) * f(i - 1, 2, 2) * R(1) * XCV(i - 1)
For j = 3 To M1
RHOM = FX(i) * RHO(i, j - 1) + FXM(i) * RHO(i - 1, j - 1)
f(i, j, 3) = f(i, j - 1, 3) + RHOM * f(i, j - 1, 1) * ARX(j - 1)
Next j
Next i
End If
'*
If LPrint(NP) Then
'*
'---CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATION
For j = 2 To M2
P(1, j) = (P(2, j) * XCVS(3) - P(3, j) * XDIF(2)) / XDIF(3)
P(L1, j) = (P(L2, j) * XCVS(L2) - P(L3, j) * XDIF(L1)) / XDIF(L2)
Next j
For i = 2 To L2
P(i, 1) = (P(i, 2) * YCVS(3) - P(i, 3) * YDIF(2)) / YDIF(3)
P(i, M1) = (P(i, M2) * YCVS(M2) - P(i, M3) * YDIF(M1)) / YDIF(M2)
Next i
P(1, 1) = P(2, 1) + P(1, 2) - P(2, 2)
P(L1, 1) = P(L2, 1) + P(L1, 2) - P(L2, 2)
P(1, M1) = P(2, M1) + P(1, M2) - P(2, M2)
P(L1, M1) = P(L2, M1) + P(L1, M2) - P(L2, M2)
PREF = P(IPREF, JPREF)
For j = 1 To M1
For i = 1 To L1
P(i, j) = P(i, j) - PREF
Next i
Next j
End If
If TIME <= (0.5 * DT) Then
Print #1, " ",
IEND = 0
Do While (IEND <> L1)
IBEG = IEND + 1
IEND = IEND + 7
IEND = Min(IEND, L1)
Print #1, " ",
For i = IBEG To IEND
Print #1, " I=", Format(i, "0000");
Next
Print #1,
If MODE <> 3 Then
For i = IBEG To IEND
Print #1, "X=", Format(X(i), "0.00E+000");
Next
Print #1,
Else
For i = IBEG To IEND
Print #1, "TH=", Format(X(i), "0.00E+000");
Next
Print #1,
End If
Loop
JEND = 0
Print #1, " ",
Do While (JEND <> M1)
JBEG = JEND + 1
JEND = JEND + 7
JEND = Min(JEND, M1)
Print #1, " ",
For j = JBEG To JEND
Print #1, "j=", Format(j, "0000");
Next
Print #1,
For j = JBEG To JEND
Print #1, "Y=", Format(Y(j), "0.00E+00");
Next
Print #1,
Loop
End If
For N = 1 To NGAM
NF = N
If LPrint(NF) Then
Print #1, " ",
Print #1, String(10, "*") + Title(NF) + String(10, "*")
IFST = 1
JFST = 1
If NF = 1 Or NF = 3 Then IFST = 2
If NF = 2 Or NF = 3 Then JFST = 2
IBEG = IFST - 7
Do
IBEG = IBEG + 7
IEND = IBEG + 6
IEND = Min(IEND, L1)
Print #1, " "
Print #1, "I=",
For i = IBEG To IEND
Print #1, Format(i, "##"); Space(10);
Next
Print #1,
Print #1, "J"
JFL = JFST + M1
For jj = JFST To M1
j = JFL - jj
Print #1, j,
For i = IBEG To IEND
Print #1, Format(f(i, j, NF), "0.00E+00"); Space(2);
Next
Print #1,
Next jj
Loop While (IEND < L1)
End If
Next N
End Sub
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -