瞬态和非线性问题的有限元方法及其FEPG实现

1、瞬态和非线性问题的有限瞬态和非线性问题的有限元方法及其元方法及其FEPGFEPG实现实现中科院研究生院 2006年6月主要内容应用FEPG求解问题回顾线性椭圆问题的算法文件瞬态问题及其处理NFE文件的基本结构NFE文件的编写方法非线性问题及其处理PDE文件的COEF功能NFE文件的COEF功能FEPG的有限元文件与元件程序的关系一些调试命令应用FEPG求解问题回顾用户编写有限元文件lPDE 文件lFBC 文件(可选)lNFE 文件lGIO 文件lGCN 文件使用公式库生成程序PDE文件简单回顾DEFI 段落l变量名、函数名声明l坐标变量名称声明l参数名声明及其赋值l形函数和积分方法选择FUNC

2、 段落l用户定义函数的函数体STIF 段落l刚度矩阵计算表达式(微分方程弱形式)LOAD 段落l载荷向量的计算表达式MASS段落DAMP段落END线性椭圆问题的算法文件ldefilstif Slmass Mlload Fltype elmdty llstep 0lequationlmatrix = SlFORC=FlSOLUTION Ulwrite(s,unod) Ulend如何由刚度矩阵计算线性代数方程组矩阵如何计算方程组的右端项结果文件的读写一些声明稳态热传导问题的NFE文件defistif Smass Mload Ftype pmdty lstep 0equationmatrix = S

3、FORC=Fsolution UWRITE(S,unod) Uenddefidisp ucoor x yshap q 4gaus qmate ek eq rhoc 0.044 0.0 1.0d0stif dist=+u/x;u/x*ek+u/y;u/y*ekload=+u*eqend线性椭圆问题的算法文件DEFI 段落段落l给出变量的符号EQUATION 段落段落l给出线性代数方程组的结构SOLUTION 段落段落END 段落段落NFE文件的DEFI段落定义代数方程组左右端表达式中用到的一些定义代数方程组左右端表达式中用到的一些变量的含义，以及开数组大小的控制信息。变量的含义，以及开数组大小的

4、控制信息。STIF给出刚度矩阵所用符号MASS给出质量矩阵所用符号DAMP给出阻尼矩阵所用符号LOAD给出荷载向量所用符号TYPE给出方程类型 （椭圆型为E、抛物型为P、波动为W）MDTY给出质量矩阵和阻尼矩阵的类型 （集中矩阵为L、分布矩阵为D）STEP给出要保存刚度矩阵、质量矩阵的步数 （非负整数）瞬态问题抛物问题双曲问题l一阶l二阶时间导数的离散瞬态问题弱形式的描述抛物方程抛物方程波动方程波动方程瞬态问题弱解形式举例非定常热传导问题瞬态问题的PDE文件defidisp ucoor x yshap q 4gaus qmate ek ef rho 0.044 0.0 1.0d0lmass q

5、 rhostif dist=+u/x;u/x*ek+u/y;u/y*ekload=+u*efend瞬态问题弱解形式举例波动方程瞬态问题的PDE文件defidisp ucoor x yshap q 4gaus qmate ek ef rho ec 0.044 0.0 1.0d0 1.0d0lmass q rholdamp q ecstif dist=+u/x;u/x*ek+u/y;u/y*ekload=+u*efend有限元离散后的方程组(用自由度表示)抛物方程波动方程M: 质量矩阵 C: 阻尼矩阵 S: 刚度矩阵时间导数的离散(向后差分)U0 为上一时刻问题的解时间导数的离散(向后差分)瞬态问

6、题的NFE文件: 抛物问题可以给出如下形式的线性代数方程组 其左端矩阵A和右端项b的计算方法A=M+S*M t b=F*M t+M*U0瞬态问题的NFE文件(抛物问题)defistif Smass Mload Ftype pmdty lstep 0equationmatrix = S*dt+MFORC=F*dt+M*U0solution UWRITE(S,unod) Uend不完备！U0哪里来？NFE文件的变量定义和赋值语句（vect, read 语句）Equationvect u0read(s,unod) u0matrix = S*dt+MFORC=F*dt+M*U0瞬态问题的NFE文件(抛

7、物问题)defistif Smass Mload Ftype pmdty lstep 0equationvect u0read(s,unod) u0matrix = S*dt+MFORC=F*dt+M*U0solution UWRITE(S,unod) Uend有限元离散后的方程组(用自由度表示)抛物方程波动方程M: 质量矩阵 C: 阻尼矩阵 S: 刚度矩阵波动方程的时间导数离散令 V=Ut，那么， Vt=Utt波动方程的时间导数离散瞬态问题的NFE文件: 波动问题可以给出如下形式的线性代数方程组 其左端矩阵A和右端项b的计算方法A=b=瞬态问题的NFE文件: 波动问题defistif Sma

8、ss Mdamp Cload Ftype wmdty lstep 0equationVECT U0,V0READ(S,unod) U0 v0matrix = S*(DT/2)*2+C*DT/2+MFORC=F*DT*DT/2+M*U0+M*V0*DT +C*U0*DT/2- S*U0*DT*DT/4solution UWRITE(S,unod) Uend不完备！V0 哪里来？NFE文件的向量运算功能SOLUTION UVECT U,U0,V0READ (S,unod) U0,V0V1=U/DT*2-U0/DT*2-V0WRITE(S,UNOD) U V0END线性问题NFE文件小结给出如何由单

9、元刚度矩阵、质量矩阵、阻尼矩阵等计算代数方程组矩阵给出如何计算方程组的右端项给出后处理的一些计算读写文件给出VECTOR定义语句NFE文件DEFI信息段NFE文件主要是给出线性代数方程组左端矩阵和右端向量的表达式，以便形成最终求解的代数方程组，NFE文件一般由四段内容组成（DEFI、COEF、EQUATION、SOLUTION），每一信息段以空行表示结束。定义代数方程组左右端表达式中用到的一些变量的含义，以及开数组大小的控制信息。STIF给出刚度矩阵所用符号MASS给出质量矩阵所用符号DAMP给出阻尼矩阵所用符号LOAD给出荷载向量所用符号TYPE给出方程类型（椭圆型为E、抛物型为P、波动为W

10、）MDTY给出质量矩阵和阻尼矩阵的类型（集中矩阵为L、分布矩 阵为D）STEP给出要保存刚度矩阵、质量矩阵的数量（非负整数）DEFI信息段信息段NFE文件EQUATION信息段给出经过线性化后的线性代数方程组的左端矩阵和右端向量的表达式。MATRIX = 给出代数方程组左端矩阵表达式给出代数方程组左端矩阵表达式FORC = 给出代数方程组右端向量表达式给出代数方程组右端向量表达式表达式与普通矩阵、向量的运算一致，但是没有矩阵乘法运算。表达式中用到的符号要加以定义或者赋值，否则系统可能随意进行取值，导致计算错误。书写NFE文件建议先写以上两个表达式，然后再写其它辅助性的段落。对于刚度矩阵、质量矩

11、阵、阻尼矩阵、载荷向量要在DEFI信息段中定义，系统自动在单元计算结束后进行相应赋值。对于系统全局变量例如时间步长dt，不需要定义和用户赋值，直接使用。对于其它变量例如已知结果等可以用向量说明VECT、或者变量说明VAR加以定义，通过READ语句赋值。NFE文件EQUATION信息段对定义的变量赋值对定义的变量赋值vect 与与 var 定义语句定义语句 vect 包含了所有求解的未知量的个数 var 对应于一个求解的未知量 如：对应于pde中 disp u v w vect x 表示包含u v w三个未知量 var x 表示只含有其中的一个未知量 var x y z 表示含有u v w三个未

12、知量NFE文件SOLUTION信息段给出求解结果存放的向量名，对结果的运算处理和结果的磁盘存储。对于非线性问题，迭代收敛性判断以及控制迭代的信息也在此段给出。SOLUTION u 给出线性代数方程组求解结果存放的向量名VECT uWRITE(s,unod) u本段可以进行向量和变量的运算本段可以进行向量和变量的运算：向量名 = 向量表达式变量名 = 变量表达式也允许用户在也允许用户在VECT和和VAR说明语句之后的任何地方插入说明语句之后的任何地方插入fortran源程序源程序：$c6 fortran源程序 （空6列插入fortran源程序）$c0 fortran源程序 （不空列插入fortr

13、an源程序）NFE文件END信息段NFE文件以空行后跟一个END表示文件的结束。NFE文件中字符不区分大小写，FEPG所需要书所需要书写的文件均不区分字符大小写。写的文件均不区分字符大小写。END语句之前可以插入fortran语句.形式为：(空行)FORTRAN begin ( 空六列)dimension id(kdgof,knode)END非线性问题的种类材料非线性材料非线性l非线性本构关系l某些材料参数的非线性几何非线性几何非线性l大变形边界非线性边界非线性l接触和碰撞非线性问题的求解非线性热传导问题非线性热传导问题的线性化采用简单迭代法进行线性化，即采用简单迭代法进行线性化，即求解过程给

14、定初值u0求解求解求解判断书写线性化后的微分方程表达式PDE文件的COEF功能disp ucoor x yshap q 4coef ungaus qmate eq 0.0stif$c6 ek=1.44d0+1.236d-3*un-4.63d-7*un*undist=+u/x;u/x*ek+u/y;u/y*ekload=+u*eqend?NFE文件的COEF功能defistif smass mload ftype emdty lstep 0COEF UN EquationVECT UNREAD(S,UNOD) UNmatrix = sforc=f solution uwrite(s,unod)

15、u endNFE文件COEF信息段COEF信息段信息段给出单元计算所需的已知函数节点值。COEF 变量函数名 变量函数名迭代残差的计算绝对误差：相对误差：在NFE文件中计算残差solution uVect u un duread(s,unod) undu=u-un$c6 ERR=0.0d0$c6 UL=0.0d0%NOD%DOF ERR = ERR+du*2 UL = UL + U*2%DOF%NOD$C6 IF (ERR .LT. 1.0E-8 .OR. ERR .LT. 1.0E-8*UL) THEN$C6 OPEN(11,FILE=END,STATUS=NEW)$C6 CLOSE(11)

16、$C6 ENDIFwrite(s,unod) u endNFE文件总结 有限元方法最终求解的是离散后的线性代数方程组，用户必须给出代数方程组的左端矩阵和右端向量的表达式，然后得到要求解的线性代数方程组，这些表达式通过算法文件 *.NFE来 给出。对于大多数方程我们可以通过线性化方法得到其线性的表达式，在FEPG系统中已经把各种线性和非线性问题的常用算法放入了算法库中.NFE文件总结DEFI 段落COEF段落EQUATION段落SOLUTION段落例题四：孔周受力问题xyFFNFE库文件简介 NFE算法库中NFE文件命名规则如下： 椭圆型方程： fLu 抛物型方程： fLutuM波动方程： fL

17、utuCtuM22在前面的两个例子中所用的线性椭圆型算法ell.nfe、线性抛物线算法par.nfe、由位移求应力的最小二乘算法str.nfe都可以在系统库中相应找到，直接使用。如果系统算法库中没有适合的算法文件可用，那么就需要用户参照库中文件书写自己的NFE文件。 NFE库文件简介 ELL.NFE 求解线性椭圆型方程的算法程序。PAR.NFE 求解线性抛物型方程的算法程序，时间离散采用Crank-Nicolson格式。PARB.NFE 求解线性抛物型方程的算法程序，时间离散采用向后差分格式。WAVE.NFE 求解线性波动方程的算法程序，时间离散采用Wave速度法格式。NEWMARK.NFE

18、求解线性波动方程的算法程序，时间离散采用Newmark格式。NELL.NFE 求解非线性椭圆型方程的算法程序。NPAR.NFE 求解非线性抛物型方程的算法程序,时间离散采用Crank-Nicolson格式。NPARB.NFE 求解非线性抛物型方程的算法程序，时间离散采用向后差分格式。NWAVE.NFE 求解非线性波动方程的算法程序，时间离散采用Wave速度法格式。NNW.NFE 求解非线性波动方程的算法程序，时间离散采用Newmark格式。STR.NFE 已知位移场求应力场的算法程序，采用最小二乘法。ell椭圆，ellipticpar抛物, parabolicwave波动（采用只求位移和速度的

19、一种算法）newmark波动（采用newmark算法，同时求解位移、速度和加速度）以上算法前面加符号n表示非线性。程序结构程序结构FEPG采用元件化技术来实现这五个步骤的自动生成，我们书写的NFE、PDE（VDE、FBC）、GCN、GIO文件中，nPDE（VDE、FBC）对应生成单元计算程序；单元计算程序；nNFE对应生成组集线性代数方程组的程序对应生成组集线性代数方程组的程序，在该程序中调用单元计，在该程序中调用单元计算程算程序，另外序，另外NFE中的中的nSOLUTION段生成结果后处理程序段生成结果后处理程序这两部分程序随问题的变化而变化，是有限元计算解决问题的核心部这两部分程序随问题的

20、变化而变化，是有限元计算解决问题的核心部分。分。其余的程序初始化，边值计算，线性代数方程组求解这三部分程序功程序初始化，边值计算，线性代数方程组求解这三部分程序功能固定，在所有有限元计算中比较通用，对不同问题几无变化，因此能固定，在所有有限元计算中比较通用，对不同问题几无变化，因此在在FEPG系统中按通用方式给出，不需要用户书写文件，在GCN文件中直接调用即可。程序结构程序结构FEPG解决用户问题的主体程序流程解决用户问题的主体程序流程（数据准备部分略）程序间关系程序间关系nSTART、BFT、SOLV程序是根据GCN文件，从系统库中提取相应的start*.epg、start*.gpg、bft.epg、bft.gpg、solv.epg、solv.gpg文件，由EPG和GPG再生成.for文件，编译连接成.obj、.exe文件。这三个文件的生成相对简单，用户可以通过修改.epg、.gpg文件来修改这几个程序。nNFE文件生成E和U程序的e%*.epg、e%*.gpg、u%*.epg、u%*.gpg文件