双线性四边形等参单元有限元程序

4
计算原理
4.1 形状函数 4 节点 4 边形等参元采用双线性插值函数，对如图 6 所示的母单元，插值函 数为：
图 6 4 节点 4 边形母单元
Ni =
ຫໍສະໝຸດ Baidu
1 (1 + ξ 0ξ )(1 + η0η ) 4
我们将位移表示为节点位移的插值形式，用矩阵写为：
[u ] = [ N ][ a ] (1)
其中： [u ] = [u , v]
(12)
将式(11)、(12)代入式(10)得：
[ε ] = [ A][G ][ a ] = [ B ][ a ] [σ ] = [ D ][ B ][ a ] (13) (14)
到此我们已将应变表示为单元节点位移的函数。进一步我们有应力矩阵： 这里[D]为弹性矩阵。 将(13)和(14)代入(8)式得：
Ue =
1 1 h[a]T ∫ ∫ [ B]T [ D][ B ] J dξ dη[a ] = [a]T [K e ][a] 2 2 −1 −1
1 1
(15)
这里单元刚度矩阵[Ke]为：
[ K e ] = h ∫ ∫ [ B]T [ D][ B] J dξ dη
−1 −1
1 1
(16)
至 此 ， 单 元 刚 度 矩 阵 的 计 算 已 经 完 成 。 函 数
∂y ⎤ ∂ξ ⎥ ⎥ ∂y ⎥ ∂η ⎥ ⎦
(4)
方程（3）逆过来写为：
⎛ ∂f ⎞ ⎛ ∂f ⎞ ⎜ ∂ξ ⎟ ⎜ ∂x ⎟ ⎟ ⎜ ⎟ = [ J ]−1 ⎜ ⎜ ∂f ⎟ ⎜ ∂f ⎟ ⎜ ∂y ⎟ ⎜ ∂η ⎟ ⎝ ⎠ ⎝ ⎠
或者：
(6)
⎛ ∂f ⎞ ⎛ ∂f ⎞ ⎜ ∂x ⎟ 1 ⎡ J − J12 ⎤ ⎜ ∂ξ ⎟ ⎟ ⎜ ⎟ = ⎢ 22 ⎥⎜ ⎜ ∂f ⎟ J ⎣ − J 21 J11 ⎦ ⎜ ∂f ⎟ ⎜ ∂y ⎟ ⎜ ∂η ⎟ ⎝ ⎠ ⎝ ⎠
界(一个单元的一条边算一个边界条件)所在的单元号，a(i,2)为第 i 个边界所在的 边在单元中的编号(编号规则如图 4 所示)。a(i,3)为第 i 个边界起始端点（按逆时
③ ④ ②
① 图 4 边界编号示意图
针方向）上的 x 方向的的压力，a(i,4)为第 i 个边界起始端点（按逆时针方向）上 的 y 方向的的压力，a(i,5)为第 i 个边界末端点（按逆时针方向）上的 x 方向的的 压力，a(i,6)为第 i 个边界末端点（按逆时针方向）上的 y 方向的的压力。如图 5 所示两个单元的结构，在 26 边受均匀横向正压力 p，则其载荷矩阵为：
图 6 程序流程图 pP= loadF(numNote,locNote,nF,h);
引入位移边界条件消除奇异
[K,P]=modifyK(pP,kK,nF,dB);
求解线性方程组
u=K\P; 节点位移输出 u=outU(u)
节点应力值
strNote=genNoteStress(numNote,locNote,nE,nF,E,v,p,u)
4.2 单元刚度矩阵及整体刚度矩阵的集成
(7)
4 节点 4 边形单元刚度矩阵从单元应变能而来：
1 U e = h ∫ [σ ]T [ε ]dA 2 e (8)
h 为结构厚度。 应变矩阵可以表示为位移的函数：
⎛ ∂u ⎞ ⎜ ⎟ ⎛ ε x ⎞ ⎜ ∂x ⎟ ⎜ ⎟ ⎜ ∂v ⎟ [ε ] = ⎜ ε y ⎟ = ⎜ ⎟ ⎜ ⎟ ⎜ ∂y ⎟ ⎝ γ xy ⎠ ⎜ ∂u ∂v ⎟ ⎜ + ⎟ ⎝ ∂y ∂x ⎠ 采用前面的微分转化公式(7)将[ε]表示为局部坐标的微分形式： ⎛ ∂u ⎞ ⎜ ∂ξ ⎟ ⎜ ⎟ ⎜ ∂u ⎟ ⎜ ∂η ⎟ ⎟ [ε ] = [ A] ⎜ ⎜ ∂v ⎟ ⎜ ∂ξ ⎟ ⎜ ⎟ ⎜ ∂v ⎟ ⎜η ⎟ ⎝ ⎠ 这里：
单元最佳应力点上的应力值
strCent=genCentStress(numNote,locNote,nE,nF,E,v,p,u)
3
主要函数简介 变量说明：locNote 为节点坐标，KE 为单元刚度矩阵，K 为总刚度矩阵，
numNote 为单元节点号逆时针，nE 为单元数，nF 为节点自由度数，E 为弹性模 量，v 为泊松比，p 为平面应力和平面应变的标志，p=1 为平面应力，p=0 为平 面应变，h 待求结构厚度。
0 1 +η 0 0 1+ ξ 0 1 −η 0 1+η −1 − ξ 0 1+ ξ
−1 − η 1−ξ 0 0
⎡u1 ⎤ ⎢ ⎥ ⎢v1 ⎥ 0 ⎤ ⎢u 2 ⎥ ⎥ ⎢v ⎥ 0 ⎥ ⎢ 2 ⎥ =[G][a] −1 − η ⎥ ⎢u3 ⎥ ⎥⎢ ⎥ 1 − ξ ⎦ ⎢v3 ⎥ ⎢u ⎥ ⎢ 4⎥ ⎢v4 ⎦ ⎥ ⎣
一 4*4 的矩阵：
⎛1 ⎜ ⎜2 ⎜5 ⎜4 ⎝
2 4 5⎞ 3 5 6⎟ 6 8 9⎟ 5 7 8⎠
⎟ ⎟
图 2 输入说明
“ noteLocation.txt ” 为 节 点 坐 标 数 据 ， a(i,j) 为 节 点 i 的 第 j 个 坐 标值，j=1 即为 X 坐标，j=2 即为 Y 坐标。如图 3 所示：
一个典型的计算结果如图 6 所示：
图 6 输出结果
其中 u 为节点位移，u(1,i)为第 i 个节点 x 方向的位移，u(2,i)为第 i 个节点 y 方向的位移。strNote 为节点应力，strNote(1,i)为第 i 个节点的 σx，strNote(2,i) 为 第 i 个节点的 σy， strNote(1,i)为第 i 个节点的 σxy。 StrCent 为单元最佳应力点应力， strCent(1,i)为第 i 个单元的 σx，strCent(1,i)为第 i 个单元的 σy，strCent(1,i)为第 i 个单元的 σxy。 2 程序结构说明 本程序的结构和一般的有限元程序没什么不同，只是没有必要的前处理过 程，需要用户手动输入。程序流程图如图 6 所示： 主程序 main()调用过程如下： 输入数据加载：
[E,v,h,p,locNote,numNote,dB]=dateInput()加载输入数据，返回弹性模量 E，泊松比 v 等。 kK=assembleK(numNote,locNote,nE,nF,E,v,h,p)为总刚度矩阵集成。包括单元刚度矩阵计 算 子 函 数 elemK=genKE(numNote,locNote,numElement,E,v,h,p) 和
图 3 节点坐标数据示例
“noteNum.txt”为单元节点号矩阵，a(i,j)为第 i 个单元的第 j 个节点的总体 编号。 单元节点按照逆时针编号。 对于四节点四边形单元该矩阵为 nE*4 的矩阵， nE 为单元数！如图 2 所示。 “displacementboundary.txt”为位移边界矩阵，该矩阵为 nF*2 的数组，nF 为 总自由度数，dB(i,1)为第 i 个节点是否位移边界判据，1 为是，0 为否,dB(i,2)为 第 i 个节点常位移值，如果该节点是位移边界则为位移值，否则为 0。 “surface load.txt”为表面载荷矩阵，该矩阵为 n*6 的矩阵。a(i,1)为第 i 个边
⎡ J 22 1 ⎢ [ A] = ⎢ 0 J ⎢ ⎣ − J 21 − J12 0 J11 0 − J 21 J 22 0 ⎤ J11 ⎥ ⎥ − J12 ⎥ ⎦
(9)
(10)
(11)
进一步，我们有：
⎛ ∂u ⎞ ⎜ ∂ξ ⎟ ⎜ ⎟ 1 −η ⎡η − 1 0 ⎜ ∂u ⎟ ⎜ ∂η ⎟ 1 ⎢ξ − 1 0 −1 − ξ ⎜ ⎟= ⎢ 0 ⎜ ∂v ⎟ 4 ⎢ 0 η − 1 ⎢ ⎜ ∂ξ ⎟ 0 ⎣ 0 ξ −1 ⎜ ⎟ ⎜ ∂v ⎟ ⎜η ⎟ ⎝ ⎠
双线性四边形等参单元有限元程序
本程序采用 matlab 编写。程序加载由用户提供的前处理数据，包括网格数据 和载荷数据。采用直接的数值运算和 matlab 符号运算两种方法（可选择）生成 单元刚度矩阵。自动集成结构刚度矩阵，选用直接解法求解线性方程组，解出节 点位移。后处理过程中，程序计算了节点应力值，对于共享节点的应力程序采用 各个单元计算值的平均，程序同时给出了单元最佳应力点的应力值！ 1 使用说明 用户需要在目录中给定的文件中按照既定的格式给出必要的前处理数据（后 面有详细说明） 。然后将 matlab 的工作目录设为我们提供的“双线性四边形等参 单元程序”目录，然后在命令行中输入 main()并按下回车键即可。如图 1 所示：
elemK=genKE2(numNote,locNote,numElement,E,v,h,p)。返回结构总刚度矩阵。 elemK=genKE(numNote,locNote,numElement,E,v,h,p)采用直接的数值计算的方法求解单 元刚度矩阵，并返回。 elemK=genKE2(numNote,locNote,numElement,E,v,h,p)采用符号运算求解单元刚度矩阵， 并返回。 p=loadF(numNote,locNote,nF,h)计算总体载荷列阵，并返回。 [K,P]=modifyK(pP,kK,nF,dB)采用乘大数的方法引入位移边界条件，消除刚度矩阵奇异 性，返回修正后的刚度矩阵和载荷列阵。 strNote=genNoteStress(numNote,locNote,nE,nF,E,v,p,u)计算节点应力， 采用围绕该节点的 单元计算应力的平均值作为该节点的应力值。 strCent=genCentStress(numNote,locNote,nE,nF,E,v,p,u) 计算单元最佳应力点即局部坐标 （0，0）点的应力值。
5 ② 3 ① 1 图 5 载荷矩阵示例图 6
4
2
[1 2 p 0 p 0; 1 2 p 0 p 0]。 “concentrated force.txt”为集中力载荷，为 1*nF 的矩阵，a(i)为第 i 个自由 度上的集中力。 另外程序在运行过程中需要输入材料的弹性模量和泊松比，结构厚度，还需 要选择平面应力或者平面应变，只需要按照程序的说明输入即可！ 注意：本程序没有提供容错处理程序，用户必须自己确认输入无误！ 1.2 输出数据
4
y = ∑ Ni * yi
i =1
(2)
等参单元中插值函数是基于母单元的局部坐标的，因此需要进行局部坐标下 微分与整体坐标微分的转化，即：
⎛ ∂f ⎞ ⎛ ∂f ⎞ ⎜ ∂ξ ⎟ ⎜ ⎟ ⎜ ⎟ = [ J ] ⎜ ∂x ⎟ ⎜ ∂f ⎟ ⎜ ∂f ⎟ ⎜ ∂y ⎟ ⎜ ∂η ⎟ ⎝ ⎠ ⎝ ⎠
(3)
图1
程序运行说明
1.1 输入数据 程序的输入数据需要在“双线性四边形等参单元程序”目录中给定的文件中 完成。每一个.txt 文件对应一个矩阵，所以数据的输入必要严格按照 matlab 矩阵 加载文件的格式完成，即文件的一行对应于矩阵行，数据之间用空格或者逗号隔 开，分号或者换行符表示进入矩阵另一列的输入。如图 2：所示输入的矩阵阵为
T
[a ] = [u1 , v1 , u2 , v2 , u3 , v3 , u4 , v4 ]T ⎡ N1 0 N 2 0 N 3 0 N 4 0 ⎤ [N ] = ⎢ ⎥ ⎣ 0 N1 0 N 2 0 N 3 0 N 4 ⎦
我们选用相同的插值函数做等参变换，即：
x = ∑ Ni * xi
i =1 4
[E,v,h,p,locNote,numNote,dB]=dateInput();
总自由度数计算
nF=2*size(locNote,2);
单元数计算
nE=size(numNote,2);
总刚度矩阵计算
kK=assembleK(numNote,locNote,nE,nF,E,v,h,p);
载荷列阵计算
[J]为雅可比矩阵：
⎡ ∂x ⎢ ∂ξ [J ] = ⎢ ⎢ ∂x ⎢ ∂η ⎣
将（2）式代入得：
J= 1 ⎡(x 3 -x 4 )(1+η )+(x 2 -x1 )(1-η ) (y 2 -y1 )(1-η )+(y3 -y 4 )(1+η ) ⎤ ⎡ J11 ⎥=⎢ 4⎢ ⎣ (x 4 -x1 )(1-ξ )+(x 3 -x 2 )(1+ξ ) (y 4 -y1 )(1-ξ )+(y3 -y 2 )(1+ξ ) ⎦ ⎣ J 21 J12 ⎤ J 22 ⎥ ⎦ (5)