第三章 线弹性静力学问题.docx
《第三章 线弹性静力学问题.docx》由会员分享，可在线阅读，更多相关《第三章 线弹性静力学问题.docx（36页珍藏版）》请在冰豆网上搜索。
第三章线弹性静力学问题
第三章线弹性静力学问题
第一节用向导AppWizard解线弹性静力学问题
本例将求解一个简支的深梁上表面受均布载荷作用下的平面应变问题，求出梁的变形和应力分布。
图3.1.1受均布压力的简支深梁
该问题的偏微分方程表达式及定解条件如下：
1）平衡方程
2）几何方程
3）本构方程
4）边界条件
第二类边界条件
具体操作步骤如下：
．设定工作路径。
点击File菜单下的WorkDir…选项，在弹出的对话框中填入工作路径。
如图3.1.2所示
图3.1.2设定路径对话框
．点击AppWizard菜单，弹出公式库列表，系统提供了数学、固体力学、岩土、热传导、渗流、流体力学等领域的各种问题的公式。
如图3.1.3所示。
点击solid，在ProjectName中填入“S”；点击“Next”进入下一步。
图3.1.3选择学科领域
．选择坐标系。
点击左边的“2dxy”，选中二维直角坐标系，如图3.1.4所示。
再点击“Next”进入下一步。
图3.1.4选择坐标系
．选择问题类型，即该问题是线弹性还是弹塑性。
点击“elastic”，选择线弹性，如图3.1.5所示。
点击“Next”进入下一步。
图3.1.5选择问题类型
．选择方程类型，即方程是平面应力问题还是平面应变问题，是静态问题还是动态问题；对于动态问题还可选择Newmark格式求解。
点击“Static_Plane_Strain_Problem”，选择选择静态平面应力问题的方程，如图3.1.6所示。
点击“Next”进入下一步。
图3.1.6选择方程类型
．选择单元类型。
如图3.1.7所示。
点中“q4”，再点击中间的“Add”按钮，加入到右边方框表示采用四节点矩形单元。
点击“Next”进入下一步。
图3.1.7选择单元类型
．选择边界单元。
对于本算例，因为是应力边界条件，所以需要边界单元。
点中“ll2”，再点击中间的“Add”按钮，加入到右边方框表示采用两节点线单元。
如图3.1.8所示。
点击“Next”进入下一步。
图3.1.8选择边界单元
．选择求解器。
系统提供有：
对称求解器（sin）、非对称求解器（nin）、Gauss-Seidel迭代求解器（gs）、逐次超松弛迭代求解器（sor）、共轭梯度法求解器（cgm）、不完全LU分解求解器（ilu）等。
点击“sin”，选择对称求解器。
如图3.1.9所示。
点击“Next”进入下一步。
图3.1.9选择求解器
．选择存储方式。
对于有限元程序自动生成系统网络版（IFEPG）仅提供用外存的存储方式。
对FEPG3.0提供外存和内存两种存储方式。
．选择计算数据。
该数据给出求解区域的有限元数据，如图3.1.10所示。
用户应注意：
点击左边文件列表选择数据时，看一下右边方框中说明的方程类型和单元类型与你前面选择的是否一致，不一致系统将给出提示，一致时才可以正确运行。
点击“sdispe_q4”，点击“Next”进入下一步。
图3.1.10选择有限元计算数据
．自动生成全部有限元程序。
点击“Run”即可生成全部有限元程序。
如图3.1.11和图3.1.12所示。
图3.1.11点击“Run”自动生成
图3.1.12全部有限元程序
计算结果
插入res01.bmp
图3.1.13算例的有限元计算结果
第二节改变施加的载荷
生成全部有限元计算程序之后，在FEPG界面的左边是WorkSpace列表，它将所有生成的文件以及生成这些文件所需的原始文件都列在其中。
FEPG施加荷载的方式有三种，一是分布荷载，作为材料属性给出，二是面力，通过边界单元材料属性给出，三是节点力，通过节点给出。
改变施加的载荷需修改前处理中的DAT文件，以下分别按三种加载分别叙述具体如何修改。
a.改变分布力的大小
在WorkSpace列表中找到“preprocessor”文件夹，在修改文件DAT文件之前先看一看MTI文件，该文件定义了有限元计算中需要的所有数据表格。
单击“preprocessor”文件夹中的s.mti文件，在界面的右边则显示出该文件。
找到mate段如下
mate
npepvfxfy
i4E10.4E10.4E10.4E10.4
mate定义为材料参数表格名；
n表示材料类型号；
pe表示弹性模量；
pv表示poisson比；
fx表示x方向的分布力；
fy表示y方向的分布力；
在了解了mate定义的材料参数意义后，再看DAT文件。
单击“preprocessor”文件夹中的s.dat文件编辑该文件，找到mate信息段内容如下
mate
1;1.0e10;0.25;0.;0.;
对应与s.mti文件，我们可以知道
1表示第一种材料，如果有多种材料可在下一行写“2”及相应的材料参数；
1.0e10表示材料的弹性模量的值；
1.25表示材料的泊松比的值；
第一个0.表示x方向分布力的大小；
第二个0.表示y方向分布力的大小；
弄清各个数值的意义我们就可以对它作相应的修改了，比如假设在x方向
有1.0e3大小的力则mate信息段改成如下形式
mate
1;1.0e10;0.25;1.0e3;0.;
下面看一下改变了载后的计算结果。
首先在s.dat上单机鼠标右键弹出如下命令
mtisdat
该mti的两个输入参数分别是s.mti和s.dat两个文件，输出是前处理程序prg.for和前处理显示程序prt.for。
运行上述命令重新生成prg.for和prt.for。
然后运行run_batch_file文件夹中的s.bat重新的计算结果如下
插入res02.bmp
b.改变面力
同样，我们先看一下s.mti文件。
在s.mti文件中找到bf段如下
bf
nfxfy
i4E10.4E10.4
bf定义为面力的参数表格名
n表示面力类型
fx面的切向方向力的分量的值
fy面的法向方向力的分量的值（参考局部标架说明）
再看s.dat中对应的行，找到bf信息段如下
bf
1;0.0;-1.0e6;
1表示第一种类型的面力
0.0表示切向方向上的力为零
-1.0e6表示法向方向上的力为1.0e6，方向和法向相反。
下面给出两种不同的面力，一是将-1.0e6改为1.0e6，这表示算例中的梁是受向上的拉力，同上面所述保存s.dat文件并用mti命令重新生成prg.for和prt.for，再运行s.bat重新进行计算，得计算结果如下
插入res03.bmp
二是将bf信息段改为如下形式
bf
1;1.0e6;0.0;
即在切向方向加上剪力，方向方向不加力。
同样重新生成prg.for并进行计算。
计算结果如下
插入res04.bmp
c.施加节点力
节点力的修改同样需在前处理中给出，让我们看一看“preprocessor”中的s.mti文件和s.dat文件。
先看s.mti中id段和disp段
idf
niduidv
i6i6i6
id表示该段定义约束信息
n表示节点编号
idux方向的约束信息存入idu数组
idvy方向的约束信息存入idv数组
dispf
nuv
i6e10.4e10.4
disp表示该段定义给定位移或荷载信息
n表示节点编号
ux方向的给定位移或荷载信息存入u数组
vy方向的给定位移或荷载信息存入v数组
再看s.dat文件中id信息段和disp信息段
id
&nn=（nx+1）*（ny+1）;
[i;1;1;i=1,nn]
1;-1;-1;
nx+1;1;-1;
&n=nn;
&nn=（nx+1）*（ny+1）定义总节点数
[i;1;1;I=1,nn]以配对的方括号表示循环，即表示先把所有节点的约束赋值为1。
FEPG系统规定1表示节点自由，-1表示节点被约束。
1;-1;-1;表示1节点在x和y两个方向都被约束住。
nx+1;1;-1;表示nx+1节点仅在y方向被约束住，x方向是自由的。
&n=nn表示把最大节点后赋值给n，用户应注意加上此行，因为这将使生成的prg.for中把约束信息写入文件的循环的上限为最大节点号，也就是保证将所有的节点约束信息写入文件。
再看disp信息段
disp
[i;0.0;0.0;i=1,（nx+1）*（ny+1）]
disp表示该段定义给定位移和荷载信息
[i;0.0;0.0;i=1,（nx+1）*（ny+1）]以配对的方括号表示循环，将所有节点在x和y两个方向的位移或荷载信息都赋零。
FEPG系统作如下判断：
如果某节点的某个方向的约束信息是-1，即该节点在该方向被约束住，则在disp中给出的相应的值表示给定位移；如果某节点的某个方向的约束信息是1，即该节点在该方向是自由的，则在disp中给出的相应的值表示给定荷载。
为了更清楚地看到施加节点荷载的结果，先将分布力和面力各方向上的荷载值置零，修改后的mate信息段和bf信息段如下
mate
1;1.0e10;0.25;0.0;0.;
bf
1;0.0;0.0;
在右上角节点沿水平方向向右施加大小为1.0e4的集中力，修改后的disp信息段如下
disp
[i;0.0;0.0;i=1,（nx+1）*（ny+1）]
（nx+1）*（ny+1）;1.0e4;0.0;
用mti命令重新生成前处理程序再计算后的如下结果
插入res05.bmp
在右上角节点沿垂直方向向上施加大小为1.0e4的集中力，修改后的disp信息段如下
disp
[i;0.0;0.0;i=1,（nx+1）*（ny+1）]
（nx+1）*（ny+1）;0.0;1.0e4;
用mti命令重新生成前处理程序再计算后的如下结果
插入res06.bmp
在左上角节点沿水平方向向左施加大小为1.0e4的集中力，修改后的disp信息段如下
disp
[i;0.0;0.0;i=1,（nx+1）*（ny+1）]
（nx+1）*ny+1;-1.0e4;0.0;
&n=（nx+1）*（ny+1）;
大家应注意此修改增加了一行赋值，即把总节点数赋给n，理由同在id中阐述的理由一样。
为何上两例没有给这样的赋值呢？
因为上两例都是给的最大节点号，所以n自然取得总节点数，否则也应有同样的赋值行。
荷载正负号的规定是和总体坐标一致为正，反之为负。
用mti命令重新生成前处理程序再计算后的如下结果
插入res07.bmp,u1放大50000被的变形图
第三节改变约束方式
上面我们已经提到id信息段给出约束信息，下面我们就举几个例子说明如何修改约束信息
假设物体左端固定，右端悬空，且在右上角受一大小为1.0e4向下的集中力，如下图所示。
对s.dat文件中id信息段做如下修改
id
&nn=（nx+1）*（ny+1）;
[i;1;1;i=1,nn]
[（i-1）*（nx+1）+1;-1;-1;i=1,ny+1]
&n=nn;
相应的节点荷载也作响应修改，修改后内容如下
disp
[i;0.0;0.0;i=1,（nx+1）*（ny+1）]
（nx+1）*（ny+1）;0.0;-1.0e4;
用mti命令重新生成前处理程序再计算后的如下结果
插入res08.bmp，变形放大5000倍的结果
假设物体左端固定，右端有一竖向支撑，且在右上角受一大小为1.0e4向下的集中力，如下图所示。
在上例的基础上只要增加右下角节点在y方向上的约束即可，修改后的id信息段如下
id
&nn=（nx+1）*（ny+1）;
[i;1;1;i=1,nn]
[（i-1）*（nx+1）+1;-1;-1;i=1,ny+1]
nx+1;1;-1;
&n=nn;
用mti命令重新生成前处理程序再计算后的如下结果
插入res09.bmp，变形放大50000倍的结果
最后一例假设物体两端固定，且在中间位置受一大小为1.0e4向下的集中力，如下图所示。
修改id信息段如下
id
&nn=（nx+1）*（ny+1）;
[i;1;1;i=1,nn]
[（i-1）*（nx+1）+1;-1;-1;i=1,ny+1]
[i*（nx+1）;-1;-1;i=1,ny+1]
&n=nn;
disp信息段也作如下修改
disp
[i;0.0;0.0;i=1,（nx+1）*（ny+1）]
（nx+1）*ny+（nx+2）/2;0.0;-1.0e4;
&n=（nx+1）*（ny+1）;
用mti命令重新生成前处理程序再计算后的如下结果
插入res10.bmp，变形放大50000倍的结果
第四节改变材料参数
前面在修改分布力是已经提到mate信息段，并解释了各参数的意义，其中前两个参数是线弹性材料的弹性模量个泊松比，下面就修改材料参数，看计算结果如何变化。
在两段固定，中间受以集中力的情况下，看不同的材料的变形情况。
下面分别给出在弹性模量E分别为1.0e9，5.0e9和1.0e10时mate信息段修改后的内容
mate
1;1.0e9;0.25;0.0;0.;
mate
1;5.0e9;0.25;0.0;0.;
mate
1;1.0e10;0.25;0.0;0.;
每次修改后保存s.dat文件，用mti命令重新生成前处理程序再计算后的得如下计算结果的比较
插入res11.bmp，变形放大5000倍的结果
第五节两种不同的材料
下面考虑两种不同材料时，如何用FEPG求解，假设求解下列两层材料深梁在集中力作用下的变形，如下图所示。
如图中的两层材料弹性模量
，
，泊松比
，
，厚度相同，右上角受垂直向下大小为1.0e4大小的力，分析在此情景下深梁的变形。
为用FEPG解算该问题，需在前处理中作相应的修改。
看一看“preprocessor”文件夹中s.mti文件的q4段，内容如下
q4
nnod1nod2nod3nod4mate
i6i6i6i6i6i6
q4表示四节点矩形单元
n表示单元编号
nod1表示单元第一节点编号存入的数组名
nod2表示单元第二节点编号存入的数组名
nod3表示单元第三节点编号存入的数组名
nod4表示单元第四节点编号存入的数组名
mate表示该单元对应的材料号
对于本算例需修改s.dat文件实现分层划分单元，上层mate的值给1，下层mate的值给2，mate对应的信息段分别给出两中材料的材料参数，修改后的具体内容如下。
对q4信息段作如下修改
q4
[[nx*（j-1）+i;（j-1）*（nx+1）+i;（j-\1）*（nx+1）+i+1;j*（nx+1）+i+1;j*（nx+1）+i;1;i=1,nx];j=1,ny/2]
[[nx*（j-1）+i;（j-1）*（nx+1）+i;（j-\1）*（nx+1）+i+1;j*（nx+1）+i+1;j*（nx+1）+i;2;i=1,nx];j=ny/2+1,ny]
mate信息段也需相应修改
mate
1;1.0e10;0.25;0.0;0.;
2;2.0e10;0.30;0.0;0.;
用mti命令重新生成前处理程序再计算后得如下结果
insertres12.bmpandres13.bmp变形放大5000倍的结果
第六节线弹性力学方程的PDE文件
FEPG通过用户填写的PDE文件生成单元子程序。
在了解如何填写PDE文件之前先回顾一下线弹性力学所对应的偏微分方程。
线弹性力学偏微分方程包括平衡方程、几何方程和本构方程。
在二维直接坐标系下，线弹性平面应变问题的方程具体可写成如下形式
平衡方程：
（1）
几何方程
（2）
本构方程
（3）
对
（1）第一个方程乘
，第二个方程乘
，将两式对体积积分并相加得
（4）
对（4）式分布积分可得
（5）
在考虑下列Cauchy公式
（6）
其中
和
分别表示应力边界上力矢量在x和y方向的分量。
将Cauchy公式和几何方程
（2）式代入（5）式得
（7）
将本构方程代入（7）式得最终的虚功方程的弱形式
（8）
其中
为一常系数
有了（8）式我们就可以开始填写PDE文件，先考虑（8）式中的体积积分项，由此构造体单元。
至于（8）中的边界积分项我们在下一节给出详细说明。
按照（8）式给出如下的PDE文件，参阅“sdispe”文件夹中的sdispe.pde文件
右边给出相应的注解
dispu,v,给出未知函数名U，V
coorx,y,给出总体坐标系下的坐标变量名X，Y
func=funa,funb,func,给出需要用到的函数名
Shap%1%2给出单元形状类型符和节点个数
gaus%3给出每个方向上积分点的个数或单元形状类型符
mass%1给出单元形状类型符
load=fufv给出载荷向量，对应于
$c6pe=prmt
（1）从PRMT数组中取出弹性模量，$c6表示插入FORTRAN源程序
$c6pv=prmt
（2）从PRMT数组中取出泊松比
$c6fu=prmt（3）从PRMT数组中取出X方向体力
$c6fv=prmt（4）从PRMT数组中取出Y方向体力
$c6fact=pe/（1.+pv）/（1.-2.*pv）由弹模和泊松比组合成系数FACT
Func给出一些函数的定义
funa=+[u/x]定义FUNA为
，[u/x]表示
funb=+[v/y]定义FUNB为
，[v/y]表示
func=+[u/y]+[v/x]定义FUNC为
stif给出刚度矩阵
dist=
+[funa;funa]*fact*（1.-pv）[funa;funa]的第一个FUNA
+[funa;funb]*fact*（pv）表示
，第二个FUNA表示
+[funb;funa]*fact*（pv）[；]表示积分。
+[funb;funb]*fact*（1.-pv）
+[func;func]*fact*（0.5-pv）
load=+[u]*fu+[v]*fv
关于PDE文件还需作如下说明
a.单元形状类型符，q表示矩形单元；
b.guas信息段后给的参数是数字表示每个方向上积分点的个数，如后跟单元形状类型符则表示采用节点积分；
c.prmt（）数组保存了我们前面所说的前处理数据中的材料参数，在PDE文件里通过赋值取出来。
在了解了PDE文件如何得来以及PDE文件各内容的意义后，让我们看一看如何由PDE文件生成单元子程序，又如何连接到单元计算的E程序中。
将sdispe.pde文件中的shap、gaus、mass信息段进行以下修改
shapq4
gausq
massq
其它内容不变。
上述修改表示将生成四节点矩形单元子程序，体积分采用节点积分。
保存sdispe.pde文件，单击FEPG界面主菜单中的“Run”，弹出FEPG命令对话框，在输入栏输入
pdegessdispeseuq4
其中pdeges为系统命令，其功能是由PDE文件生成GES文件和单元子程序，它的第一个参数为PDE文件主名，第二个参数为将要生成的GES文件也是单元子程序文件的主名。
第二个参数也可不写，此时表示生成的GES文件或单元子程序文件和PDE文件有相同的主名。
另PDE文件和GES文件的关系是：
将PDE文件中的shap段和gaus段用形函数库和高斯积分库作替换即得GES文件。
上述命令表示由sdispe.pde文件生成seuq4.ges和seuq4.for。
输入上述命令后回车，FEPG即可生成用户所需的单元子程序。
把新生成的单元子程序连接到单元计算的E程序中，使用FEPG的let命令，具体做法是在“sdispe”文件夹上单击鼠标右键弹出let命令，选择并单击该命令即可实现把单元子程序连接到单元计算的E程序中。
第七节边界单元
在上面讲解线弹性力学方程PDE文件时，我们得到了线弹性力学虚功方程的弱形式。
上面所说的PDE文件的给出了虚功方程弱形式中体积积分的描述，这里我们将考虑如何描述边界积分项。
FEPG是通过添加边界单元来描述边界积分的，边界单元的作用相当于加力的边界条件或弹簧边界，这就等效于偏微分方程的第二类或第三类边界。
下面我们看边界单元的PDE文件。
dispu,v,给出未知函数名u，v
coorx,给出单元局部坐标系下的坐标变量名
Shapl2给出单元形状类型符和节点个数，l表示线单元，
2表示两节点单元
gaus2给出每个方向上积分点的个数
$c6fx=prmt
（1）从前处理中取出切向表面力，
$c6表示插入FORTRAN源程序
$c6fy=prmt
（2）从前处理中取出法向表面力
Stif给出刚度矩阵
dist=[u;u]*0.0此算例为力的边界，所以这里刚度矩阵等于零
load=+[u]*fx+[v]*fy给出载荷向量，对应于
end结束符
同样使用pdeges命令由PDE文件生成相应的单元子程序。
假设上面的文件名为loadbnd.pde，单击FEPG界面主菜单中的“Run”，弹出FEPG命令对话框，在输入栏输入
pdegesloadbndseull2
按回车，系统即可由loadbnd.pde文件生成seull2.ges文件和边界单元的单元计算的FORTRAN源程序文件seull2.for。
把新生成的单元子程序连接到单元计算的E程序中，同样使用FEPG的let命令，具体做法是在“sdispe”文件夹上单击鼠标右键弹出let命令，选择并单击该命令即可实现把单元子程序连接到单元计算的E程序中。
下面我们再更进一步考虑有弹簧边界的情景，即有等效于偏微分方程的第三类边界条件的边界。
此时边界力矢量
和
和边界上的位移有关，假设
（9）
其中
，
分别为切向和方向的弹簧刚度。
将
和
的表达式代入（8）式得有第三类边界的弱虚功方程
（10）
不难看出（10）式比（8）多一项，此即由第三类边界引起的刚度矩阵的修正。
下面我们再给出有第三类边界的边界单元的PDE文件的填写
dispu,v,给出未知函数名u，v
coorx,给出单元局部坐标系下的坐标变量名
Shapl2给出单元形状类型符和节点个数，l表示线单元，
2表示两节点单元
gaus2给出每个方向上积分点的个数
$c6fx=prmt
（1）从前处理中取出切向表面力，
$c6表示插入FORTRAN源程序
$c6fy=prmt
（2）从前处理中取出法向表面力
$c6ax=prmt（3）从前处理中取出弹簧切向刚度
$c6ay=prmt（4）从前处理中取出弹簧法向刚度
Stif给出刚度矩阵
dist=[u;u]*ax+[v;v]*ay给出刚度矩阵的修正，对应于
load=+[u]*fx+[v]*fy给出载荷向量，对应于
end结束符
用FEPG的pdeges命令就可由此文件生成有弹簧边界的边界单元子程序，再用let命令将单元子程序连接到E程序。
这里还有一点需要说明的是需在前处理中多定义两个参数，具体做法是：
a.修改前处理中s.mti文件bf段，修改
升级会员 