collide.vb

一个lbm的程序

Public Class Collide

    Friend Shared uv, eu, exf(3) As Double

    Friend Shared Function getfe0(ByVal i As Int32, ByVal j As Int32) As Double()
        Dim f(q) As Double
        Dim k As Byte
        uv = u(i, j) * u(i, j) + v(i, j) * v(i, j)
        For k = 0 To q
            If k = 0 Then
                f(k) = wi(0) * rho(i, j) * (1 - 1.5 * uv)
            Else
                eu = ei(k, 0) * u(i, j) + ei(k, 1) * v(i, j)
                f(k) = wi(k) * rho(i, j) * (1 + 3 * eu + 4.5 * eu * eu - 1.5 * uv)
            End If
        Next
        Return f
    End Function

    Friend Shared Function getfe(ByVal i As Int32, ByVal j As Int32) As Double()
        Dim ux, vx As Double
        Dim ft(q) As Double
        Dim k As Byte
        exf = getexfAt(i, j)
        ux = u(i, j) - (exf(0) + exf(2)) / rho(i, j)
        vx = v(i, j) - g / w / rho(i, j) - (exf(1) + exf(3)) / rho(i, j)
        uv = ux * ux + vx * vx
        For k = 0 To q
            If k = 0 Then
                ft(k) = wi(0) * rho(i, j) * (1 - 1.5 * uv)
            Else
                eu = ei(k, 0) * ux + ei(k, 1) * vx
                ft(k) = wi(k) * rho(i, j) * (1 + 3 * eu + 4.5 * eu * eu - 1.5 * uv)
            End If
        Next
        Return ft
    End Function

    Friend Shared Function getfe2(ByVal i As Int32, ByVal j As Int32, ByVal k As Byte) As Double
        uv = u(i, j) * u(i, j) + v(i, j) * v(i, j)
        eu = ei(k, 0) * u(i, j) + ei(k, 1) * v(i, j)
        Return wi(k) * ra * (1 + 3 * eu + 4.5 * eu * eu - 1.5 * uv)
    End Function

    Friend Shared Function getfe0(ByVal i As Int32, ByVal j As Int32, ByVal k As Byte) As Double
        Dim uv, eu As Double
        uv = u(i, j) * u(i, j) + v(i, j) * v(i, j)
        If k = 0 Then
            Return wi(0) * rho(i, j) * (1 - 1.5 * uv)
        Else
            eu = ei(k, 0) * u(i, j) + ei(k, 1) * v(i, j)
            Return wi(k) * rho(i, j) * (1 + 3 * eu + 4.5 * eu * eu - 1.5 * uv)
        End If
    End Function

    Shared Sub getmacro()
        Dim i, j As Int32
        Dim k As Byte
        Dim f() As Double
        For i = 0 To m
            For j = 0 To n
                If flag(i, j) >= CONST_INTERFACE Then
                    rho(i, j) = 0
                    u(i, j) = 0
                    v(i, j) = 0
                    f = field(i, j).f
                    For k = 0 To q
                        rho(i, j) += f(k)
                        u(i, j) += ei(k, 0) * f(k)
                        v(i, j) += ei(k, 1) * f(k)
                    Next
                    u(i, j) /= rho(i, j)
                    v(i, j) /= rho(i, j)
                End If
            Next
        Next
    End Sub

    Shared Sub collision()
        Dim i, j As Int32
        Dim k As Byte
        getmacro()
        Dim f(), feq() As Double
        For i = 0 To m
            For j = 0 To n
                If flag(i, j) >= CONST_INTERFACE Then
                    f = field(i, j).f
                    feq = getfe(i, j)
                    For k = 0 To q
                        f(k) = (1 - w) * f(k) + w * feq(k)
                    Next
                End If
            Next
        Next
    End Sub

    Friend Shared Function getexfAt(ByVal i As Int32, ByVal j As Int32) As Double()
        Dim k As Byte
        Dim ans(3) As Double
        Dim myv As Double = 0.001
        For k = 1 To q
            ans(0) += (1 - eps(i + ei(k, 0), j + ei(k, 1))) * myv * ei(k, 0)
            ans(1) += (1 - eps(i + ei(k, 0), j + ei(k, 1))) * myv * ei(k, 1)
            If flag(i + ei(k, 0), j + ei(k, 1)) = CONST_WALL Then
                ans(2) += gw(k) * ei(k, 0)
                ans(3) += gw(k) * ei(k, 1)
            End If
        Next
        ans(0) = beta * eps(i, j) * rho(i, j) * ans(0)
        ans(1) = beta * eps(i, j) * rho(i, j) * ans(1)
        ans(2) = eps(i, j) * rho(i, j) * ans(2)
        ans(3) = eps(i, j) * rho(i, j) * ans(3)
        Return ans
    End Function

End Class