www.pudn.com > Numercial-heat-transfer-LBM.rar > Form1.frm, change:2012-04-20,size:12078b


VERSION 5.00 
Begin VB.Form Form1  
   Caption         =   "Form1" 
   ClientHeight    =   8025 
   ClientLeft      =   60 
   ClientTop       =   345 
   ClientWidth     =   12885 
   LinkTopic       =   "Form1" 
   ScaleHeight     =   8025 
   ScaleWidth      =   12885 
   StartUpPosition =   3  '窗口缺省 
   Begin VB.TextBox Text6  
      Height          =   495 
      Left            =   0 
      TabIndex        =   10 
      Text            =   "1" 
      Top             =   4560 
      Width           =   1215 
   End 
   Begin VB.TextBox Text5  
      Height          =   495 
      Left            =   0 
      TabIndex        =   9 
      Text            =   "1000" 
      Top             =   3240 
      Width           =   1215 
   End 
   Begin VB.TextBox Text4  
      Height          =   735 
      Left            =   0 
      TabIndex        =   6 
      Top             =   6000 
      Width           =   1575 
   End 
   Begin VB.TextBox Text3  
      Height          =   6615 
      Left            =   1920 
      MultiLine       =   -1  'True 
      TabIndex        =   5 
      Top             =   1080 
      Width           =   10695 
   End 
   Begin VB.TextBox Text2  
      Height          =   615 
      Left            =   0 
      TabIndex        =   2 
      Text            =   "10" 
      Top             =   1800 
      Width           =   1335 
   End 
   Begin VB.TextBox Text1  
      Height          =   615 
      Left            =   0 
      TabIndex        =   1 
      Text            =   "10" 
      Top             =   360 
      Width           =   1215 
   End 
   Begin VB.CommandButton Command1  
      Caption         =   "计算" 
      Height          =   735 
      Left            =   0 
      TabIndex        =   0 
      Top             =   6840 
      Width           =   1695 
   End 
   Begin VB.Label Label6  
      Caption         =   "Pr" 
      Height          =   255 
      Left            =   0 
      TabIndex        =   12 
      Top             =   4080 
      Width           =   975 
   End 
   Begin VB.Label Label5  
      Caption         =   "Re" 
      Height          =   255 
      Left            =   0 
      TabIndex        =   11 
      Top             =   2760 
      Width           =   975 
   End 
   Begin VB.Label Label4  
      Caption         =   "各点密度值" 
      Height          =   495 
      Left            =   3240 
      TabIndex        =   8 
      Top             =   240 
      Width           =   2175 
   End 
   Begin VB.Label Label3  
      Caption         =   "迭代次数" 
      Height          =   255 
      Left            =   0 
      TabIndex        =   7 
      Top             =   5520 
      Width           =   1335 
   End 
   Begin VB.Label Label2  
      Caption         =   "y方向格子数" 
      Height          =   375 
      Left            =   0 
      TabIndex        =   4 
      Top             =   1200 
      Width           =   1335 
   End 
   Begin VB.Label Label1  
      Caption         =   "x方向格子数" 
      Height          =   255 
      Left            =   0 
      TabIndex        =   3 
      Top             =   0 
      Width           =   1215 
   End 
End 
Attribute VB_Name = "Form1" 
Attribute VB_GlobalNameSpace = False 
Attribute VB_Creatable = False 
Attribute VB_PredeclaredId = True 
Attribute VB_Exposed = False 
Private Sub Command1_Click() 
Dim NI, NJ As Integer 
NI = Val(Text1.Text) 
NJ = Val(Text2.Text) 
ReDim feq(8, NI, NJ), f(8, NI, NJ), f_(8, NI, NJ), fnext(8, NI, NJ), geq(8, NI, NJ), g(8, NI, NJ), g_(8, NI, NJ), gnext(8, NI, NJ), u(NI, NJ), v(NI, NJ), w(NI, NJ), rou(NI, NJ), r(NI, NJ), T(NI, NJ) As Double 
Dim i, j, k, x, y As Integer 
Dim times As Long 
Dim Re, Pr, u0, T0, taof, rou0, diff, sumr, sumr0 As Double 
 
'赋初值 
times = 0 
u0 = 0.1 
T0 = 1 
Re = Val(Text5.Text) 
Pr = Val(Text6.Text) 
taof = 3 * u0 * NJ / Re + 0.5 
taog = (taof - 0.5) / Pr + 0.5 
sumr0 = 0 
sumr = 0 
For i = 0 To NI 
  For j = 0 To NJ 
    u(i, j) = 0 
    v(i, j) = 0 
    rou(i, j) = 1 
    T(i, j) = 0 
    sumr0 = sumr0 + rou(i, j) 
    If j = NJ And i <> 0 And i <> NI Then 
    u(i, NJ) = u0 
    T(i, NJ) = T0 
    End If 
    feq(0, i, j) = 4 / 9 * rou(i, j) * (1 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
    feq(1, i, j) = 1 / 9 * rou(i, j) * (1 + 3 * u(i, j) + 4.5 * u(i, j) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
    feq(2, i, j) = 1 / 9 * rou(i, j) * (1 + 3 * v(i, j) + 4.5 * v(i, j) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
    feq(3, i, j) = 1 / 9 * rou(i, j) * (1 - 3 * u(i, j) + 4.5 * u(i, j) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
    feq(4, i, j) = 1 / 9 * rou(i, j) * (1 - 3 * v(i, j) + 4.5 * v(i, j) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
    feq(5, i, j) = 1 / 36 * rou(i, j) * (1 + 3 * (u(i, j) + v(i, j)) + 4.5 * (u(i, j) + v(i, j)) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
    feq(6, i, j) = 1 / 36 * rou(i, j) * (1 + 3 * (-u(i, j) + v(i, j)) + 4.5 * (-u(i, j) + v(i, j)) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
    feq(7, i, j) = 1 / 36 * rou(i, j) * (1 - 3 * (u(i, j) + v(i, j)) + 4.5 * (u(i, j) + v(i, j)) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
    feq(8, i, j) = 1 / 36 * rou(i, j) * (1 + 3 * (u(i, j) - v(i, j)) + 4.5 * (u(i, j) - v(i, j)) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
    For k = 0 To 8 
    geq(k, i, j) = feq(k, i, j) * T(i, j) 
    f(k, i, j) = feq(k, i, j) 
    g(k, i, j) = geq(k, i, j) 
    Next 
  Next 
Next 
Do 
    times = times + 1 
    '内部碰撞 
    For i = 1 To NI - 1 
      For j = 1 To NJ - 1 
        For k = 0 To 8 
        f_(k, i, j) = f(k, i, j) - (f(k, i, j) - feq(k, i, j)) / taof 
        g_(k, i, j) = g(k, i, j) - (g(k, i, j) - geq(k, i, j)) / taog 
        Next 
      Next 
    Next 
    '边界碰撞 
    For j = 1 To NJ - 1                                    '左右 
      For k = 0 To 8 
       f_(k, 0, j) = feq(k, 0, j) + (1 - 1 / taof) * (f(k, 1, j) - feq(k, 1, j)) 
       f_(k, NI, j) = feq(k, NI, j) + (1 - 1 / taof) * (f(k, NI - 1, j) - feq(k, NI - 1, j)) 
       g_(k, 0, j) = geq(k, 0, j) + (1 - 1 / taog) * (g(k, 1, j) - geq(k, 1, j)) 
       g_(k, NI, j) = geq(k, NI, j) + (1 - 1 / taog) * (g(k, NI - 1, j) - geq(k, NI - 1, j)) 
      Next 
    Next 
    For i = 1 To NI - 1                                    '上下 
      For k = 0 To 8 
        f_(k, i, NJ) = feq(k, i, NJ) + (1 - 1 / taof) * (f(k, i, NJ - 1) - feq(k, i, NJ - 1)) 
        f_(k, i, 0) = feq(k, i, 0) + (1 - 1 / taof) * (f(k, i, 1) - feq(k, i, 1)) 
        g_(k, i, NJ) = geq(k, i, NJ) + (1 - 1 / taog) * (g(k, i, NJ - 1) - geq(k, i, NJ - 1)) 
        g_(k, i, 0) = geq(k, i, 0) + (1 - 1 / taog) * (g(k, i, 1) - geq(k, i, 1)) 
      Next 
    Next 
       For k = 0 To 8                                      '角 
         f_(k, 0, 0) = feq(k, 0, 0) + (1 - 1 / taof) * (f(k, 1, 1) - feq(k, 1, 1)) 
         f_(k, NI, 0) = feq(k, NI, 0) + (1 - 1 / taof) * (f(k, NI - 1, 1) - feq(k, NI - 1, 1)) 
         f_(k, 0, NJ) = feq(k, 0, NJ) + (1 - 1 / taof) * (f(k, 1, NJ - 1) - feq(k, 1, NJ - 1)) 
         f_(k, NI, NJ) = feq(k, NI, NJ) + (1 - 1 / taof) * (f(k, NI - 1, NJ - 1) - feq(k, NI - 1, NJ - 1)) 
         g_(k, 0, 0) = geq(k, 0, 0) + (1 - 1 / taog) * (g(k, 1, 1) - geq(k, 1, 1)) 
         g_(k, NI, 0) = geq(k, NI, 0) + (1 - 1 / taog) * (g(k, NI - 1, 1) - geq(k, NI - 1, 1)) 
         g_(k, 0, NJ) = geq(k, 0, NJ) + (1 - 1 / taog) * (g(k, 1, NJ - 1) - geq(k, 1, NJ - 1)) 
         g_(k, NI, NJ) = geq(k, NI, NJ) + (1 - 1 / taog) * (g(k, NI - 1, NJ - 1) - geq(k, NI - 1, NJ - 1)) 
      Next 
    '梳理 
    For i = 1 To NI - 1                                  '内部 
      For j = 1 To NJ - 1 
          fnext(0, i, j) = f_(0, i, j) 
          fnext(1, i, j) = f_(1, i - 1, j) 
          fnext(2, i, j) = f_(2, i, j - 1) 
          fnext(3, i, j) = f_(3, i + 1, j) 
          fnext(4, i, j) = f_(4, i, j + 1) 
          fnext(5, i, j) = f_(5, i - 1, j - 1) 
          fnext(6, i, j) = f_(6, i + 1, j - 1) 
          fnext(7, i, j) = f_(7, i + 1, j + 1) 
          fnext(8, i, j) = f_(8, i - 1, j + 1) 
          gnext(0, i, j) = g_(0, i, j) 
          gnext(1, i, j) = g_(1, i - 1, j) 
          gnext(2, i, j) = g_(2, i, j - 1) 
          gnext(3, i, j) = g_(3, i + 1, j) 
          gnext(4, i, j) = g_(4, i, j + 1) 
          gnext(5, i, j) = g_(5, i - 1, j - 1) 
          gnext(6, i, j) = g_(6, i + 1, j - 1) 
          gnext(7, i, j) = g_(7, i + 1, j + 1) 
          gnext(8, i, j) = g_(8, i - 1, j + 1) 
      Next 
    Next 
     
    '宏观参数 
    For i = 1 To NI - 1 
      For j = 1 To NJ - 1 
        rou(i, j) = fnext(0, i, j) + fnext(1, i, j) + fnext(2, i, j) + fnext(3, i, j) + fnext(4, i, j) + fnext(5, i, j) + fnext(6, i, j) + fnext(7, i, j) + fnext(8, i, j) 
        u(i, j) = (fnext(1, i, j) + fnext(5, i, j) + fnext(8, i, j) - fnext(3, i, j) - fnext(6, i, j) - fnext(7, i, j)) / rou(i, j) 
        v(i, j) = (fnext(2, i, j) + fnext(5, i, j) + fnext(6, i, j) - fnext(4, i, j) - fnext(7, i, j) - fnext(8, i, j)) / rou(i, j) 
        T(i, j) = (gnext(0, i, j) + gnext(1, i, j) + gnext(2, i, j) + gnext(3, i, j) + gnext(4, i, j) + gnext(5, i, j) + gnext(6, i, j) + gnext(7, i, j) + gnext(8, i, j)) / rou(i, j) 
        For k = 0 To 8 
          f(k, i, j) = fnext(k, i, j) 
          g(k, i, j) = gnext(k, i, j) 
        Next 
      Next 
    Next 
    For i = 1 To NI - 1 
      rou(i, 0) = rou(i, 1) 
      rou(i, NJ) = rou(i, NJ - 1) 
    Next 
    For j = 0 To NJ 
      rou(0, j) = rou(1, j) 
      rou(NI, j) = rou(NI - 1, j) 
      rou(NI, NJ) = 1 
    Next 
    For i = 0 To NI 
    '更新平衡态 
      For j = 0 To NJ 
        feq(0, i, j) = 4 / 9 * rou(i, j) * (1 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
        feq(1, i, j) = 1 / 9 * rou(i, j) * (1 + 3 * u(i, j) + 4.5 * u(i, j) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
        feq(2, i, j) = 1 / 9 * rou(i, j) * (1 + 3 * v(i, j) + 4.5 * v(i, j) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
        feq(3, i, j) = 1 / 9 * rou(i, j) * (1 - 3 * u(i, j) + 4.5 * u(i, j) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
        feq(4, i, j) = 1 / 9 * rou(i, j) * (1 - 3 * v(i, j) + 4.5 * v(i, j) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
        feq(5, i, j) = 1 / 36 * rou(i, j) * (1 + 3 * (u(i, j) + v(i, j)) + 4.5 * (u(i, j) + v(i, j)) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
        feq(6, i, j) = 1 / 36 * rou(i, j) * (1 + 3 * (-u(i, j) + v(i, j)) + 4.5 * (-u(i, j) + v(i, j)) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
        feq(7, i, j) = 1 / 36 * rou(i, j) * (1 - 3 * (u(i, j) + v(i, j)) + 4.5 * (u(i, j) + v(i, j)) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
        feq(8, i, j) = 1 / 36 * rou(i, j) * (1 + 3 * (u(i, j) - v(i, j)) + 4.5 * (u(i, j) - v(i, j)) ^ 2 - 1.5 * (u(i, j) ^ 2 + v(i, j) ^ 2)) 
        For k = 0 To 8 
        geq(k, i, j) = feq(k, i, j) * T(i, j) 
        Next 
        sumr = sumr + rou(i, j) 
      Next 
    Next 
diff = (Abs(sumr - sumr0)) / sumr 
sumr0 = sumr 
sumr = 0 
Loop Until diff < 10 ^ (-6) 
Text4.Text = Val(times) 
'输出 
'文本框输出各点密度 
For j = NJ To 0 Step -1 
  For i = 0 To NI 
  w(i, j) = (u(i, j) ^ 2 + v(i, j) ^ 2) ^ 0.5 
    Text3.Text = Text3.Text & Format(rou(i, j), "0.0000") & vbTab 
  Next 
    Text3.Text = Text3.Text & vbCrLf 
Next 
'观察速度 
Open App.Path & "\u.dat" For Output As #1 
    For j = 0 To NJ 
      For i = 0 To NI - 1 
        Print #1, Format(w(i, j), "0.0000") & vbTab; 
      Next 
        Print #1, Format(w(NI, j), "0.0000") 
        Print vbCrLf 
    Next 
Close #1 
'输出x轴1/2处ux分布 
Open App.Path & "\Ux.dat" For Output As #2 
  x = NI / 2 
  For j = 0 To NJ 
    Print #2, j & vbTab & Format(u(x, j), "0.0000") 
    Print vbCrLf 
  Next 
Close #2 
'输出y轴1/2处vy分布 
Open App.Path & "\Vx.dat" For Output As #3 
  y = NJ / 2 
  For i = 0 To NI 
    Print #3, i & vbTab & Format(v(i, y), "0.0000") 
    Print vbCrLf 
  Next 
Close #3 
'输出问题T 
Open App.Path & "\T.dat" For Output As #4 
   For j = 0 To NJ 
    For i = 0 To NI - 1 
      Print #4, Format(T(i, j), "0.0000") & vbTab; 
    Next 
      Print #4, Format(T(NI, j), "0.0000") 
      Print vbCrLf 
   Next 
Close #4 
     
 
End Sub