幂率流体问题单元程序的理论文本

1、第五章 幂律流体问题单元程序的理论文本在FEPG3.0的模板库（APPWIZARD）中共有十种坐标系的单元程序，我们在此主要介绍三维直角坐标系下的单元程序。其它坐标系下的单元程序与此类同。对于幂律流体问题，在此采用欧拉法描述流场。51 稳态问题1 平衡方程 2 几何方程 3 本构方程 其中表示有效应力；表示有效应变率；为系数和指数。4 边界条件第一类边界条件第二类边界条件第三类边界条件5 运用牛顿迭代法求速率虚功方程如下：将本构方程代入上式，可得为了简化公式的推导，可将上面的方程简写为 将上式做一下变换 设得 由，得 此幂律流体问题采用体积不可压缩条件。 为此在虚功方程中增加两项。 定义泛函F

2、如下其中表示位移和压强。虚功方程等价于，由牛顿迭代法如下 得 上式两端同时增加一项，得 设，得 由泛函F的定义，可得 (1.2) 同样可得 于是可给出方程(1.1)右端的具体表达式 (1.3)6 将上面的方程写入系统文件中（PDE文件或GES文件）在FUNC信息段写入广义应变函数的定义如下func$cv exx=un/x$cv exy=un/y$cv exz=un/z$cv eyy=vn/y$cv eyx=vn/x$cv eyz=vn/z$cv ezz=wn/z$cv ezx=wn/x$cv ezy=wn/y$cv exn=exx$cv eyn=eyy$cv ezn=ezz$cv exyn=e

3、xy+eyx$cv exzn=exz+ezx$cv eyzn=eyz+ezy$cv emid=exn*exn+eyn*eyn+ezn*ezn$cv emid=emid+(exyn*exyn+exzn*exzn+eyzn*eyzn)/2$cv emid=emid*2./3+1.0e-10$cv fact=4./3*pk*pc*emid*(pc-2)$cv ef1=(pc-1)*2*exn$cv ef2=(pc-1)*2*eyn$cv ef3=(pc-1)*2*ezn$cv ef4=(pc-1)*2*eyzn$cv ef5=(pc-1)*2*exzn$cv ef6=(pc-1)*2*exynfun

4、a=+u/x funb=+v/y func=+w/z fund=+v/z+w/y fune=+w/x+u/z funf=+u/y+v/x fung=+u/x*exn+v/y*eyn+w/z*ezn+v/z*eyzn/2+w/y*eyzn/2+w/x*exzn/2+u/z*exzn/2+u/y*exyn/2+v/x*exyn/2 funh=+u/x+v/y+w/z其中un，vn，wn是上一个迭代步求得的速率，u，v，w是未知量。对照式(1.1) 和(1.2)在STIF信息段写入分布式刚度矩阵如下stifdist =+funa;funa*fact*(emid)+funb;funb*fact*(em

5、id)+func;func*fact*(emid)+fund;fund*fact*(emid/2)+fune;fune*fact*(emid/2)+funf;funf*fact*(emid/2)+fung;fung*fact*(pc-1)*4./3)+funh;p+p;funh-p/x;p/x*det*(2./3.)*1.0e-6/pk-p/y;p/y*det*(2./3.)*1.0e-6/pk-p/z;p/z*det*(2./3.)*1.0e-6/pk其中增加的后三行的作用是保证系数矩阵分解能无条件进行；且能保证解的收敛性。对照式(1.1) 和(1.3)在LOAD信息段写入载荷向量如下loa

6、d = +u*fu+v*fv+w*fw+funa*fact*ef1*(emid)+funb*fact*ef2*(emid)+func*fact*ef3*(emid)+fund*fact*ef4*(emid/2)+fune*fact*ef5*(emid/2)+funf*fact*ef6*(emid/2)方程右端项的边界力部分应该填写在边界单元的GES文件中，如下load=+u*fx+v*fy+w*fz7 注意事项如需在节点上施加集中力，应在前处理数据文件的DISP表格中填写。如需在边界面上施加分布力，应在此边界面上制作边界单元（包括前出理数据文件和单元GES文件）。模板中的边界单元均采用局部坐标

7、系，局部坐标系的产生办法如下：由单元的第一个节点到第二个节点的矢量作为局部坐标系的X轴，垂直于边界单元且与边界单元的节点顺序满足右手螺旋法则的矢量作为局部坐标系的Z轴，同时与X轴和Z轴垂直，且能与X轴和Z轴构成右手坐标系的矢量作为Y轴。因此边界单元的GES文件中的是指局部坐标系下的边界力。NFE文件具体实现了迭代过程，即首先由上一个迭代步的位移来形成刚度矩阵和载荷向量，然后求解得出当前迭代步的位移，如此反复，直至收敛。52 波动问题（时间离散采用Crank-Nicolson格式）1 平衡方程 2 几何方程 3 本构方程 其中表示有效应力；表示有效应变率；为系数和指数。4 边界条件 第一类边界条

8、件 第二类边界条件 第三类边界条件 5 运用牛顿迭代法求位移虚功方程如下： 将本构方程代入上式，可得 为了简化公式的推导，可将上面的方程简写为 将上式做一下变换 设得 由，得 此幂律流体问题采用体积不可压缩条件。 为此在虚功方程中增加两项。 定义泛函F如下其中表示位移和压强。虚功方程等价于，由牛顿迭代法如下 得 上式两端同时增加一项，得 设，得 由泛函F的定义，可得 (2.2) 同样可得 于是可给出方程(1.1)右端的具体表达式 (2.3)6 将上面的方程写入系统文件中（PDE文件或GES文件）在FUNC信息段写入广义应变函数的定义如下func$cv exx=un/x$cv exy=un/y$

9、cv exz=un/z$cv eyy=vn/y$cv eyx=vn/x$cv eyz=vn/z$cv ezz=wn/z$cv ezx=wn/x$cv ezy=wn/y$cv exn=exx$cv eyn=eyy$cv ezn=ezz$cv exyn=exy+eyx$cv exzn=exz+ezx$cv eyzn=eyz+ezy$cv emid=exn*exn+eyn*eyn+ezn*ezn$cv emid=emid+(exyn*exyn+exzn*exzn+eyzn*eyzn)/2$cv emid=emid*2./3+1.0e-10$cv fact=4./3*pk*pc*emid*(pc-2)

10、$cv ef1=(pc-1)*2*exn$cv ef2=(pc-1)*2*eyn$cv ef3=(pc-1)*2*ezn$cv ef4=(pc-1)*2*eyzn$cv ef5=(pc-1)*2*exzn$cv ef6=(pc-1)*2*exynfuna=+u/x funb=+v/y func=+w/z fund=+v/z+w/y fune=+w/x+u/z funf=+u/y+v/x fung=+u/x*exn+v/y*eyn+w/z*ezn+v/z*eyzn/2+w/y*eyzn/2+w/x*exzn/2+u/z*exzn/2+u/y*exyn/2+v/x*exyn/2 funh=+u/x

11、+v/y+w/z其中un，vn，wn是上一个迭代步求得的速率，u，v，w是未知量。 对照式(2.1)和式(2.2)在STIF信息段写入分布式刚度矩阵如下stifdist =+funa;funa*fact*(emid)+funb;funb*fact*(emid)+func;func*fact*(emid)+fund;fund*fact*(emid/2)+fune;fune*fact*(emid/2)+funf;funf*fact*(emid/2)+fung;fung*fact*(pc-1)*4./3)+funh;p+p;funh-p/x;p/x*det*(2./3.)*1.0e-6/pk-p/y

12、;p/y*det*(2./3.)*1.0e-6/pk-p/z;p/z*det*(2./3.)*1.0e-6/pk其中增加的后三行的作用是保证系数矩阵分解能无条件进行；且能保证解的收敛性。 对照式(2.3)在MASS信息段写入集中式质量矩阵如下 mass lump = rouu1 +rouu2 +rouu3 +rouu4 +rouu5 +rouu6 +rouu7 +rouu8 rouv1 +rouv2 +rouv3 +rouv4 +rouv5+rouv6+rouv7+rouv8rouw1+rouw2+rouw3+rouw4+rouw5+rouw6+rouw7+rouw8 其中 对照式(2.1)和

13、式(2.3)在LOAD信息段写入载荷向量如下load = +u*fu+v*fv+w*fw+funa*fact*ef1*(emid)+funb*fact*ef2*(emid)+func*fact*ef3*(emid)+fund*fact*ef4*(emid/2)+fune*fact*ef5*(emid/2)+funf*fact*ef6*(emid/2)式(2.3)中的边界力部分应该填写在边界单元的GES文件中，如下load=+u*fx+v*fy+w*fz7 填写NFE文件 式(2.1)进行空间离散化后，写成矩阵形式可得 其中M为质量矩阵，为加速度向量，为速度向量，S为刚度矩阵，F为载荷向量。 我

14、们采用Crank-Nicolson的时间离散格式 由 可将上式改写为 对照上式，在EQUATION信息段填写方程组总体刚度矩阵如下 matrix = S*DT/2+M 在EQUATION信息段填写方程组右端向量如下 FORC=F*DT+M*U1-S*U1*DT/2 其中DT表示，U1表示。一 注意事项 与时间有关的问题，如需在节点上施加集中力，应在BFT.FOR文件中填 写。 与时间有关的问题如需设定初值，应该填写两个DISP表，第二个表放初 值，第一个表无意义。 如需在边界面上施加分布力，应在此边界面上制作边界单元（包括前处理 数据文件和单元GES文件）。 模板中的边界单元均采用局部坐标系，

15、局部坐标系的产生办法如下：由单元的第一个节点到第二个节点的矢量作为局部坐标系的X轴，垂直于边界单元且与边界单元的节点顺序满足右手螺旋法则的矢量作为局部坐标系的Z轴，同时与X轴和Z轴垂直，且能与X轴和Z轴构成右手坐标系的矢量作为Y轴。因此边界单元的GES文件中的是指局部坐标系下的边界力。 NFE文件具体实现了迭代过程，即在每一时间步首先由上一个迭代步的 位移来形成刚度矩阵和载荷向量，然后求解得出当前迭代步的位移，如此反 复，直至收敛。 用户可以通过修改文本文件TIME0来确定计算的初始时刻、结束时刻和时间步长。53 已知速率求应力的单元程序对于幂律流体问题，稳态问题和波动问题的已知速率求应力的单

16、元程序是一样的，均从幂律本构方程出发，采用最小二乘法。1 本构方程 其中表示有效应力；表示有效应变率；为系数和指数。 应力所做的虚功 将上式变化为 设得 由，得 可以得出 2 最小二乘法由上式得 3 将上面的方程写入系统文件中（PDE文件或GES文件）在STIF信息段写入计算应力的程序及分布式刚度矩阵如下$cv exx=u/x$cv exy=u/y$cv exz=u/z$cv eyy=v/y$cv eyx=v/x$cv eyz=v/z$cv ezz=w/z$cv ezx=w/x$cv ezy=w/y$cv exn=exx$cv eyn=eyy$cv ezn=ezz$cv exyn=exy+ey

17、x$cv exzn=exz+ezx$cv eyzn=eyz+ezy$cv emid=exn*exn+eyn*eyn+ezn*ezn$cv emid=emid+(exyn*exyn+exzn*exzn+eyzn*eyzn)/2$cv emid=emid*2./3+1.0e-10$cv fact=4./3*pk*pc*emid*(pc-2)$cv funa=+u/x$cv funb=+v/y$cv func=+w/z$c6 fsa=funa*(emid)$c6 fsa=fsa*fact+p$c6 fsb=funb*(emid)$c6 fsb=fsb*fact+p$c6 fsc=func*(emid

18、)$c6 fsc=fsc*fact+p$cv fund=+v/z+w/y$cv fune=+w/x+u/z$cv funf=+u/y+v/x$c6 fsd=fund*(emid/2)$c6 fsd=fsd*fact$c6 fse=fune*(emid/2)$c6 fse=fse*fact$c6 fsf=funf*(emid/2)$c6 fsf=fsf*factdist = sa;sa*0.0在MASS信息段写入集中式质量矩阵如下masslump = 1.sa1+1.sa2+1.sa3+1.sa4+1.sa5+1.sa6+1.sa7+1.sa8 1.sb1+1.sb2+1.sb3+1.sb4+1.sb5+1.sb6+1.sb7+1.sb8 1.sc1+1.sc2+1.sc3+1.sc4+1.sc5+1.sc6+1.sc7+1.sc8 1.sd1+1.sd2+1.sd3+1.sd4+1.sd5+1.sd6+1.sd7+1.s