'===========================================
'Structural Analysis Program For Plan Frame
'===========================================
Option Explicit
Public nn As Integer, ne As Integer
Public nd As Integer, ndf As Integer
Public nf As Integer, npj As Integer
Public npe As Integer, n As Integer
Public al(50) As Double, t(6, 6) As Double
Public x(40) As Double, y(40) As Double
Public jl(50) As Integer, jr(50) As Integer
Public ii(6) As Integer, ea(50) As Double
Public ei(50) As Double, c(6, 6) As Double
Public r(120, 120) As Double, p(120) As Double
Public pe(120) As Double, ibd(20) As Integer
Public bd(20) As Double, mj(20) As Integer
Public qj(20, 3) As Double, ff(6) As Double
Public f(6) As Double, dis(6) As Double
Public mf(30) As Integer, ind(30) As Integer
Public aq(30) As Double, bq(30) As Double
Public q1(30) As Double, q2(30) As Double
Public kec(30) As Integer
'===========================================
'Main Program
'===========================================
Sub frame()
Open "d:\my data\ww.txt" For Input As #1
Open "d:\my data\wt.txt" For Output As #2
Call input1
Call wstiff
Call load
Call bound
Call gauss
Call nqm
Close 1
Close 2
End Sub
'===========================================
'SUB-1 Read and Print Initial Data
'===========================================
Sub input1()
Dim i As Integer, j As Integer
Dim inti As Integer, intj As Integer, ie As Integer
Dim dx As Double, dy As Double
Print #2, "Structural Analysis"
Print #2, "*********************************"
Print #2,
Print #2, "Input Data"
Print #2, "==============="
Print #2,
Print #2, "Structural Control Data"
Print #2, "-------------------------"
Print #2, "nn"; Spc(5); "nf"; Spc(4); "nd"; Spc(3); "ndf"; Spc(4); "ne"; Spc(3); "npj"; Spc(3); "npe"; Spc(4); "n"
Input #1, nn, nf, nd, ndf, ne, npj, npe
n = 3 * (nn - nf)
Print #2, nn; Spc(3); nf; Spc(3); nd; Spc(3); ndf; Spc(3); ne; Spc(3); npj; Spc(3); npe; Spc(3); n
Print #2,
Print #2, "Nodal Coordinates"
Print #2, "------------------------"
Print #2, "Node"; Spc(3); "x"; Spc(5); "y"
For i = 1 To nn
Input #1, i, x(i), y(i)
Print #2, i; Spc(3); x(i); Spc(3); y(i)
Next i
Print #2,
Print #2, "Element Information"
Print #2, "-------------------"
Print #2, "Ele.No"; Spc(3); "jl"; Spc(3); "jr"; Spc(7); "ea"; Spc(8); "ei"; Spc(15); "al"; Spc(13); "kec"
For i = 1 To ne
Input #1, i, jl(i), jr(i), ea(i), ei(i), kec(i)
Next i
For inti = 1 To ne
If jl(inti) >= jr(inti) Then Stop
Next inti
ie = 1
Do While ie <= ne
i = jl(ie)
j = jr(ie)
dx = x(j) - x(i)
dy = y(j) - y(i)
al(ie) = Sqr(dx * dx + dy * dy)
ie = ie + 1
Loop
For i = 1 To ne
Print #2, i; Spc(6); jl(i); Spc(3); jr(i); Spc(3); ea(i); Spc(3); ei(i); Spc(3); al(i); Spc(3); kec(i)
Next i
If npj <> 0 Then
Print #2,
Print #2, "Nodal Load"
Print #2, "--------------------"
Print #2, "i"; Spc(6); "mj"; Spc(4); "xd"; Spc(6); "yd"; Spc(6); "md"
For intj = 1 To npj
Input #1, intj, mj(intj), qj(intj, 1), qj(intj, 2), qj(intj, 3)
Print #2, intj; Spc(3); mj(inti); Spc(3); qj(inti, 1); Spc(5); qj(intj, 2); Spc(5); qj(intj, 3)
Next intj
End If
If npe <> 0 Then
Print #2,
Print #2, "Element loads"
Print #2, "--------------------"
Print #2, Spc(1); "i"; Spc(5); "mf"; Spc(3); "ind"; Spc(4); "aq"; Spc(7); "bq"; Spc(13); "q1"; Spc(10); "q2"
For i = 1 To npe
Input #1, i, mf(i), ind(i), aq(i), bq(i), q1(i), q2(i)
Print #2, i; Spc(3); mf(i); Spc(3); ind(i); Spc(3); aq(i); Spc(5); bq(i); Spc(5); q1(i); Spc(5); q2(i)
Next i
End If
If ndf <> 0 Then
Print #2,
Print #2, "Boundary Conditions"
Print #2, "-------------------------"
Print #2, "i"; Spc(3); "ibd"; Spc(3); "bd"
For j = 1 To ndf
Input #1, j, ibd(j), bd(j)
Print #2, j; Spc(3); ibd(j); Spc(3), bd(j)
Next j
End If
End Sub
'==============================================
'SUB-2 Assemble Structural Stiffness Matrix [R]
'==============================================
Sub wstiff()
Dim i As Integer, j As Integer
Dim ie As Integer, kk As Integer
Dim l As Integer, ct(6, 6) As Double
Dim 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 Until ie > ne
Call stiff(ie)
kk = kec(ie)
If kk <> 1 Then
Call matc(ie, kk, ct)
For i = 1 To 6
For j = 1 To 6
t(i, j) = c(i, j)
Next j
Next i
For i = 1 To 6
For j = 1 To 6
c(i, j) = 0#
Next j
Next i
For i = 1 To 6
For j = 1 To 6
For l = 1 To 6
c(i, j) = c(i, j) + t(i, l) * ct(l, j)
Next l
Next j
Next i
End If
Call locat(ie)
For k1 = 1 To 6
i = ii(k1)
If i > n Then GoTo 30
For k2 = k1 To 6
j = ii(k2)
If j > n Then GoTo 20
r(i, j) = r(i, j) + c(k1, k2)
20 Next k2
30 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 Element Stiffness Matrix in Global Coordinate System [C]
'=====================================================================
Sub stiff(ie)
Dim i As Integer, j As Integer
Dim b1 As Double, b2 As Double
Dim b3 As Double, b4 As Double
Dim s1 As Double, s2 As Double
Dim s3 As Double, s4 As Double
Dim s5 As Double, s6 As Double
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)
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 {II}
'=========================================
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, k As Integer
For i = 1 To n
p(i) = 0#
Next i
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
End If
For i = 1 To n
p(i) = p(i) + pe(i)
Next i
End Sub
'========================
'SUB-6 Set up Effective Nodal Loads
'========================
Sub Eload()
Dim i As Integer, k As Integer
Dim k1 As Integer, kk As Integer
Dim j As Integer, ct(6, 6) As Double
Dim k2 As Integer
Dim l As Integer, ll As Integer
For i = 1 To n
pe(i) = 0#
Next i
For i = 1 To npe
k = mf(i)
Call trans(k)
Call locat(k)
Call efix(i)
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
kk = kec(k)
If kk <> 1 Then
Call matc(k, kk, ct)
For j = 1 To 6
ff(j) = 0#
For l = 1 To 6
ff(j) = ff(j) + ct(l, j) * f(l)
Next l
Next j
For ll = 1 To 6
f(ll) = ff(ll)
Next ll
End If
For k1 = 1 To 6
j = ii(k1)
If j > n Then GoTo 30
pe(j) = pe(j) - f(k1)
Next k1
30 Next i
End Sub
'==============================
'SUB-7 Set up Fixed End Forces
'==============================
Sub efix(i)
Dim j As Integer, k As Integer
Dim sl As Double, a As Double
Dim b As Double, p1 As Double
Dim p2 As Double, b1 As Double
Dim b2 As Double, b3 As Double
Dim c1 As Double, c2 As Double
Dim d1 As Double, d2 As Double
Dim c3 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 ^ 2 - a ^ 2
Select Case ind(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 * 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 * c2 ^ 2 + c2 ^ 2 * sl - 2# * c2 ^ 3 / 15#) / 36# / sl ^ 2
ff(5) = -0.5 * p2 * 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 * (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) = -p1 * d2 / 2# / sl
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 a As Double, i As Integer
Dim k As Integer
If ndf <> 0 Then
a = 1E+20
For i = 1 To ndf
k = ibd(i)
r(k, k) = a
p(k) = a * bd(i)
Next i
End If
End Sub
'=======================
'SUB-10 Solve Equilibrium Equation
'=======================
Sub gauss()
Dim n1 As Integer, k As Integer
Dim k1 As Integer, i As Integer
Dim j As Integer, 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) - p(k) * c
j = k1
Do While j <= n
r(i, j) = r(i, j) - r(k, j) * c
j = j + 1
Loop
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) = p(k) - r(k, j) * p(j)
Next j
p(k) = p(k) / r(k, k)
Next i
Print #2,
Print #2, "Output Data"
Print #2, " "
Print #2, "Nodal Displacements"
Print #2, "---------------------------------"
Print #2, "Node No."; Spc(3); "u"; Spc(34); "v"; Spc(34); "fai"
For i = 1 To nn
Print #2, Spc(2); i; Spc(8); p(3 * i - 2); Spc(8); p(3 * i - 1); Spc(8); p(3 * i)
Next i
End Sub
'========================================
'SUB-11 Calculate the Forces of Elements
'========================================
Sub nqm()
Dim ie As Integer, i As Integer
Dim j As Integer, kk As Integer
Dim k As Integer, ct(6, 6) As Double
Dim l As Integer, jj As Integer
Print #2,
Print #2, "Element.No.and Member end Force:"
Print #2, "Ele.No."; Spc(3); "n(l)"; Spc(3); "q(l)"; Spc(3); "m(l)"; Spc(3); "n(r)"; Spc(3); "q(r)"; Spc(3); "m(r)"
ie = 1
Do While ie <= ne
Call stiff(ie)
Call locat(ie)
For i = 1 To 6
dis i = 0#
j = ii(i)
If j > n Then GoTo 10
dis(i) = p(j)
10 Next i
kk = kec(ie)
If kk <> 1 Then
Call matc(ie, kk, ct)
For i = 1 To 6
For j = 1 To 6
t(i, j) = c(i, j)
Next j
Next i
For i = 1 To 6
For j = 1 To 6
c(i, j) = 0#
For l = 1 To 6
c(i, j) = c(i, j) + t(i, l) * ct(l, j)
Next l
Next j
Next i
End If
For i = 1 To 6
ff(i) = 0#
For j = 1 To 6
ff(i) = ff(i) + c(i, j) * dis(j)
Next j
Next i
Call trans(ie)
For i = 1 To 6
f(i) = 0#
For j = 1 To 6
f(i) = f(i) + t(i, j) * ff(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)
If kk <> 1 Then
For j = 1 To 6
dis(j) = ff(j)
Next j
For j = 1 To 6
ff(j) = 0#
For jj = 1 To 6
ff(j) = ff(j) + ct(jj, j) * dis(jj)
Next jj
Next j
End If
For j = 1 To 6
f(j) = f(j) + ff(j)
Next j
End If
i = i + 1
Loop
End If
Print #2, Spc(2); ie; Spc(5); f(1); Spc(5); f(2); Spc(5); f(3); Spc(5); f(4); Spc(5); f(5); Spc(5); f(6)
ie = ie + 1
Loop
End Sub
'======
'SUB-12
'======
Sub matc(ie, kk, ct)
Dim i As Integer, j As Integer
Dim sl As Double, cl As Double
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
If i = j Then
ct(i, j) = 1#
Else
ct(i, j) = 0#
End If
Next j
Next i
sl = cy / al(ie)
cl = cx / al(ie)
If kk = 2 Then
ct(3, 1) = -3# * sl / 2#
ct(3, 2) = 3# * cl / 2#
ct(3, 3) = 0#
ct(3, 4) = -ct(3, 1)
ct(3, 5) = -ct(3, 2)
ct(3, 6) = -0.5
Else
ct(3, 1) = -sl
ct(3, 2) = cl
ct(3, 3) = 0#
ct(3, 4) = sl
ct(3, 5) = -cl
ct(3, 6) = 0#
ct(6, 1) = -sl
ct(6, 2) = cl
ct(6, 3) = 0#
ct(6, 4) = sl
ct(6, 5) = -cl
ct(6, 6) = 0#
End If
End Sub