结构矩阵原理分析与程序设计程序代码.docx

1、结构矩阵原理分析与程序设计程序代码Option ExplicitDim nn As Integer, ne As Integer, nd As Integer, ndf As IntegerDim nf As Integer, npj As Integer, npe As Integer, n As IntegerDim al(50) As Double, t(6, 6) As Double, x(40) As Double, y(40) As DoubleDim jl(50) As Integer, jr(50) As Integer, ea(50) As Double, ei(50) As

2、 DoubleDim c(6, 6) As Double, r(120, 120) As Double, p(120) As Double, pe(120) As DoubleDim ibd(20) As Integer, ii(6) As Integer, bd(20) As Double, ff(6) As DoubleDim mj(20) As Integer, qj(20, 3) As Double, f(6) As Double, dis(6) As DoubleDim mf(30) As Integer, ind(30) As Integer, aq(30) As Double,

3、bq(30) As DoubleDim q1(30) As Double, q2(30) As DoubleSub frame()Open d:my dataww.txt For Input As #1Open d:my datawt.txt For Output As #2Call input1Call wstiffCall loadCall boundCall gaussCall nqmClose 1Close 2End SubSub input1()Dim inti As Integer, intj As Integer, i As Integer, j As Integer, k As

4、 IntegerDim dx, dy As DoublePrint #2, Plane Frame Structural AnalysisPrint #2, *Print #2, Input DataPrint #2, =Print #2,Print #2, Structural Control DataPrint #2, -Print #2, nn; Spc(3); ne; Spc(3); nf; Spc(3); nd; Spc(3); ndf; Spc(2); npj; Spc(2); npe; Spc(3); nInput #1, nn, ne, nf, nd, ndf, npj, np

5、en = 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); nPrint #2,Print #2, Nodal CoordinatesPrint #2, -Print #2, Node; Spc(2); x; Spc(5); yi = nnFor inti = 1 To i Input #1, inti, x(inti), y(inti) Print #2, inti; Spc(2); x(inti); Spc(3); y(i

6、nti)Next intiPrint #2,Print #2, Element InformationPrint #2, -Print #2, Ele.No; Spc(2); jl; Spc(2); jr; Spc(5); ea; Spc(7); ei; Spc(5); ali = neFor inti = 1 To i Input #1, inti, jl(inti), jr(inti), ea(inti), ei(inti)Next intiFor inti = 1 To i If jl(inti) = jr(inti) Then StopNext intiFor inti = 1 To

7、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 intiPrint #2,k = npjIf k 0 Then Print #2, Nodal Load Print #2, - Print #2, i; Spc(3);

8、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 intiEnd IfPrint #2,i = npeIf i 0 Then Print #2, Element loads Print #2, - Pri

9、nt #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)

10、Next intiEnd If Print #2,j = ndfIf j 0 Then Print #2, Boundary Conditions 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 intiEnd IfEnd SubSub wstiff()Dim i As Integer, j As Integer, ie As

11、 Integer, k1 As Integer, k2 As IntegerFor i = 1 To n For j = 1 To n r(i, j) = 0 Next j Next i ie = 1Do 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 +

12、1LoopFor i = 2 To n For j = 1 To (i - 1) r(i, j) = r(j, i) Next jNext iEnd SubSub stiff(ie)Dim i As Integer, j As IntegerDim cx As Double, cy As Double, b1 As Double, b2 As Double, b3 As Double, b4 As DoubleDim s1 As Double, s2 As Double, s3 As Double, s4 As Double, s5 As Double, s6 As Doublei = jl(

13、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) 3b3 = 6# * ei(ie) / al(ie) 2b4 = 2# * ei(ie) / al(ie)s1 = b1 * cx 2 + b2 * cy 2s2 = (b1 - b2) * cx * cys3 = b3 * cys4 = b1 * cy 2 + b2 * cx 2s5 = b3 * cxs6 = b4c(1, 1) = s1c(1, 2) = s2c(1,

14、 3) = s3c(1, 4) = -s1c(1, 5) = -s2c(1, 6) = s3c(2, 2) = s4c(2, 3) = -s5c(2, 4) = -s2c(2, 5) = -s4c(2, 6) = -s5c(3, 3) = 2# * s6c(3, 4) = -s3c(3, 5) = s5c(3, 6) = s6c(4, 4) = s1c(4, 5) = s2c(4, 6) = -s3c(5, 5) = s4c(5, 6) = s5c(6, 6) = 2# * s6For i = 2 To 6 For j = 1 To (i - 1) c(i, j) = c(j, i) Next

15、 jNext iEnd Sub=SUB-4 set up element location vector (II)=Sub locat(ie)Dim i As Integer, j As Integer i = jl(ie)j = jr(ie)ii(1) = 3 * i - 2ii(2) = 3 * i - 1ii(3) = 3 * iii(4) = 3 * j - 2ii(5) = 3 * j - 1ii(6) = 3 * jEnd Sub=SUB-5 set up total nodal vevtor p=Sub load()Dim i As Integer, j As Integer,

16、k As Integeri = 1Do While i 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 iEnd IfIf npe 0 Then Call eload i = 1 Do While i = n p(i) = p(i) + pe(i) i = i + 1 LoopEnd IfEnd Sub=SUB-6 set up element effective load=Sub eload()Dim i As Integer,

17、 j As Integer, k As Integer, k1 As Integer, k2 As Integer, k3 As Integeri = 1Do While i = n pe(i) = 0# i = i + 1Loopj = 1Do 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 k2Next k1For k3 = 1 To 6 i =

18、 ii(k3)If i = n Then pe(i) = pe(i) - f(k3) End IfNext k3j = j + 1LoopEnd Sub=sub-7 set up fixed-end force of element=Sub efix(i)Dim j As Integer, k As IntegerDim sl As Double, a As Double, b As Double, p1 As Double, p2 As DoubleDim b1 As Double, b2 As Double, b3 As Double, c1 As Double, c2 As Double

19、, c3 As DoubleDim d1 As Double, d2 As Double For j = 1 To 6 ff(j) = 0# Next j k = mf(i) sl = al(k) a = aq(i) b = bq(i) p1 = q1(i) p2 = q2(i) b1 = sl - (a + b) / 2# b2 = b - a b3 = (a + b) / 2# c1 = sl - (2# * b + a) / 3# c2 = b2 c3 = (2# * b + a) / 3# d1 = b 3 - a 3 d2 = b * b - a * a Select Case in

20、d(i) Case 1 ff(2) = -p1 * (sl - a) 2 * (1# + 2# * a / sl) / sl 2 ff(3) = p1 * a * (sl - a) 2 / sl 2 ff(5) = -p1 - ff(2) ff(6) = -p1 * a 2 * (sl - a) / sl 2 Case 2 ff(2) = -p1 * b2 * (12# * b1 2 * sl - 8# * b1 3 + b2 2 * sl - 2# * b1 * b2 2) / (4# * sl 3) ff(3) = p1 * b2 * (12# * b3 * b1 2 - 3 * b1 *

21、 b2 2 + b2 2 * sl) / 12# / sl 2 ff(5) = -p1 * b2 - ff(2) ff(6) = -p1 * b2 * (12# * b3 2 * b1 + 3# * b1 * b2 2 - 2# * b2 2 * sl) / 12# / sl 2 Case 3 ff(2) = -p2 * c2 * (18 * c1 2 * sl - 12 * c1 3 + c2 2 * sl - 2 * c1 * c2 2 - 4 * c2 3 / 45) / 12 / sl 3 ff(3) = p2 * c2 * (18# * c3 * c1 2 - 3# * c1 * c

22、2 2 + c2 2 * sl - 2# * c2 3 / 15#) / 36# / sl 2 ff(5) = -0.5 * p1 * c2 - ff(2) ff(6) = -p2 * c2 * (18# * c3 2 * c1 + 3 * c1 * c2 2 - 2 * c2 2 * sl + 2 * c2 3 / 15#) / 36# / sl 2 Case 4 ff(2) = -6# * p1 * a * (sl - a) / sl 3 ff(3) = p1 * (sl - a) * (3# * a - sl) / sl 2 ff(5) = -ff(2) ff(6) = p1 * a *

23、 (2# * sl - 3# * a) / sl 2 Case 5 ff(2) = -p1 * (3# * sl * d2 - 2# * d1) / sl 3 ff(3) = p1 * (2# * d2 + (b - a) * sl - d1 / sl) / sl ff(5) = -ff(2) ff(6) = p1 * (d2 - d1 / sl) / sl Case 6 ff(1) = -p1 * (1# - a / sl) ff(4) = -p1 * a / sl Case 7 ff(1) = -p1 * (b - a) * (1# - (b + a) / (2# * sl) ff(4)

24、= -p1 * d2 / 2# / sl Case 8 ff(3) = -a * (p1 - p2) * ei(k) / b ff(6) = -ff(3) End SelectEnd Sub = sub-8 = 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#

25、 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=Sub bound()Dim i As Integer, j As Integer, k As IntegerDim a As DoubleIf ndf 0 Then For j = 1 To ndf a = 1E+20 For i = 1 To ndf k = ibd(i) r(k, k) = a p

26、(k) = a * bd(i) Next i Next j End IfEnd Sub=SUB-10 Solve Equilibrium Equations=Sub gauss()Dim i As Integer, j As Integer, k As Integer, k1 As Integer, n1 As IntegerDim c As Doublen1 = n - 1For k = 1 To n1 k1 = k + 1 For i = k1 To n c = r(k, i) / r(k, k) p(i) = p(i) - p(k) * c For j = k1 To n r(i, j)

27、 = r(i, j) - r(k, j) * c Next j Next iNext kp(n) = p(n) / r(n, n)For i = 1 To n1 k = n - i k1 = k + 1 For j = k1 To n p(k) = p(k) - r(k, j) * p(j) Next j p(k) = p(k) / r(k, k) Next iFor 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 k1Next kPrint #2,Print #2, Output DataPrint #2, =Print #2,Print #2, Nodal DisplacementPrint #2, -Print #2, Node No., u, v, faiFor i = 1 To nn Print #2, i, p(3 * i - 2), p(3 * i - 1), p(3 * i)NextEnd Sub=SUB-11 Calculate