form1.frm

一维交错格点法vb源程序, 取所有文本内容,然后以回车换行符进行字串分割.

VERSION 5.00
Begin VB.Form Form1 
   AutoRedraw      =   -1  'True
   Caption         =   "Form1"
   ClientHeight    =   6900
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   9495
   LinkTopic       =   "Form1"
   ScaleHeight     =   6900
   ScaleWidth      =   9495
   StartUpPosition =   3  '窗口缺省
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   495
      Left            =   7560
      TabIndex        =   0
      Top             =   3240
      Width           =   1335
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
'交错格点法，上下边界条件均为水位
'T为时间间隔，JP为每小时里的时层数，CN为糙率，M为计算时段的小时数，
'p为断面数，p1为间距数，
Const T = 60, JP = 60, P1 = 20, CN = 0.037, M = 36, P = 21

Dim Z0(P), Z1(P), Z2(P), W0(P), W1(P), W2(P), W3(P), ZP(2)

Dim U(2, P), Z(2, P), ZB(2, M), X(P1), R(P)
Dim II(P), JJ(P)
Dim Mat_A1(3, 3), Mat_U(P, 3), Mat_Z(P, 3), Mat_A2(3, 3), Mat_R(P)



Private Sub Command1_Click()
'定义两个时层的流量Q，水位Z，河宽B，截面面积A，边界断面每个小时节点的水位ZB，断面间距X
'定义断面形状的参数，Z0,Z1,Z2,W0,W1,W2,W3分别为断面中两个梯形的高和宽，ZP(2)是为了首末断面插值而设
IP = 1

For i = 1 To P
    For j = 1 To 3
        Mat_U(i, j) = 10
        Mat_Z(i, j) = 10
    Next j
    Mat_R(i) = 0
Next i

Mat_A1(1, 1) = 1 / 3: Mat_A1(1, 2) = 1 / 6: Mat_A1(1, 3) = 0
Mat_A1(2, 1) = 1 / 6: Mat_A1(2, 2) = 2 / 3: Mat_A1(2, 3) = 1 / 6
Mat_A1(3, 1) = 1 / 6: Mat_A1(3, 2) = 1 / 3: Mat_A1(3, 3) = 0
Mat_A2(1, 1) = -1 / 2: Mat_A2(1, 2) = 1 / 2: Mat_A2(1, 3) = 0
Mat_A2(2, 1) = -1 / 2: Mat_A2(2, 2) = 0: Mat_A2(2, 3) = 1 / 2
Mat_A2(3, 1) = -1 / 2: Mat_A2(3, 2) = 1 / 2: Mat_A2(3, 3) = 0
Open "in.txt" For Input As #1
Open "out1.txt" For Output As #2
'Open "out2.txt" For Output As #3
'首末断面第m个小时的水位
For i = 1 To 2
    For j = 1 To M
        Input #1, ZB(i, j)
    Next j
Next i

'令初始速度0，如果格式稳定，误差将在计算中消除
For i = 1 To P
    U(1, i) = 5
    U(0, i) = 5
Next i
'内插每个断面初始水位


'间距
For i = 1 To P1
    Input #1, X(i)
Next i
For i = 1 To P1
    X(i) = 1000
Next i
'断面形状参数
For i = 1 To P
    Input #1, Z0(i), Z1(i), Z2(i)
  '  Print Z0(i)
Next i

'For j = 1 To M
'    ZB(1, j) = ZB(1, j) - Z0(1)
'    ZB(2, j) = ZB(2, j) - Z0(21)
'Next j
For i = 1 To P
    Z(1, i) = ZB(1, 1) + (ZB(2, 1) - ZB(1, 1)) * (i - 1) / P1
    Z(0, i) = Z(1, i)
    'Print #2, Z(1, i);
Next i
For i = 1 To P
    Input #1, W0(i), W1(i), W2(i), W3(i)
Next

For j = 1 To P
   ' II(j) = (Z0(j) - Z0(j + 1)) / 2
   II(j) = 0
   ' Print II(j)
Next j
'II(1) = II(1) / 2
'II(P) = II(P1) / 2



'Dim it As Integer
'计算各个断面面积和宽度
 'For I = 1 To P
'    it = I
'    Call F_and_B(A(1, I), B(1, I), Z(1, I), it)
'Next I
'Print #2, "断面";
For i = 1 To 21
    Print #2, Tab(10 * i); "#" & (i);
Next i
Print #2,
Dim flag1, flag2, n
flag2 = 1
n = 1
For i = 1 To 2
    ZP(i) = ZB(i, 1)
Next i
Do While n <= JP * M - JP - 1
    For i = 1 To P
        For j = 1 To 3
            Mat_U(i, j) = 0
            Mat_Z(i, j) = 0
        Next j
        Mat_R(i) = 0
    Next i
    If flag2 < 0 Then
        flag2 = 1
        For j = 1 To P
            Z(0, j) = Z(1, j)
            U(0, j) = U(1, j)
        Next j
        n = n + 1
        If n Mod JP = 0 Then Print n / JP, "hour"
    
        '内插该时层的水位值，因为边界条件只有每个小时的水位值，而时间间隔为60s
        For i = 1 To 2
            ZP(i) = ZB(i, IP) + (n - (IP - 1) * JP) / JP * (ZB(i, IP + 1) - ZB(i, IP))
            'Print ZP(1)
        Next i
        If n Mod JP = 0 Then IP = IP + 1
    
    End If
   Dim r1 As Double
  ' For i = 1 To P
  '      it = i
  '      Call SubR(R(i), Z(1, i), it)
        'Print #2, R(i); Tab(30); Z(1, i)
   'Next i
  ' Print #2, "------"
  'For i = 1 To P
  '      r1 = R(i)
  '      r1 = r1 ^ (4 / 3)
   '     JJ(i) = CN * CN * U(1, i) * U(1, i) / r1 'Sgn(U(1, i)) *
        'Print #2, JJ(i)
  'Next i
   
        
    
   'Print #2, R(18)
    'Dim it As Integer, R As Single
     '计算各个断面面积和宽度
        
    Call Matrix(U, Mat_U)
    
    Call Matrix(Z, Mat_Z)
 ' For i = 1 To P
 '       For j = 1 To 3
 '       Print #2, Mat_U(i, j),
 '       Next j
 '       Print #2,
 '  Next i
 '  Print #2, "---------"
   flag1 = 1
    Call Val_of_R(flag1, flag2, Mat_R)
  '  For i = 1 To P
  '    Print #2, Mat_R(i)
  '  Next i
  '  Print #2, "-------"
    
   Call Guass(Z, flag1, flag2, Mat_R)
  ' For i = 1 To P
  '    Print #2, Z(2, i)
  '  Next i
   'Print #2, "-------"
    
    flag1 = 2
    Call Val_of_R(flag1, flag2, Mat_R)
    'For i = 1 To P
    '   Print #2, Mat_R(i)
    'Next i
    'Print #2, "-------"
      
  '
    Call Guass(U, flag1, flag2, Mat_R)
    For i = 1 To P
      Print #2, U(2, i)
   Next i
   Print #2, "-------"
    For j = 1 To P
        
        
        Z(1, j) = Z(2, j)
        U(1, j) = U(2, j)
        'Print #2, Z(1, j)
    Next j
   
        
    flag2 = flag2 - 1
If n Mod JP = 0 And flag2 < 0 Then
    For i = 1 To P
        Print #2, Tab(10 * i); Format(Z(2, i) + Z0(i), "##.##");
       'Print #2, Tab(10 * i); Format(U(2, i), "##.##");
    Next i
    Print #2,
End If
    
        
       
Loop
        
 
Close #1
Close #2

        
        
    
End Sub

Private Sub SubR(R, Z, it)

Dim a, b, c, d, zh
zh = Z - Z0(it)
If zh <= Z1(it) - Z0(it) Then
'If Z <= Z1(it) - Z0(it) Then
   b = W0(it) + (W1(it) - W0(it)) / (Z1(it) - Z0(it)) * zh
   a = (W0(it) + b) / 2 * zh
   
Else
   b = W2(it) + (W3(it) - W2(it)) / (Z2(it) - Z1(it)) * (zh - Z1(it) + Z0(it))
   a = (W2(it) + b) / 2 * (zh - Z1(it) + Z0(it)) + (W0(it) + W1(it)) / 2 * (Z1(it) - Z0(it))
   
End If

R = a / b
End Sub

Private Sub Matrix(Xa, Y)



Dim X1, X2
X1 = Xa(1, 1) + 0.5 * Xa(1, 2) / 3
X2 = 0.5 * Xa(1, 1) + Xa(1, 2) / 3
Y(1, 1) = -X1
Y(1, 2) = X1
Y(2, 1) = -X2
Y(2, 2) = X2
For i = 2 To P1
    X1 = (Xa(1, i) + 0.5 * Xa(1, i + 1)) / 3
    X2 = (0.5 * Xa(1, i) + Xa(1, i + 1)) / 3
    Y(i, 2) = Y(i, 2) - X1
    Y(i, 3) = Y(i, 3) + X1
    Y(i + 1, 1) = Y(i + 1, 1) - X2
    Y(i + 1, 2) = Y(i + 1, 2) + X2
Next i
'Print Y(2, 1)
End Sub
Private Sub Val_of_R(flag1, flag2, Mat_R1)
Dim I1, dt, dt1
For i = 1 To P

    Mat_R1(i) = 0
Next i

If flag2 = 1 Then dt = T / 2
If flag2 = 0 Then dt = T

If flag1 = 1 Then
    For i = 1 To 2
        I1 = i + 19
        Mat_R1(1) = Mat_R1(1) + Mat_A1(1, i) * Z(0, i) - dt / X(1) * (Mat_U(1, i) * Z(1, i) + Mat_Z(1, i) * U(1, i))
        Mat_R1(21) = Mat_R1(21) + Mat_A1(3, i) * Z(0, I1) - dt / X(20) * (Mat_U(21, i) * Z(1, I1) + Mat_Z(21, i) * U(1, I1))
    Next i
    
    For i = 2 To P1
       dt1 = dt / X(i)
        For j = 1 To 3
            Mat_R1(i) = Mat_R1(i) + Mat_A1(2, j) * Z(0, i + j - 2) - dt1 * (Mat_U(i, j) * Z(1, i + j - 2) + Mat_Z(i, j) * U(1, i + j - 2))
        Next j
    Next i
ElseIf flag1 = 2 Then
    For i = 1 To 2
        I1 = i + 19
        Mat_R1(1) = Mat_R1(1) + Mat_A1(1, i) * U(0, i) - dt / X(1) * (9.8 * Mat_A2(1, i) * Z(1, i) + Mat_U(1, i) * U(1, i)) - 0.5 * dt * 9.8 * (II(i) - JJ(i))
        Mat_R1(21) = Mat_R1(1) + Mat_A1(3, i) * U(0, I1) - dt / X(20) * (9.8 * Mat_A2(3, i) * Z(1, I1) + Mat_U(21, i) * U(1, I1)) - 0.5 * dt * 9.8 * (II(I1) - JJ(I1))
    Next i
    
    For i = 2 To P1
        dt1 = dt / X(i)
        For j = 1 To 3
            Mat_R1(i) = Mat_R1(i) + Mat_A1(2, j) * U(0, i + j - 2) - dt1 * (9.8 * Mat_A2(2, j) * Z(1, i + j - 2) + Mat_U(i, j) * U(1, i + j - 2)) - 0.5 * dt * 9.8 * (II(i) - JJ(i))
        Next j
    Next i
End If


End Sub
Private Sub Guass(U_or_Z, flag1, flag2, Mat_R)
Dim GA(P), GR(P), Q

If flag2 = 1 And flag1 = 1 Then
    
   Mat_R(1) = ZP(1) * Mat_A1(1, 1) * 10 ^ 10
    Mat_R(P) = ZP(2) * Mat_A1(3, 2) * 10 ^ 10
End If
GA(1) = Mat_A1(1, 1)
If flag2 = 1 And flag1 = 1 Then GA(1) = GA(1) * 10 ^ 10
'GA(1) = GA(1) * 10 ^ 10
GR(1) = Mat_R(1)
GA(2) = Mat_A1(2, 2) - Mat_A1(2, 1) * Mat_A1(1, 2) / GA(1)
GR(2) = Mat_R(2) - Mat_A1(2, 1) * Mat_R(1) / GA(1)
For i = 3 To P1
    Q = Mat_A1(2, 1) / GA(i - 1)
    GA(i) = Mat_A1(2, 2) - Q * Mat_A1(2, 3)
    GR(i) = Mat_R(i) - Q * GR(i - 1)
Next i
GA(P) = Mat_A1(3, 2) - Mat_A1(3, 1) * Mat_A1(2, 3) / GA(P1)
If flag2 = 1 And flag1 = 1 Then GA(P) = Mat_A1(3, 2) * 10 ^ 10 - Mat_A1(3, 1) * Mat_A1(2, 3) / GA(P1)
'GA(P) = Mat_A1(3, 2) * 10 ^ 10 - Mat_A1(3, 1) * Mat_A1(2, 3) / GA(P1)
GR(P) = Mat_R(P) - Mat_A1(3, 1) * GR(P1) / GA(P1)
U_or_Z(2, P) = GR(P) / GA(P)
For i = 20 To 2 Step -1
    U_or_Z(2, i) = GR(i) / GA(i) - Mat_A1(2, 3) * U_or_Z(2, i + 1) / GA(i)
Next i
U_or_Z(2, 1) = GR(1) / GA(1) - Mat_A1(1, 2) * U_or_Z(2, 2) / GA(1)
'If flag1 = 1 Then
 
 '   U_or_Z = z
'ElseIf flag1 = 2 Then
'    U_or_Z = U
'End If
'Print U_or_Z(2, 20)
'Print GR(P) / GA(P)
'Print ZP(2)
End Sub