生成刚度方程_解线性方程组.bas

钢筋混凝土框架结构设计计算源程序

Attribute VB_Name = "Module5"
Option Explicit
DefInt I-N
DefSng X-Y
Dim i As Integer, j As Integer, k As Integer
Type 杆端内力                        '杆端内力Member innner force自定义数据类型，存放单元的杆端力/固端力
    N1 As Single
    Q1 As Single
    M1 As Single
    N2 As Single
    Q2 As Single
    M2 As Single
End Type
'杆端内力(i) 是一个数组，其中每一元素包含如上6个数
'存放一个单元的杆端力,i即单元号ngle
'方程编号数组 Neq(结点号,位移号1/2/3)  Nequ-方程总数
Public 杆端力(1000) As 杆端内力
Public 固端力 As 杆端内力
Public Neq(1500, 3), Nequ          '方程编号数组 Neq(结点号,位移号1/2/3)  Nequ-方程总数
Public D(2000) As Double
Public A(1000, 1000) As Double
'D() 荷载列向量  A( , ) 结构刚度矩阵
Dim X1, Y1, x2, y2

Sub 方程编号()
    Dim i, j, n
    For i = 1 To Nodg                '假定为铰结点
        Neq(i, 1) = 1: Neq(i, 2) = 1: Neq(i, 3) = 0
    Next i
    For i = 1 To Ncell         '不是铰结点，是刚结点
        If NGN(i, 1) > 0 Then Neq(NGN(i, 1), 3) = 1
        If NGN(i, 2) > 0 Then Neq(NGN(i, 2), 3) = 1
    Next i
    For i = 1 To Nodg         '支座被约束方向，方程编号=0
        For j = 1 To 3
            If XYM(i, j) = "0" Then
                Neq(i, j) = 0
            End If
        Next j
    Next i
    n = 1
    For i = 1 To Nodg      '逐结点、同一结点xyφ方向依次编号
        For j = 1 To 3
            If Neq(i, j) = 1 Then
                Neq(i, j) = n
                n = n + 1
            End If
        Next j
    Next i
    Nequ = n - 1           '方程总数
End Sub
Sub 形成荷载向量()
    Dim fg(6) As Single, NV(6)
    Dim Ne, L As Single, sinA As Single, cosA As Single
    Dim Nd, i, n
    Erase D()
    For Ne = 1 To Ncell
        If Ld(Ne).Pq Like "[qptQPT]" Then
            '计算单元Ne 的长度L、倾角的Sin & Cos值，端点座标
            Call SetT(Ne, L, sinA, cosA, X1, Y1, x2, y2)
            Call 求固端力(Ld(Ne), L, NGN(Ne, 1), NGN(Ne, 2))
            With 固端力      '整体座标下的固端力
                fg(1) = .N1 * cosA - .Q1 * sinA
                fg(2) = .N1 * sinA + .Q1 * cosA
                fg(4) = .N2 * cosA - .Q2 * sinA
                fg(5) = .N2 * sinA + .Q2 * cosA
                fg(3) = .M1
                fg(6) = .M2
            End With
            Call VEC(Ne, NV())       '单元Ne的定位向量
            For i = 1 To 6           '对号入座，形成荷载向量
               If NV(i) > 0 Then
                   D(NV(i)) = D(NV(i)) - fg(i)
               End If
           Next i
        End If
    Next Ne
    For Nd = 1 To Nodg              ' 迭加上结点荷载
        For i = 1 To 3
            If XYM(Nd, i) <> "0" And Val(XYM(Nd, i)) <> 0 Then
                n = Neq(Nd, i)
                D(n) = D(n) + XYM(Nd, i)
            End If
        Next i
    Next Nd
End Sub

Sub 形成整体刚度矩阵()
    Dim ES(6, 6) As Single, NV(6) As Integer
    Dim i, j, Ne
    Erase A()
    For Ne = 1 To Ncell
        Call VEC(Ne, NV())         '单元Ne的定位向量NV(1:6)
        Call ESTF(Ne, ES())        '单元Ne的刚度矩阵ES(1:6,1:6)
        For i = 1 To 6
            If NV(i) > 0 Then         'NV(i)-i行对应的方程编号
                For j = 1 To 6
                    If NV(j) > 0 Then 'NV(j)-j列对应的方程编号
                        A(NV(i), NV(j)) = A(NV(i), NV(j)) + ES(i, j)
                    End If            '对号入座
                Next j
            End If
        Next i
    Next Ne
End Sub


Sub Gauss消去法解方程()
    Dim i, j, k, n
    n = Nequ
    For i = 1 To n - 1         '消元次数为N-1次
        For j = i + 1 To n
            For k = i + 1 To n
                 A(j, k) = A(j, k) - A(i, k) * A(j, i) / A(i, i)
            Next k
            D(j) = D(j) - A(j, i) * D(i) / A(i, i)
        Next j
    Next i
    For i = n To 1 Step -1     '回代
        For j = i + 1 To n
            D(i) = D(i) - A(i, j) * D(j)
        Next j
        D(i) = D(i) / A(i, i)
    Next i
End Sub