mdln3plane.bas

用visual basic编写的有限元程序！

Attribute VB_Name = "mdlN3Plane"
Public ArrAK() As Double
Public ArrBK() As Double
Public ArrAR() As Double
Public ArrAQ() As Double

Function Solu()
    Rem 以下合成总刚矩阵
    If Top_Node = 0 Or Top_Elem = 0 Then Exit Function
    
    ReDim ArrAK(1 To 2 * Top_Node, 1 To 2 * Top_Node)
    ReDim ArrBK#(1 To 2 * Top_Node, 1 To 2 * Top_Node)
    ReDim ArrAR#(1 To 2 * Top_Node, 1 To 1)
    ReDim ArrAQ#(1 To 2 * Top_Node, 1 To 1)
    
    Dim I&, R%, S%
    
    Rem 计算材料特性矩阵
    For I = 1 To Top_Mati
        Call Mati_Matrix(Matis(I))
    Next
    
    Rem 合成总刚矩阵
    Dim R0&, S0&, R1&, S1&
    For I = 1 To Top_Elem
        Call Cal_Ke(Elems(I), Elem2s(I))  '单刚矩阵
        For R = 1 To 3: For S = 1 To 3
            R1 = 2 * Elems(I).Nodes(R): S1 = 2 * Elems(I).Nodes(S)
            R0 = R1 - 1: S0 = S1 - 1
            ArrAK(R0, S0) = ArrAK(R0, S0) + Elem2s(I).ArrK(2 * R - 1, 2 * S - 1)
            ArrAK(R0, S1) = ArrAK(R0, S1) + Elem2s(I).ArrK(2 * R - 1, 2 * S - 0)
            ArrAK(R1, S0) = ArrAK(R1, S0) + Elem2s(I).ArrK(2 * R - 0, 2 * S - 1)
            ArrAK(R1, S1) = ArrAK(R1, S1) + Elem2s(I).ArrK(2 * R - 0, 2 * S - 0)
        Next: Next
    Next
    
    Rem 总刚矩阵备份
    Call MOPER_MULT1(ArrBK, ArrAK, 1)
    
    Rem 以下添加载荷
    For I = 1 To Top_Node
        If Nodes(I).F_Code And (2 ^ CN_X) Then
            ArrAR(2 * I - 1, 1) = Nodes(I).F(CN_X)
        End If
        If Nodes(I).F_Code And (2 ^ CN_Y) Then
            ArrAR(2 * I, 1) = Nodes(I).F(CN_Y)
        End If
    Next
    
    Rem 以下添加约束条件
    For I = 1 To Top_Node
        If Nodes(I).U_Code > 0 Then
            Call Constrains(ArrBK, ArrAR, I, Nodes(I))
        End If
    Next
    ' MsgArr ArrBK
    
    Rem 计算
    Call MOPER_SOLV(ArrAQ, ArrBK, ArrAR)
    Call MOPER_MULT(ArrAR, ArrAK, ArrAQ)
End Function

Public Sub Post1()
    Dim I&, J&
    Rem 计算节点位移
    For I = 1 To Top_Node
        Node2s(I).U(0) = ArrAQ(2 * I - 1, 1)
        Node2s(I).U(1) = ArrAQ(2 * I, 1)
    Next
    Rem 计算节点载荷
    For I = 1 To Top_Node
        Node2s(I).R(0) = ArrAR(2 * I - 1, 1)
        Node2s(I).R(1) = ArrAR(2 * I, 1)
    Next
    'MsgArr ArrAR
End Sub

Function Constrains(AKe#(), ARe#(), N&, Node1 As Ans_Node)
    Rem 加位移约束
    Dim M&, I&, NX&, NY&, KX#, KY#, TEMPX#, TEMPY#
    M = UBound(AKe, 1): NX = 2 * N - 1: NY = 2 * N
    KX = Cos(Node1.Theta(0)): KY = Sin(Node1.Theta(0))
    If bBit(Node1.U_Code, CN_X) And bBit(Node1.U_Code, CN_Y) Then   '全部约束
        If Node1.Theta(0) <> 0 Then '把非零角度转为零角度
            TEMPX = Node1.U(CN_X) * KX - Node1.U(CN_Y) * KY
            TEMPY = Node1.U(CN_X) * KY + Node1.U(CN_Y) * KX
            Node1.U(CN_X) = TEMPX: Node1.U(CN_Y) = TEMPY
            Node1.Theta(0) = 0
        End If
        For I = 1 To M
            AKe(NX, I) = 0: AKe(NY, I) = 0
        Next
        AKe(NX, NX) = 1
        AKe(NY, NY) = 1
        ARe(NX, 1) = Node1.U(CN_X)
        ARe(NY, 1) = Node1.U(CN_Y)
    ElseIf bBit(Node1.U_Code, CN_X) Then 'UX方向约束
        For I = 1 To M
            AKe(NY, I) = AKe(NY, I) * KX - AKe(NX, I) * KY
            AKe(NX, I) = 0
        Next
        ARe(NY, 1) = ARe(NY, 1) * KX - ARe(NX, 1) * KY
        AKe(NX, NX) = KX
        AKe(NX, NY) = KY
        ARe(NX, 1) = Node1.U(CN_X)
    ElseIf bBit(Node1.U_Code, CN_Y) Then  'UY方向约束
        For I = 1 To M
            AKe(NX, I) = AKe(NX, I) * KX + AKe(NY, I) * KY
            AKe(NY, I) = 0
        Next
        ARe(NX, 1) = ARe(NX, 1) * KX + ARe(NY, 1) * KY
        AKe(NY, NX) = -KY
        AKe(NY, NY) = KX
        ARe(NY, 1) = Node1.U(CN_Y)
    End If
End Function

Function GetRi(ByVal N%)
    Rem 获结点合外力
    Dim I%
    Rx = 0: Ry = 0
    For I = 0 To MN - 1
        Rx = Rx + AKe.arr(2 * N, I) * AQe.arr(I, 0)
        Ry = Ry + AKe.arr(2 * N + 1, I) * AQe.arr(I, 0)
    Next
End Function

Private Sub Cal_Ke(Ea As Ans_Elem, Eb As Ans_Elem2)
    Dim I&, J&, K&
    Rem 取节点
    I = Ea.Nodes(1): J = Ea.Nodes(2): K = Ea.Nodes(3)
    
    Rem 1* 计算常数: ai,aj,ak,bi,bj,bk,ci,cj,ck
    Eb.a(1) = NX(J) * NY(K) - NX(K) * NY(J)
    Eb.a(2) = NX(K) * NY(I) - NX(I) * NY(K)
    Eb.a(3) = NX(I) * NY(J) - NX(J) * NY(I)
    Eb.b(1) = NY(J) - NY(K)
    Eb.b(2) = NY(K) - NY(I)
    Eb.b(3) = NY(I) - NY(J)
    Eb.C(1) = NX(K) - NX(J)
    Eb.C(2) = NX(I) - NX(K)
    Eb.C(3) = NX(J) - NX(I)
    
    Rem 2* 计算面积: Area
    Eb.Area = (Eb.a(1) + Eb.a(2) + Eb.a(3)) / 2#
    Dim R%, S%
    
    Rem 3* 填写几何矩阵: [B]
    For R = 1 To 3
        Eb.ArrB(1, 2 * R - 1) = Eb.b(R): Eb.ArrB(1, 2 * R) = 0
        Eb.ArrB(2, 2 * R - 1) = 0:       Eb.ArrB(2, 2 * R) = Eb.C(R)
        Eb.ArrB(3, 2 * R - 1) = Eb.C(R): Eb.ArrB(3, 2 * R) = Eb.b(R)
    Next
    Call MOPER_MULT1(Eb.ArrB, Eb.ArrB, 0.5 / Eb.Area)
    
    Eb.ArrD(1, 1) = Eb.E / (1 - Eb.U ^ 2): Eb.ArrD(1, 2) = Eb.E * Eb.U / (1 - Eb.U ^ 2)
    Eb.ArrD(2, 2) = Eb.E / (1 - Eb.U ^ 2): Eb.ArrD(2, 1) = Eb.ArrD(1, 2)
    Eb.ArrD(3, 3) = Eb.E / (1 + Eb.U ^ 2) / 2#
    
    Rem 4* 计算刚度矩阵: [K]
    Dim ArrBt#(1 To 6, 1 To 3), ArrR1#(1 To 3, 1 To 6)
    Call MOPER_TURN(ArrBt, Eb.ArrB)
    Call MOPER_MULT(ArrR1, Matis(Ea.MAT).ArrD, Eb.ArrB)
    Call MOPER_MULT(Eb.ArrK, ArrBt, ArrR1)
    Call MOPER_MULT1(Eb.ArrK, Eb.ArrK, Reals(Ea.REL).R(1) * Eb.Area)
End Sub