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

第五章 幂律流体问题单元程序的理论文本 在FEPG3.0的模板库（APPWIZARD）中共有十种坐标系的单元程序，我们在此主要介绍三维直角坐标系下的单元程序。其它坐标系下的单元程序与此类同。对于幂律流体问题，在此采用欧拉法描述流场。 5．1 稳态问题 平衡方程 几何方程 本构方程 其中表示有效应力；表示有效应变率；为系数和指数。 边界条件 第一类边界条件 第二类边界条件 第三类边界条件 运用牛顿迭代法求速率 虚功方程如下： 将本构方程代入上式，可得 为了简化公式的推导，可将上面的方程简写为 将上式做一下变换 设得 由，得 此幂律流体问题采用体积不可压缩条件。 为此在虚功方程中增加两项。 定义泛函F如下 其中表示位移和压强。虚功方程等价于，由牛顿迭代法如下 得 上式两端同时增加一项，得 设，得 由泛函F的定义，可得 (1.2) 同样可得 于是可给出方程(1.1)右端的具体表达式 (1.3) 将上面的方程写入系统文件中（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=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) $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*exyn funa=+[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信息段写入分布式刚度矩阵如下 stif dist = +[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;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信息段写入载荷向量如下 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