基于MATLAB平台的悬臂梁静力弹性分析

谭 宸

(同济大学土木工程学院,上海 200092)

0 引言

有限单元法是当今工程分析中获得广泛应用的数值方法[1],采用有限元方法可以求解具有复杂的几何边界的弹性问题。曹子龙等[2]建立了求解圆孔薄板弹性问题的杂交应力单元,推导出圆孔薄板问题的极坐标下应力插值矩阵,可较高精度的求解孔边附近应力。陈立胜等[3]采用有限单元法初步对弹性平面问题进行了初步尝试。张明哲等[4]开发的计算软件平台有助于学生很好的理解在弹塑性问题中改进无单元Galerkin方法的应用。目前,伴随着有限元软件功能的不断完善,对于刚接触有限元软件的学生而言,往往对其一知半解,在有限元分析中着重于物理参数的调节[5]。

为了加深学生们对于有限元方法基本原理的理解,本文通过基于MATLAB平台,编写结构化的有限元代码,可以应用于基本的结构工程问题。该代码的结构易于理解,并可由学生轻松扩展。本文开发的代码计算悬臂梁变形结果与商业有限元软件ABAQUS计算结果一致,开发的一个结构良好、易于使用的计算悬臂梁结构弹性变形的有限单元分析程序具有很好的适用性。

1 弹性问题的有限单元格式

工程中许多问题,通常是以未知场函数在域内应满足的偏微分方程形式提出。根据弹性力学基本方程和与之等效的变分原理,可以从基本方程的强形式推导出有限元格式,进而可获得局部单元刚度矩阵。

本小节以二维平面域中弹性问题为例,考虑边界条件,对其偏微分方程的有限元格式进行推导。从面积域中取一无限小单元,其受力示意图如图1所示。

从图1可以看出,该单元的平衡方程可表示为:

(1)

(2)

用笛卡尔张量符号可以进一步将弹性力学平衡方程改写为:

σji,j+bi=0 (i,j=1,2)

(3)

通过加权余量法可写出平衡微分方程的等效积分弱形式。等效积分弱形式可以通过分部积分得到:

(4)

作为平衡微分方程的等效积分弱形式,在导出过程中并未涉及物理方程,因此,不仅仅可以适用于线弹性问题,还可以用于非线性弹性或者弹塑性问题。

将权函数wi替换成δui,可以得到有限单元法的控制方程:

(5)

改写成矩阵形式可得:

(6)

本文采用4结点矩形单元来处理二维平面弹性力学问题,4结点矩形单元如图2所示。

对于每个角点,构造出它的形函数为:

(7)

进一步,根据弹性力学中位移与应变的关系,可将单元应变用结点位移向量来进行表示:

{ε}=εxεyγxyìîíïïïïüþýïïïï=N1(x,y)x0N2(x,y)x0N3(x,y)x0N4(x,y)x00N1(x,y)y0N2(x,y)y0N3(x,y)y0N4(x,y)yN1(x,y)yN1(x,y)xN2(x,y)yN2(x,y)xN3(x,y)yN3(x,y)xN4(x,y)yN4(x,y)xéëêêêêêêêùûúúúúúúúüþýïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïïï[B] u1 v1 u2 v2 u3 v3 u4 v4u1v1u2v2u3v3u4v4ìîíïïïïïïïïïïïüþýïïïïïïïïïïï(8)

将形函数代入可得梯度矩阵[B]:

(9)

局部单元刚度矩阵K:

(10)

代入B矩阵,进而可以得到:

[K]= u1 v1 u2 v2 u3 v3 u4 v4k11k12k13k14k15k16k17k18k21k22k23k24k25k26k27k28k31k32k33k34k35k36k37k38k41k42k43k44k45k46k47k48k51k52k53k54k55k56k57k58k61k62k63k64k65k66k67k68k71k72k73k74k75k76k77k78k81k82k83k84k85k86k87k88éëêêêêêêêêêêêùûúúúúúúúúúúúu1v1u2v2u3v3u4v4(11)

其中,

。

2 程序架构

该程序基于MATLAB编程语言开发,依托于有限单元法解决弹性问题,可以高效便捷的对悬臂梁受集中荷载作用下的应力应变进行分析。程序通过先对几何域进行单元网格划分,建立节点和对单元进行编号,随后输入边界条件等模型数据,对其进行采用有限单元法的弹性分析后得到处理结果,最后将结果云图以可视化的界面进行展示。程序可以分为前处理程序、有限元弹性分析计算程序和可视化后处理程序三个部分,其流程图和MATLAB函数如图3所示。

前处理程序主要是通过建立几何模型,生成节点和网格,将分析域离散化,组建全局坐标系统,同时引入边界条件。输入参数后预处理必备的结构参数值,前处理关键环节需要定义节点、生成每个单元坐标信息、输入载荷信息和边界条件,为有限元法分析做准备。计算分析程序进行有限单元法的计算,最主要的是进行矩阵计算和求解。计算分析程序通过求解全局刚度矩阵获得每个节点的位移{u},进而评估每个单元上的应力和应变。后处理程序模块对每个单元建立节点位移和应力向量,进而绘制加载后模型的位移、应力云图,实现计算结果的可视化。

2.1 前处理程序

前处理程序是数值计算中非常重要的环节,影响后续数值计算的精度和效率。前处理可分为如下部分:输入材料形状大小、材料性质参数、生成节点坐标等单元信息以及载荷和边界条件定义等。

前处理程序的首要是生成节点坐标,包括节点的编号及相应的点的坐标值。节点位置信息定义为函数square_node_array,生成节点信息后,需要获取每个单元的节点信息。节点的空间信息存储在element数组中。为便于后续边界条件的处理,方便使用,按照单元中节点逆时针编号进行存储记录,在完成节点的计算、单元节点信息的输入后,编写supportcond函数引入边界条件。对于悬臂梁在端部受集中荷载而言,其荷载位移边界条件的matlab代码如下:

function [topEdge,topEdge1,dispNodes,dispNodes1,leftNodes1]=supportcond(numx,numy)

nnx=numx+1;

nny=numy+1;

uln=nnx*(nny-1)+1;

urn=nnx*nny;

lrn=nnx;

lln=1;

topEdge=[uln:1:(urn-1);(uln+1):1:urn]′;

topEdge1=topEdge;

botEdge=[lln:1:(lrn-1);(lln+1):1:lrn]′;

rightEdge=(lrn:nnx:(urn))′;

botNodes=unique(botEdge);

topNodes=unique(topEdge);

rightNodes=unique(rightEdge);

leftNodes=rightNodes-(nnx-1);

dispNodes=botNodes;

rightNodes1=rightNodes(2:end);

leftNodes1=leftNodes(2:end);

dispNodes1=leftNodes;

end

2.2 有限元计算分析程序

其中,B矩阵的matlab代码如下:

function Bfem=Bmatrix(pt,iel)

global node element

sctr=element(iel,:);

nn=length(sctr);

[N,dNdxi]=shape_func(pt);

J0=node(sctr,:)′*dNdxi;

invJ0=inv(J0);

dNdx=dNdxi*invJ0;

Bfem=zeros(3,2*nn);

Bfem(1,1:2:2*nn)=dNdx(:,1)′;

Bfem(2,2:2:2*nn)=dNdx(:,2)′;

Bfem(3,1:2:2*nn)=dNdx(:,2)′;

Bfem(3,2:2:2*nn)=dNdx(:,1)′;

end

2.3 后处理可视化程序

后处理程序获取有限元计算分析程序得到的位移等信息,并重新组织成标准后处理结果文件的形式,利用matlab可视化技术再现网格信息,将变形前后的网格展示出来。同时,分析计算结果,如位移云图、应力云图的结果以图像形式显示出来。本系统的后处理程序包括数据处理程序、图像显示程序。

3 数值算例计算与验证

3.1 算例计算

对二维受集中力的悬臂梁采用弹性问题的有限单元法进行小变形数值分析。通过对该悬臂梁进行自主编程建模,建模核心代码流程如第2小节所示。利用前处理程序模块在梁内布置节点,进行单元编号,记录每个单元的节点信息。将输入节点信息和荷载位移边界条件输入到有限元分析计算程序中,进而计算获得节点的位移和应力、应变等数据。将所得结果与ABAQUS软件分析结果进行比较,验证了本文提出的弹性问题有限单元法程序的有效性,可进一步加深学生们采用有限元法分析问题编制程序的理解。

悬臂梁自由端受集中荷载示意图如图4所示。梁的几何尺寸为L=10 m,D=1 m,梁的厚度为t=0.1 m。梁受集中荷载P=0.1 N。材料的弹性模量为E=100 kPa,泊松比为v=0.25,按照平面应力问题进行求解。

采用前处理程序进行结点生成,单元网格划分,水平划分单元数为100,竖向划分单元数为10,结果如图5所示。

将节点信息和荷载条件输入到有限元分析计算程序中进行计算,获得变形后的悬臂梁有限单元网格如图6所示。

悬臂梁的位移云图如图7所示。

从图7中可以看出,右端部顶处竖向结点位移为0.039 8 m。水平应力σxx与竖向应力σyy云图见图8。

3.2 方法验证

基于ABAQUS有限元的数值模拟方法是检验本文程序编制的正确性的有效途径之一。本小节通过建立与3.1算例相同的悬臂梁结构有限元模型,并对模型赋予材料参数、施加集中荷载边界条件,验证本文方法的正确性。

在ABAQUS中,建立的悬臂梁三维结构几何参数为L=10 m,D=1 m,t=0.1 m。悬臂梁材料弹性模量E=100 kPa,泊松比为v=0.25。有限元模型如图9所示。在有限元模型中,所施加的载荷为力载荷。模型网格划分如图10所示。

ABAQUS计算获得的悬臂梁的位移云图如图11所示,从图11中可以看出,右端部顶处竖向结点位移为0.040 m。水平应力σxx与竖向应力σyy云图见图12。

通过将本文编写的悬臂梁有限单元法程序与ABAQUS计算的结果进行对比,可以看出,所开发的基于MATLAB平台的悬臂梁弹性分析程序可以较好的预测梁的小变形,所得计算结果与有限元软件计算结果一致。

4 结论

本文以求解悬臂梁受集中荷载作用下变形为算例,介绍了有限单元法求解弹性问题的算法。通过将商业有限元分析软件结果与本文程序编制结果对比,两者一致,验证了本文程序的有效性。该程序基于MATLAB平台编制而成,程序简单实用,计算效率比采用有限元软件ABAQUS分析弹性问题的效率更高。本文展示的有限元求解弹性问题的方法通用性强,通过该程序,可以进一步加深学生们采用有限元法分析问题编制程序的理解,具有重要意义和参考价值。通过对该程序的学习,可以举一反三,在该程序的基础上修改荷载和位移边界条件,进一步求解其他的弹性问题。