module1.bas

较好的工程计算程序功能较强 计算桥梁的跨度等

Attribute VB_Name = "ModuleCompute"
Global CurrentI As Integer '正在定义第CurrentI跨

Global N As Integer '跨数
Global L() As Double '每跨长度
Global EI0 As Double 'EI的基础值
Global EI() As Double  '每跨的EI与dEI的比
Global i0() As Double 'i=EI/l
Global Flag0 As Boolean, flagN As Boolean '结点类型，固结为真，铰接为假

Global Fc() As Integer '每跨上荷载的类型 0为无荷载，1为均布荷载,2为跨中集中荷载
Global M01() As Double, M02() As Double  '每跨左端和右端的固端弯矩
Global q() As Double '跨上荷载的大小
Global PosA() As Double '荷载距离左端相对位置，从0到1

Global K() As Double '刚度矩阵
Global Theta() As Double '转角
Global P() As Double '等效结点荷载

Global m() As Double '结点弯矩，上段纤维受拉为正


'求刚度矩阵K
Public Sub GetK(dK() As Double, di0() As Double, _
          dEI0 As Double, dEI() As Double, dL() As Double)
  Dim i As Integer
  For i = 1 To N
    di0(i) = dEI0 * dEI(i) / dL(i)
  Next
  
    dK(1, 1) = 4 * di0(1)
    dK(1, 2) = 2 * di0(1)
    dK(2, 1) = 2 * di0(1)
  For i = 2 To N
    dK(i, i) = 4 * (di0(i) + di0(i - 1))
    dK(i, i + 1) = 2 * di0(i)
    dK(i + 1, i) = 2 * di0(i)
  Next
    dK(N + 1, N + 1) = 4 * di0(N)
    
End Sub

'求等效结点荷载
Public Sub GetP(dP() As Double, dL() As Double, _
       dFc() As Integer, dQ() As Double, dPosA() As Double)
 
 Dim i As Integer


 
 For i = 1 To N
   Select Case dFc(i)
   Case 0
    M01(i) = 0: M02(i) = 0
   Case 1
    temp = (6 - 8 * dPosA(i) + 3 * dPosA(i) ^ 2)
    M01(i) = -1 * dQ(i) * dPosA(i) ^ 2 * dL(i) ^ 2 * temp / 12
    M02(i) = dQ(i) * dPosA(i) ^ 3 * dL(i) ^ 2 * (4 - 3 * dPosA(i)) / 12
   Case 2
    M01(i) = -1# * dQ(i) * dL(i) * dPosA(i) * (1 - dPosA(i)) ^ 2
    M02(i) = dQ(i) * dL(i) * dPosA(i) ^ 2 * (1 - dPosA(i))
   End Select
 Next
 dP(1) = -1 * M01(1)
 For i = 2 To N
   dP(i) = -1 * (M01(i) + M02(i - 1))
 Next
 dP(N + 1) = -1 * M02(N)
 
End Sub

'解方程 K*theta=P
  
Public Function Formula(dTheta() As Double, dK() As Double, dP() As Double) As Boolean
 ' 函数返回值表示计算是否成功


    Dim i, j As Integer '临时循环变量
    '检查对称性
    For i = 1 To N + 1
     For j = 1 To N + 1
      If K(i, j) <> K(j, i) Then
        respond = MsgBox("刚度矩阵不对称！", vbOKOnly, "错误")
        Formula = False
        Exit Function
      End If
     Next
    Next
    '检查K是否为连续梁的形式
   
      Call Flian(dTheta(), dK(), dP())
      Formula = True
    
    
    
End Function

Public Sub Flian(dTheta() As Double, dK() As Double, dP() As Double)
    '三条对角线定义为三个数组
    Dim tempP() As Double
    ReDim tempP(1 To N + 1)
    
    Dim K1() As Double, K2() As Double, K3() As Double
    ReDim K1(1 To N), K2(1 To N + 1), K3(2 To N + 1)
    Dim a() As Double, b() As Double, d As Integer  'theta(d) = a(i)*theta(d)+b
    ReDim a(1 To N + 1), b(1 To N + 1)
    Dim i As Integer, j As Integer, nd As Integer
    

    For i = 1 To N + 1: tempP(i) = dP(i): Next
    K1(1) = K(1, 2):    K2(1) = K(1, 1)
    For i = 2 To N
      K1(i) = K(i, i + 1)
      K2(i) = K(i, i)
      K3(i) = K(i, i - 1)
    Next
    K2(N + 1) = K(N + 1, N + 1):    K3(N + 1) = K(N + 1, N)
    
    '处理支座条件
    If Flag0 Then
      dK(1, 1) = 1: dK(1, 2) = 0: dK(2, 1) = 0: tempP(1) = 0
      dTheta(1) = 0: d = 2
      
    Else
      d = 1
    End If
    If flagN Then
      dK(N + 1, N + 1) = 1: dK(N + 1, N) = 0: dK(N + 1, N) = 0: tempP(N + 1) = 0
      dTheta(N + 1) = 0: nd = N
    Else
      nd = N + 1
    End If
    '解方程
    If N > 1 Then
        a(d) = 1: b(d) = 0
        a(d + 1) = -1 * K2(d) / K1(d): b(d + 1) = tempP(d) / K1(d)
        For i = d + 1 To nd - 1
         a(i + 1) = -1 * (K3(i) * a(i - 1) + K2(i) * a(i)) / K1(i)
         b(i + 1) = (tempP(i) - K3(i) * b(i - 1) - K2(i) * b(i)) / K1(i)
        Next
        dTheta(d) = (tempP(nd) - K3(nd) * b(nd - 1) - K2(nd) * b(nd)) / (K3(nd) * a(nd - 1) + K2(nd) * a(nd))
    
        For i = d + 1 To nd
            dTheta(i) = a(i) * dTheta(d) + b(i)
        Next
    Else
        If Flag0 = False And flagN = False Then
           a(2) = -1 * K(1, 1) / K(1, 2)
           b(2) = tempP(1) / K(1, 2)
           dTheta(1) = (tempP(2) - b(2) * K(2, 2)) / (K(2, 1) + K(2, 2) * a(2))
           dTheta(2) = a(2) * dTheta(1) + b(2)
        Else
           dTheta(1) = tempP(1) / K(1, 1)
           dTheta(2) = tempP(2) / K(2, 2)
        End If
    End If
  End Sub
  
  '计算杆端弯矩
  Public Sub GetM(dM() As Double, di0() As Double, dTheta() As Double, dP() As Double)
  
  Dim i As Integer
  
  If Flag0 Then
      dM(1) = dP(1) + (-2) * di0(1) * dTheta(2)
  Else
      dM(1) = 0
  End If
  
  For i = 2 To N
    dM(i) = -1 * M01(i) - (4 * di0(i) * dTheta(i) + 2 * di0(i) * dTheta(i + 1))
  Next
  
  If flagN Then
     dM(N + 1) = -1 * dP(N + 1) + 2 * di0(N) * dTheta(N)
  End If

  End Sub