结构矩阵分析代码

'main program

'====================================== Sub frame()

Open \"e:\\PAD\\fr.txt\" For Input As 1 Open \"e:\\PAD\\fw.txt\" For Output As 2

Call intput1 Call wstiff Call load Call bound Call gauss Call nqm

Close 1 Close 2

Label1.Caption = \"adadsfdsff\" End Sub

'================================ 'sub-1read and print initial data

'================================

'111111111111111111111111111111111111111111111111111111111111 Sub input1()

Dim inti As Integer, intj As Integer, i As Integer, j As Integer, k As Integer Dim dx, dy As Double

Print #2, \"plane frame structural analysis\" Print #2, \"*************************\" Print #2, \"input data\" Print #2, \"=====\"

Print #2, \"structural control data\" Print #2, \"---------------\"

Print #2, \"nn\"; Spc(3); \"ne\"; Spc(3); \"nf\"; Spc(3); \"nd\"; Spc(3); \"ndf\"; Spc(2); \"npj\"; Spc(2); \"npe\"; Spc(3); \"n\"

Input #1, nn, ne, nf, nd, ndf, npj, npe n = 3 * (nn - nf)

Print #2, \"nn\"; Spc(2); \"ne\"; Spc(2); \"nf\"; Spc(2); \"nd\"; Spc(2); \"ndf\"; Spc(2); \"npj\"; Spc(2); \"npe\"; Spc(2); \"n\"

Print #2,

Print #2, \"nodal coordinates\" Print #2, \"---------------\"

Print #2, \"node\"; Spc(2); \"x\"; Spc(5); \"y\" i = nn

For inti = 1 To i

Input #1, inti, x(inti), y(inti)

Print #2, inti, Spc(2); x(inti); Spc(3); y(inti) Next inti

Print #2,

Print #2, \"element information\" Print #2, \"-------------------\"

Print #2, \"ele.no.\"; Spc(2); \"jl\"; Spc(2); \"jr\"; Spc(5); \"ea\"; Spc(7); \"ei\"; Spc(5); \"al\" i = ne

For inti = 1 To i

Input #1, inti, jl(inti), jr(inti), ea(inti), ei(inti) Next inti

For inti = 1 To i

If jl(inti) >= jr(inti) Then Stop Next inti

For inti = 1 To i j = jl(inti) k = jr(inti) dx = x(k) - x(j) dy = y(k) - y(j)

al(inti) = Sqr(dx * dx + dy * dy)

Print #2, Spc(3); inti; Spc(4); jl(inti); Spc(3); jr(inti); Spc(2); ea(inti); Spc(2); ei(inti); Spc(2); al(inti) Next inti Print #2, k = npj

If k <> 0 Then

Print #2, \"Nodal Load\" Print #2, \"-----------------\"

Print #2, \"i\"; Spc(3); \"mj\"; Spc(2); \"xd\"; Spc(2); \"yd\"; Spc(2); \"md\" For inti = 1 To k

Input #1, inti, mj(inti), qj(inti, 1), qj(inti, 2), qj(inti, 3)

Print #2, inti; Spc(2); mj(inti); Spc(2); qj(inti, 1); Spc(2); qj(inti, 2); Spc(2); qj(inti, 3) Next inti End If Print #2, i = npe

If i <> 0 Then

Print #2, \"Element load\"

Print #2, \"----------------------\"

Print #2, \"i\"; Spc(3); \"mf\"; Spc(2); \"ind\"; Spc(2); \"aq\"; Spc(2); \"bq\"; Spc(2); \"q1\"; Spc(3); \"q2\"

For inti = 1 To i

Input #1, inti, mf(inti), ind(inti), aq(inti), bq(inti), q1(inti), q2(inti)

Print #2, inti; Spc(2); mf(inti); Spc(3); ind(inti); Spc(2); aq(inti); Spc(2); bq(inti); Spc(2); q1(inti); Spc(3); q2(inti) Next inti End If Print #2, j = ndf

If j <> 0 Then

Print #2, \"Boundary Condition\" Print #2, \"-----------------------\"

Print #2, \"i\"; Spc(3); \"ibd\"; Spc(3); \"bd\" For inti = 1 To j

Input #1, inti, ibd(inti), bd(inti)

Print #2, inti; Spc(3); ibd(inti); Spc(3); bd(inti) Next inti End If

End Sub

'============================================== 'SUB-2 Assemble Structural Stiffness Matrix{R}装配结构刚度

'============================================== Sub wstiff()

Dim i As Integer, j As Integer, ie As Integer, k1 As Integer, k2 As Integer

For i = 1 To n For j = 1 To n r(i, j) = 0 Next j Next i ie = 1

Do While ie <= ne Call stiff(ie) Call locat(ie) For k1 = 1 To 6 i = ii(k1)

If i <= n Then

For k2 = k1 To 6

j = ii(k2)

If j <= n Then

r(i, j) = r(i, j) + c(k1, k2) End If Next k2 End If Next k1 ie = ie + 1 Loop

For i = 2 To n

For j = 1 To (i - 1) r(i, j) = r(j, i) Next j Next i End Sub

'================================================= 'SUB-3 Set Up Stiffness Matrix[C] (Global Coordinate System) '建立了刚度矩阵

'=================================================

Sub stiff(ie)

Dim i As Integer, j As Integer

Dim cx As Double, cy As Double, b1 As Double, b2 As Double, b3 As Double, b4 As Double Dim s1 As Double, s2 As Double, s3 As Double, s4 As Double, s5 As Double, s6 As Double

i = jl(ie) j = jr(ie)

cx = (x(j) - x(i)) / al(ie) cy = (y(j) - y(i)) / al(ie) b1 = ea(ie) / al(ie)

b2 = 12# * ei(ie) / al(ie) ^ 3 b3 = 6# * ei(ie) / al(ie) ^ 2 b4 = 2# * ei(ie) / al(ie)

s1 = b1 * cx ^ 2 + b2 * cy ^ 2 s2 = (b1 - b2) * cx * cy s3 = b3 * cy

s4 = b1 * cy ^ 2 + b2 * cx ^ 2 s5 = b3 * cx s6 = b4 c(1, 1) = s1

c(1, 2) = s2 c(1, 3) = s3

'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ c(1, 4) = -s1 c(1, 5) = -s2 c(1, 6) = s3 c(2, 2) = s4 c(2, 3) = -s5 c(2, 4) = -s2 c(2, 5) = -s4 c(2, 6) = -s5 c(3, 3) = 2# * s6 c(3, 4) = -s3 c(3, 5) = s5 c(3, 6) = s6 c(4, 4) = s1 c(4, 5) = s2 c(4, 6) = -s3 c(5, 5) = s4 c(5, 6) = s5 c(6, 6) = 2# * s6 For i = 2 To 6

For j = 1 To (i - 1) c(i, j) = c(j, i) Next j Next i End Sub

'SUB-4 set up element location vector{11}设置元素位置矢量

Sub locat(ie)

Dim i As Integer, j As Integer i = jl(ie) j = jr(ie)

ii(1) = 3 * i - 2 ii(2) = 3 * i - 1 ii(3) = 3 * i ii(4) = 3 * j - 2 ii(5) = 3 * j - 1 ii(6) = 3 * j End Sub

'SUB-5 set up total nodal vector {p}

Sub load()

Dim i As Integer, j As Integer, k As Integer End Sub i = 1

Do While i <= n p(i) = 0# i = i + 1 Loop

If npj > 0 Then

For i = 1 To npj k = mj(i)

p(3 * k - 2) = qj(i, 1) p(3 * k - 1) = qj(i, 2) p(3 * k) = qj(i, 3) Next i End If

If npe <> 0 Then Call eload i = 1

Do While i <= n

p(i) = p(i) + p(ei) i = i + 1 Loop End If

End Sub

'=====================

'USB-6 set up element effective load设置元素的有效载荷

'===================== Sub eload()

Dim i As Integer, j As Integer, k As Integer, k1 As Integer, k2 As Integer, k3 As Integer i = 1

Do While i <= n pe(i) = 0# i = i + 1 Loop j = 1

Do While j <= npe k = mf(j) Call trans(k) Call locat(k) Call efix(j)

For k1 = 1 To 6 f(k1) = 0#

For k2 = 1 To 6

f(k1) = f(k1) + t(k2, k1) * ff(k2) Next k2 Next k1

For k3 = 1 To 6 i = ii(k3)

If i <= n Then

pe(i) = pe(i) - f(k3) End If Next k3 j = j + 1

'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Loop

End Sub

'=====================

'SUB-7 set up fixed-end force of element '=====================

Sub efix(i)

Dim j As Integer, k As Integer

Dim s1 As Double, a As Double, p1 As Double, p2 As Double

Dim b1 As Double, b2 As Double, b3 As Double, c1 As Double, c2 As Double, c3 As Double Dim d1 As Double, d2 As Double

For j = 1 To 6 ff(j) = 0# Next j k = mf(i) s1 = a1(k) a = aq(i) b = bq(i) p1 = q1(i) p2 = q2(i)

b1 = s1 - (a + b) / 2#

b2 = b - a

b3 = (a + b) / 2#

c1 = s1 - (2# * b + a) / 3# c2 = b2

c3 = (2# * b + a) / 3# d1 = b ^ 3 - a ^ 3 d2 = b * b - a * a Select Case ind(i) Case 1

ff(2) = -p1 * (s1 - a) ^ 2 * (1# + 2# * a / s1) / s1 ^ 2 ff(3) = p1 * a * (s1 - a) ^ 2 / s1 ^ 2 ff(5) = -p1 - ff(2)

ff(6) = -p1 * a ^ 2 * (s1 - a) / s1 ^ 2 Case 2

ff(2) = -p1 * b2 * (12# * b1 ^ 2 * s1 - 8# * b1 ^ 3 + b2 ^ 2 * s1 - 2# * b1 * b2 ^ 2) / (4# * s1 ^ 3)

ff(3) = p1 * b2 * (12# * b3 * b1 ^ 2 - 3 * b1 * b2 ^ 2 + b2 ^ 2 * s1) / 12# / s1 ^ 2 ff(5) = -p1 * b2 - ff(2)

ff(6) = -p1 * b2 * (12# * b3 ^ 2 * b1 + 3# * b1 * b2 ^ 2 - 2# * b2 ^ 2 * s1) / 12# / s1 ^ 2 Case 3

ff(2) = -p2 * c2 * (18 * c1 ^ 2 * s1 - 12 * c1 ^ 3 + c2 ^ 2 * s1 + 2 * c1 * c2 ^ 2 - 4 * c2 ^ 3 / 45) / 12 / s1 ^ 3

ff(3) = p2 * c2 * (18# * c3 * c1 ^ 2 - 3# * c1 * c2 ^ 2 + c2 ^ 2 * s1 - 2# * c2 ^ 3 / 15#) / 36# / s1 ^ 2

ff(5) = -0.5 * p1 * c2 - ff(2)

ff(6) = -p2 * c2 * (18# * c3 ^ 2 * c1 + 3 * c1 * c2 ^ 2 - 2 * c2 ^ 2 * s1 + 2 * c2 ^ 3 / 15#) / 36# / s1 ^ 2 Case 4

ff(2) = -6# * p1 * a * (s1 - a) / s1 ^ 3

ff(3) = p1 * (s1 - a) * (3# * a - s1) / s1 ^ 2 ff(5) = -ff(2)

ff(6) = p1 * a * (2# * s1 - 3# * a) / s1 ^ 2 Case 5

ff(2) = -p1 * (3# * s1 * d2 - 2# * d1) / s1 ^ 3 ff(3) = p1 * (2# * d2 + (b - a) * s1 - d1 / s1) / s1 ff(5) = -ff(2)

ff(6) = p1 * (d2 - d1 / s1) / s1 Case 6

ff(1) = -p1 * (1# - a / s1) ff(4) = -p1 * a / s1 Case 7

ff(1) = -p1 * (b - a) * (1# - (b - a) / (2# * s1)) ff(4) = -p1 * d2 / 2# / s1 Case 8

ff(3) = -a * (p1 - p2) * ei(k) / b ff(6) = -ff(3) End Select End Sub

'====================

'sub-8 set up coordinate transfer matrix[t] '==================== Sub trans(ie)

Dim i As Integer, j As Integer Dim cx As Double, cy As Double i = jl(ie) j = jr(ie)

cx = (x(j) - x(i)) / al(ie) cy = (y(j) - y(i)) / al(ie) For i = 1 To 6 For j = 1 To 6 t(i, j) = 0# Next j Next i

For i = 1 To 4 Step 3 t(i, i) = cx t(i, i + 1) = cy t(i + 1, i) = -cy t(i + 1, i + 1) = cx t(i + 2, i + 2) = 1# Next i

End Sub

'====================

'sub-9 introduce support conditions '==================== Sub bound()

Dim i As Integer, j As Integer, k As Integer Dim a As Double

If ndf <> 0 Then For j = 1 To ndf a = 1E+20

For i = 1 To ndf k = ibd(i) r(k, k) = a

p(k) = a * bd(i) Next i Next j

End If

End Sub

'= = = = = = = = = = = = = = = = = 'SUB-10 Solve Equilibrium Equations '= = = = = = = = = = = = = = = = = Sub gauss()

Dim i As Integer, j As Integer, k As Integer, k1 As Integer, n1 As Integer Dim c As Double

n1 = n - 1

For k = 1 To n1 k1 = k + 1 For i = k1 To n

c = r(k, i) / r(k, k) p(i) = p(i) - (k) * c For j = k1 To n

r(i, j) = r(i, j) - r(k, j) * c Next j Next i Next k

p(n) = p(n) / r(n, n) For i = 1 To n1 k = n - i k1 = k + 1

For j = k1 To n

p (k) - r(k, j) * p(j) Next j

p(k) = p(k) / r(k, k) Next i

For k = 1 To ne Call locat(k) For k1 = 1 To 6 i = ii(k1) If i > n Then p(i) = 0# End If Next k1 Next k Print #2,

Print #2, \"Output Data\"

Print #2, \"===============\" Print #2,

Print #2, \"Nodal Displacement\" Print #2, \"----------------\"

Print #2, \"Node No.\For i = 1 To nn

Print #2, i, p(3 * i - 2), p(3 * i - 1), p(3 * i) Next

End Sub

'= = = = = = = = = = = = = = = = = = = = = = = =

'SSUB-11 Calculate Member-End Forces Of Elements '= = = = = = = = = = = = = = = = = = = = = = = = Sub nqm()

Dim i As Integer, j As Integer, ie As Integer, k As Integer Dim fe(6) As Double, f1(6) As Double

Print #2, \"Element No.&member-end force\" Print #2, \"= = = = = = = = = = = =\" Print #2,

Print #2, \"Ele No.\Print #2, \"-----------------------------------\" ie = 1

Do While ie <= ne Call trans(ie) Call stiff(ie) Call locat(ie) For i = 1 To 6 dis(i) = 0# j = ii(i) If j <= n Then dis(i) = p(j) End If Next i

For i = 1 To 6 fe(i) = 0# For j = 1 To 6

fe(i) = fe(i) + c(i, j) * dis(j) Next j Next i

For i = 1 To 6 f(i) = 0#

For j = 1 To 6

f(i) = f(i) + t(i, j) * fe(j)

Next j Next i

If npe < 0# Then i = 1

Do While i <= npe k = mf(i) If k = ie Then Call efix(i) For j = 1 To 6 f(j) = f(j) + ff(j) Next j End If i = i + 1 Loop End If

Print #2, ie, f(1), f(2), f(3), f(4), f(5), f(6) ie = ie + 1 Loop

Label1.Caption = \"sadfadsfd\" End Sub