有限元源程序

Public E As Integer, tt As Integer, I As Integer

Private Sub Command1_Click()

' 平面刚架结构计算程序（FEM） '

FileName$ = InputBox("

输入数据原始文件名:") FileName$ = "F:\VB2\

工程计算\DATA1.TXT"

Open FileName$ For Input As #1

' 读入节点和单元信息

Input #1, NJ, NE, NP, NF, NI

For E = 1 To NE: Input #1, JM(E, 1), JM(E, 2): Next E

For E = 1 To NE: Input #1, DH(E, 1), DH(E, 2), DH(E, 3), DH(E, 4), DH(E, 5), DH(E,

6): Next E

For E = 1 To NE: Input #1, GP(E, 1), GP(E, 2), GP(E, 3), GP(E, 4), GP(E, 5): Next E

If NP <>0 Then

For I = 1 To NP: Input #1, JP(I), PJ(I): Next I

End If

If NF <>0 Then

For I = 1 To NF: Input #1, PF(I, 1), PF(I, 2), PF(I, 3), PF(I, 4), PF(I, 5): Next I

End If

Close #1

‘返回节点和单元信息，并形成节点力列阵'

FileName$ = InputBox("

输入数据结果文件名

:") FileName$ = "F:\VB2\

工程计算

\DATA2.TXT"

Open FileName$ For Output As #2

Print #2, "节点和单元信息"

Print #2, NJ, NE, NP, NF, NI

For E = 1 To NE: Print #2, JM(E, 1), JM(E, 2): Next E

For E = 1 To NE: Print #2, DH(E, 1), DH(E, 2), DH(E, 3), DH(E, 4), DH(E, 5), DH(E,

6): Next E

For E = 1 To NE: Print #2, GP(E, 1), GP(E, 2), GP(E, 3), GP(E, 4), GP(E, 5): Next E

For I = 1 To NI: P(I) = 0: Next I

If NP <>0 Then

For I = 1 To NP: Print #2, JP(I), PJ(I): Next I

For I = 1 To NP: II = JP(I): P(II) = P(II) + PJ(I): Next I

End If

If NF <>0 Then

For I = 1 To NF: Print #2, PF(I, 1), PF(I, 2), PF(I, 3), PF(I, 4), PF(I, 5): Next I

For I = 1 To NF: E = PF(I, 1): tt = PF(I, 2)

Call zbzh(E, t())

If tt = 1 Then Call zyhz1(I, FF())

If tt = 2 Then Call zyhz2(I, FF())

If tt = 3 Then Call zyhz3(I, FF())

If tt = 4 Then Call zyhz4(I, FF())

For L1 = 1 To 6

PE(L1) = 0

For M1 = 1 To 6

PE(L1) = PE(L1) - t(M1, L1) * FF(M1)

Next M1

Next L1

For j = 1 To 6

II = DH(E, j)

If II >0 Then P(II) = P(II) + PE(j)

Next j

Next I

End If

' 调用单元刚度矩阵子程序，组集结构刚度矩阵

For I = 1 To NI

For j = 1 To NI

K(I, j) = 0

Next j

Next I

For E = 1 To NE

Call dygd(E, ke())

For H = 1 To 6

jh = DH(E, H)

If jh <>0 Then

For L = 1 To 6

JL = DH(E, L)

If JL <>0 Then

K(jh, JL) = K(jh, JL) + ke(H, L)

End If

Next L

End If

Next H

Next E

' 调用高斯消元法求解总刚方程，得到节点位移列阵

n = NI

For K1 = 1 To n - 1: For I = K1 + 1 To n: C = -K(I, K1) / K(K1, K1): P(I) = P(I) + P(K1)

* C

For j = K1 + 1 To n: K(I, j) = K(I, j) + K(K1, j) * C: Next j: Next I

Next K1

D(n) = P(n) / K(n, n)

For I = n - 1 To 1 Step -1

For j = I + 1 To n: P(I) = P(I) - K(I, j) * D(j): Next j

D(I) = P(I) / K(I, I)

Next I

' 输出节点位移

Print "节点位移计算结果："

Print #2, "节点位移计算结果："

For K1 = 1 To NI

Print K1, D(K1)

Print #2, "第"K1"个节点位移="D(K1)

Next K1

' 计算各单元杆端内力（轴力、剪力和弯矩）

Print "单元杆端内力（轴力、剪力和弯矩）计算结果："

Print #2, "单元杆端内力（轴力、剪力和弯矩）计算结果："

For E = 1 To NE: Call dygd(E, ke()): Call zbzh(E, t())

For II = 1 To 6: JJ = DH(E, II)

If JJ = 0 Then

DE(II) = 0

Else

DE(II) = D(JJ)

End If

Next II

For I = 1 To 6

FE(I) = 0

For j = 1 To 6

For K1 = 1 To 6

FE(I) = FE(I) + t(I, j) * ke(j, K1) * DE(K1)

Next K1: Next j: Next I

If NF <>0 Then

For I = 1 To NF

If PF(I, 1) = E Then

tt = PF(I, 2)

If tt = 1 Then Call zyhz1(I, FF())

If tt = 2 Then Call zyhz2(I, FF())

If tt = 3 Then Call zyhz3(I, FF())

If tt = 4 Then Call zyhz4(I, FF())

For j = 1 To 6: FE(j) = FE(j) + FF(j): Next j

End If

Next I

End If

If GP(E, 2) = 0 Then

Print #2, "单元："E, "节点："JM(E, 1)"--->"JM(E, 2), "轴力单元"

Print #2, "N1="-FE(1)"N2="FE(4)

Print "单元："E, "节点："JM(E, 1)"--->"JM(E, 2), "轴力单元"

Print "N1="-FE(1)"N2="FE(4)

End If

If GP(E, 2) <>0 Then

If GP(E, 3) = 0 Then

Print #2, "单元："E, "节点："JM(E, 1)"--->"JM(E, 2), "弯曲单元"

轻型门式刚架计算软件是一种用于进行轻型门式刚架结构设计和计算的计算机程序，它能够自动执行许多不同类型的结构分析，包括静力学、动力学和非线性分析等。此外，它还可以自动绘制出相应的结构荷载、节点位移、剪力、弯矩等分布图。

轻型门式刚架结构是一种广泛应用于物流仓库、工业厂房、大型商场等建筑结构中的组合式框架结构，它由一些水平和垂直的杆件、连接件和支撑件组成。与传统的钢结构相比，轻型门式刚架结构具有自重轻、强度高、施工快等优点。因此，在现代建筑设计和施工中得到了广泛应用。

轻型门式刚架计算软件可以帮助设计人员快速、准确地计算和设计轻型门式刚架结构，并生成相应的结构图纸和计算报告。它可以提高设计效率，减少设计成本，同时确保设计质量和结构安全性。

你提出的问题我之前刚好做过，使用有限元方法来进行桁架结构分析。

Matlab编程实现平面杆单元分析

首先，明确Matlab程序要实现的5个重要模块分别为：单元刚度矩阵的求解、单元组装、节点位移的求解、单元应力的求解、节点力的求解。下面给出这5个模块的实现。

1． 单元刚度矩阵求解

定义函数Bar2D2Node_Stiffness，该函数计算单元的刚度矩阵，输入d性模量E，横截面积A，两个节点坐标输出单元刚度矩阵k(4X4)。具体代码如下：

function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2)

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))

x=acos((x2-x1)/L)

C=cos(x)

S=sin(x)

k=E*A/L*[C*C C*S -C*C -C*S C*S S*S -C*S -S*S

-C*C -C*S C*C C*S -C*S -S*S C*S S*S]

2． 单元组装

定义函数Bar2D2Node_Assembly，该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j。输出整体刚度矩阵KK，具体代码如下：

function z = Bar2D2Node_Assembly(KK,k,i,j)

DOF(1)=2*i-1

DOF(2)=2*i

DOF(3)=2*j-1

DOF(4)=2*j

for n1=1:4

for n2=1:4

KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2)

end

end

z=KK

3． 节点位移的求解

定义函数Bar2D2Node_Disp(KK,num,p)，该函数输入KK为总体刚度矩阵；num为活动自由度编号数组；p为活动自由度方向上的节点力；输出节点位移列阵。具体代码如下：

function u = Bar2D2Node_Disp(KK,num,p)

k=KK(num,num)

u=k\p

4． 单元应力的求解

定义函数函数Bar2D2Node_Stress(E,x1,y1,x2,y2,u)，该函数计算单元的应力输入d性模量E，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2）单元节点位移矢量u，返回单元应力标量stress 。具体代码如下：

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))

x=acos((x2-x1)/L)

C=cos(x)

S=sin(x)

stress=E/L*[-C -S C S]*u

5． 计算节点力

定义函数Bar2D2Node_Forces(KK,q)，该函数用于计算节点力，KK为刚度矩阵，q为节点位移阵列

function P= Bar2D2Node_Forces(KK,q)

q=zeros(8,1)

q(num)=u

P=KK*q

至此，基于Matlab的杆单元有限元分析的程序设计已经完成，遇到实际问题时可以直接调用这些函数就可以解决问。

经典算例

如图所示的结构，各个杆的d性模量和横截面积都为 ,  。试基于MATLAB平台求解该结构的节点位移、单元应力以及支反力。

四杆桁架结构

解答：对该问题进行有限元分析的过程如下

（1） 结构的离散化与编号

对该结构进行自然离散，节点编号和单元编号如上图所示

（2）计算各单元的刚度矩阵(基于国际标准单位)

输入d性模量E、横截面积A，各点坐标。然后分别针对单元1，2，3和4，调用4次Bar2D2Node_Stiffness，就可以得到单元的刚度矩阵。

对应的主程序中代码：

E=2.95e11A=0.0001x1=0y1=0x2=0.4y2=0x3=0.4y3=0.3x4=0y4=0.3

k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2)

k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3)

k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3)

k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3)

（3） 建立整体刚度方程

由于该结构共有4个节点，因此，设置结构总的刚度矩阵为KK(8×8)，先对KK清零，然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。相关主程序代码为：

KK=zeros(8,8)

KK=Bar2D2Node_Assembly (KK,k1,1,2)

KK=Bar2D2Node_Assembly (KK,k2,2,3)

KK=Bar2D2Node_Assembly (KK,k3,1,3)

KK=Bar2D2Node_Assembly (KK,k4,4,3)

（4）边界条件的处理及刚度方程的求解

由图可以看出，被约束的自由度有：节点1的x,y方向自由度，节点2的y方向自由度，4节点的x、y方向两个自由度。则活动自由度编号为3，5，6.活动自由度对应的节点载荷F3=20000N，F5=0N，F6=25000N，采用高斯消去法进行求解，对应的代码为：

num=[3,5,6]%可活动的自由度编号

p=[200000-25000]

u=Bar2D2Node_Disp(KK,num,p)

（5）支反力的计算

在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。这部分对应的主程序的代码如下：

q=zeros(8,1)

q(num)=u%节点位移阵列

P=Bar2D2Node_Forces(KK,q)

（6）单元应力的计算

先从整体位移列阵q中提取出单元的位移列阵，然后，调用计算单元应力的函数Bar2D2Node_ElementStress，就可以得到各个单元的应力分量。

u1=[q(1)q(2)q(3)q(4)]

stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1)

u2=[q(3)q(4)q(5)q(6)]

stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2)

u3=[q(1)q(2)q(5)q(6)]

stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3)

u4=[q(7)q(8)q(5)q(6)]

stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)

（7）计算结果的整理

通过主程序的运行得计算结果。

主程序

%计算各单元的刚度矩阵(以国际标准单位)

E=2.95e11

A=0.0001

x1=0

y1=0

x2=0.4

y2=0

x3=0.4

y3=0.3

x4=0

y4=0.3

k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2)

k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3)

k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3)

k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3)

%建立整体刚度方程

%由于该结构共有4个节点，因此，结构总的刚度矩阵为KK(8×8)，先对K清零，然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。

KK=zeros(8,8)

KK=Bar2D2Node_Assembly (KK,k1,1,2)

KK=Bar2D2Node_Assembly (KK,k2,2,3)

KK=Bar2D2Node_Assembly (KK,k3,1,3)

KK=Bar2D2Node_Assembly (KK,k4,4,3)

%边界条件的处理及刚度方程求解

num=[3,5,6]%可活动的自由度编号

p=[200000-25000]

u=Bar2D2Node_Disp(KK,num,p)

%支反力的计算

q=zeros(8,1)

q(num)=u%节点位移阵列

P=Bar2D2Node_Forces(KK,q)

%各单元的应力计算

u1=[q(1)q(2)q(3)q(4)]

stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1)

u2=[q(3)q(4)q(5)q(6)]

stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2)

u3=[q(1)q(2)q(5)q(6)]

stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3)

u4=[q(7)q(8)q(5)q(6)]

stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)

说的可能有些罗嗦，注意其中有5个function，和最后一个主程序，计算的时候直接运行主程序就可以了。希望能帮助到你。

---