当前位置:文档之家› 有限元教材-第十章 有限元程序设计

有限元教材-第十章 有限元程序设计

有限元教材-第十章 有限元程序设计
有限元教材-第十章 有限元程序设计

第十章有限元程序设计

有限元方法作为一门系统的技术,仅学会了它的基本理论是远远不够的,只有形成完整的计算程序,问题才最终得到了解决。完成这样的有限元程序设计是一项工作量很大的工程。本章就是要结合简单的有限元教学程序FEMED,简要介绍有限元程序设计技术。FEMED是专为有限元程序设计教学编制的程序,它不包含复杂的前后处理功能,可进行平面问题及平面桁架的线弹性静力分析,在程序结构上与大型程序类似,具有计算单元的任意扩充功能,在方程的组集和求解上也采用了较为流行的变带宽存储方式。

有限元程序大致可分为两类,第一类是专用程序,主要用于研究或教学,一般这类程序规模较小,前后处理功能较弱。用于研究的程序能够解一些特殊的问题,满足研究工作的需要。而教学程序则是为了学生了解有限元的主要结构和设计方法设计的,程序比较简单,FEMED就属于这类程序。第二类是大型通用程序,是大型结构分析的得力工具,目前国际上流行的大约有2000多种。常用的有NASTRAN、MARC、ANSYS、ADINA和ABAQUS等。这类程序一般前后处理功能比较强,有友好的界面,能进行大型计算,但往往无法完成具有特殊要求的计算。通过本章的学习,使读者初步掌握有限元编程的基本方法,具有开发特殊功能的专用程序或为通用程序开发具有特殊功能的计算模块的能力。

§10.1有限元程序的基本结构

有限元程序一般包括三项基本内容:前处理、结构分析和后处理。早期有限元分析软件的研究重点在于推导新的高效率求解方法和高精度的单元,随着数值分析方法的逐步完善,尤其是计算机内存和运算速度的飞速发展,整个计算系统用于求解运算的时间越来越少,加之求解问题的日益大型化和复杂化,使得数据准备和运算结果的表现问题日益突出。因此目前几乎所有的商业化有限元程序系统都有功能很强的前后处理模块,这直接关系到分析软件的可推广性。它是商用有限元软件不可或缺的部分,但它不是有限元的中心部分,在本书中不作详细介绍。

10.1.1 前处理

在结构分析之前,必须完成的工作就是前处理的任务,它主要包括以下内容:结构造型、选取单元类型,、网格剖分,从而形成节点数、节点编码、节点坐标、单元数、单元节点编码等。另外根据所解题目的不同,还需读入不同的材料参数、边界约束条件及载荷工况等。有限元前处理程序的功能就是为用户提供一种工具,使其能尽可能方便地完成上述工作。而这一问题解决的好坏直接关系到程序使用者需付出的工作量,而网格剖分质量还直接影响计算精度。很多程序都建立了对用户非常友好的图形用户界面,使用户能以可视图形方式直观快速地进行结构造型和网格自动划分,生成有限元分析所需数据。

在上述工作中工作量最大的就是结构造型和网格剖分部分,目前越来越倾向于由专用软件完成上述工作,所以当今所有的商业化有限元软件开发商都开发了和著名的CAD软件的接口,从而极大地提高了设计水平和效率。今天,工程师可以在集成的CAD和FEA软件环境中快捷地解决一个在以前无法应付的复杂工程分析问题。

在FEMED程序中,只包含了很简单的节点、网格及给定位移的自动生成。所有前处理功能是

由INPUT及相关子程序完成的。

10.1.2 结构分析

结构分析部分是有限元程序的主体部分,主要需计算插值函数矩阵[N],插值函数的导数矩阵[B],进而进行数值积分得到单元刚度矩阵,组集总刚度矩阵及总载荷向量,解线性代数方程组,得到节点位移。对于动力问题则需计算质量矩阵,求解特征值问题。而对于动力响应问题及非线性问题,则需进行多增量步的计算。在得到节点位移的基础上还需进行应力的计算。

在本章中重点介绍的就是结构分析中的静力分析部分,在FEMED程序中主要包含的也就是这部分。通过这部分的学习,希望读者能掌握有限元程序的基本编程技巧。

10.1.3 后处理

在大型分析软件中,程序的后处理功能也是非常重要的,在有限元结构分析完成后,通过强有力的后处理图形功能,可以给出各应力分量、等效应力等的分布云图,结构的变形图,振型图和实

图10.1

时动画等,使得使用者对各物理量在整个结构中的变化情况有一个全面的认识,也可检查网格剖分和所加载荷及约束是否有误。否则面对输出成千上万个计算结果数据,往往使人如坠雾中,不得要领,需花费大量的时间进行数据的整理和解释。

由于是简单的教学程序,在FEMED程序中不包含后处理功能,仅仅是将计算结果输出。

图10.1给出了一般求解静力学问题的程序框图,FEMED程序的结构基本如图所示,但没有包含带宽优化部分。

§10.2 数组的半动态分配

在有限元程序中需开辟许多不同功能的数组,而这些数组的大小则与所解题目的大小及类型有关,如节点坐标数组、节点位移数组等与节点数相关,而单元节点编号数组等则与单元数相关。此外,上述数组的大小还和自由度维数等相关。对于一个有限元程序而言,应具有求解各类问题及各

种规模问题的机动性,因此求解不同的问题,应定义不同大小的数组。FORTRAN 语言规定数组的大小必须预先给定,在子程序中可以定义大小可变化的可调数组,但其内存分配必须由该子程序的上一级模块确定。

最简单的办法就是每个数组都定义得足够大,使得要求解的问题都能被满足,但这显然会带来内存的巨大浪费,因此,一般采用数组半动态分配的方法来解决上述问题,FEMED 也采用了这一方法。

在主控程序PCONTR 中开辟了两个大数组A(*)和M(*),分别存储实型量和整型量,将各数组中的数顺序存放在大数组中,在读入控制信息NJ (节点数)、NE (单元数)、NUMMAT (材料数)、NDF (节点自由度数)、NDM (空间维数)、NEN (单元节点数)等之后,就可进行存储单元的分配。在图10.2中以大数组A 的存储为例,先存放节点坐标,占用NJ ×NDM 个存贮单元,再存放材料参数,占用NUMMA T ×20个存贮单元……

图10.2

在上述过程中,最重要的是计算出各数组存放的起始位置又称作地址,如NXC 为XC 数组在A 中的地址,ND 为D 数组在A 中的地址……,这项工作由POINT 子程序完成。需要分别计算出大数组A 和M 中所存数组的所有地址,最后计算出A 和M 数组需要的总长度,与预设的数组长度作比较,若没有超界计算继续进行,若超界则计算停止,显示错误信息。

需要指出的是,在上述数组的半动态分配时,将一系列存储容量不定的一维、二维数组均用地址对应关系分配到事先容量固定的大数组A 和M 中去,对于大数组A 和M 来说存储量是不变的,但各数组在大数组中的存储分配却是动态的。我们把这种关系称为数组的半动态分配,这一过程是通过哑实结合完成的。这些数组称为可调数组,即其数组长度是可调的。在子程序的数组说明中,一维数组写成下列形式是等价的,对数组长度没有限制:

F(*),F(1),F(NJ),

对于二维数组,则只可一个方向自由变动。由于FORTRAN 规定二维数组先变行后变列,因此行的长度必须是确定的。,在子程序的数组说明中,下列形式是等价的:

XC(2,*),XC(2,1),XC(2,NJ),XC(NDM,*),(NDM=2)

注意,只有对上一模块中已定义过的数组才可作如上的数组说明,否则,则必须准确确定数组的长度。

§10.3 构造单元刚度矩阵及总刚度矩阵的组集

构造单元刚度矩阵是由单元模块完成的,不同的单元调用不同的模块,在FEMED 中提供了平面桁架单元(TRUSS )和平面单元(PLANE )两个单元计算模块。读者可根据需要随意开发单元类型,通过ELEMLIB 子程序接入程序中。针对不同的单元,在不同的单元模块中构造不同的[B]、[D]矩阵,单元刚度矩阵可表示为:

V

dV =?e e T K B DB

(10.1)

对于常应变单元,上述被积函数为常数,直接用被积函数值乘积分域的大小即可。而对于等参元等非常应变单元,则需对上式做GAUSS 积分,最终得到一个二维的单元刚度矩阵。

将单元刚度矩阵正确组集为总刚度矩阵,这一工作是由ADDSTF 子程序完成的。下面将举例说明总刚度矩阵的组集方法。

例1:平面桁架结构如图示,共有5个节点,6个单元。每个单元形成一个4×4的单元刚度矩阵,总刚度阵为10×10。

图10.3

以1号单元为例,其左节点为2,右节点为4,其4×4的单元刚度矩阵的各个参数应迭加在(10.2)式总刚度矩阵中的图示位置上。需要指出的是,总刚度矩阵的一个位置上往往有多个单元刚度矩阵做贡献,如与3号节点对应的(5,5),(6,5),(6,6)位置上就有2、3、4和6号单元都作了贡献。

全部单元组集完毕后,总刚度矩阵如式(10.2),形成一对称带状阵,即仅对角元附近的元素不为零,而远离对角线的元素皆为零。解题越大带状性质越明显。

???

??

?

???

??

?

??

?

???????

????????????????????????

?????????????????

=0000000044

4342410033323100222111][称

K

(10.2)

若储存的是总刚度矩阵的上三角元素,则带宽为每行从对角元,到最后一个非零元素的个数;若储存的是总刚度矩阵的下三角元素,则从第一个非零元素到对角元元素的个数称为这行的带宽。每行的带宽各不相同,且与节点编号的方式有关,在大型商业软件中,一般都包含带宽优化模块,就是为了减少总体带宽。

§10.4 总刚度矩阵的变带宽存储

总刚度阵占用的空间是有限元占用空间的绝大部分,一个软件的解题规模往往是受此限制。所以,一直以来人们在这方面想了很多的办法。很显然,将总刚度阵作为一个二维数组,完全储存起来是最直接和最自然的方法,但这也是占用空间最大、最不合理的方法,一般很少采用。

考虑到总刚度阵具有带状稀疏对称的特征,使用半带宽存储无疑是一个最自然的改进,它利用一个二维数组,只存储每行从对角元到这一行的最后一个非零元,如图10.4所示

??

?

????

??

??

?

?

??

?

?????

???

???

????????

?

?

???

?????

??????

??????

???

???

??????

?

???????

?

????????????????????????0 000000000000000000000000000000000000称对

图10.4

由10.4的示意图知,原10×10的矩阵现在可用10×4的矩阵存储,存储空间减少了很多,但若其中有个别行的带宽特别大,则这种方法依然有存贮空间的大量浪费。为解决上述问题,一般多采用变带宽存储又称一维压缩存储的方法,这样刚度阵中大量的零既不需要存储也不参与计算,大大减少了存储空间和计算量,是目前使用较为广泛的方法。

在有限元的发展历史中波前法也是一种不可不提的计算方法,当计算机的内存非常有限时,即使采用变带宽存储也无法求解较大规模的问题,波前法曾经起到了不可忽视的作用。由于波前法无需组集总刚,所以对内存要求非常低,曾经被广泛地使用,但它也有非常明显的缺点,内外存交换非常频繁,导致计算速度受到很大的限制。随着计算机技术突飞猛进地发展,计算机的内存有了本质的提高,所以波前法已经很少被采用了。本书重点介绍变带宽存储的方法。

变带宽存储就是考虑到总刚度矩阵的对称性及带状稀疏的特点,用一维数组将总刚度矩阵下三角(或上三角)中的非零元素储存起来,如图10.5中虚线以上部分,注意一行中夹在中间的零必须储存。

图10.5

设A(i,j)为具有带状对称性质的二维数组,若将其存储在一维数组B(k)中,则必须开辟一个存放对角元地址的数组JP(i),指明每一个对角元在一维数组中的位置。如上图中矩阵的JP(i)应为:1,3,6,9,11,15,18,21,24,27。若第i 行,第j 列的数据存在一维数组的第k 个位置上,则k 与i 、j 之间的对应关系为

k=JP(i)-i+j

(10.3)

第i 行的第一个非零元素在m i 列,则

m i =i-JP(i)+JP(i-1)+1 (i ≠1)

(10.4)

???????????????

?

???????????????????

???????????????????????

0000000000000000000000称

当i=1时,很明显m 1=1。有了(10.3、4)两式,我们完全能够对一维数组很方便地进行方程求解,对正定对称方程A(i,j)的求解方法都可以使用,只是对A(i,j)的操作都需先计算k ,然后对B(k)做同样的操作即可。另外,由于m i 为第i 行的第一个非零元素,所以对第i 行的操作只需从第m i 列开始即可。

对于上述10 10的矩阵,我们需要存储27个元素,节省有限,而对于大规模的问题则节省的内存和计算工作量都将是非常显著的。

另外需要指出的是,带宽的大小与编码方式有关,我们将中间有单元相连的节点称为相关节点,尽量减小相关节点间的编号差,则可减小带宽。减小带宽对于减小内存占用和减小计算工作量都是很有意义的。在图10.6的二维结构中,(a )、(b )两种编码方式的半带宽就相差很多,(a )中的最大半带宽为8,而(b)中的最大半带宽则为16。这会使内存占用和计算工作量都相差很多。

例2:试分别给出图10.6中两种编码的带宽、对角元地址及总刚长度。

解:j 节点对应的是2j-1行和2j 行,两节点若包含在同一单元中则称相关节点,设与j 节点相关的节点中,编号最小的为i ,则

2j-1行的带宽d=2(j-i)+1。2j 行的带宽d=2(j-i)+2 对于(

对于(

(a )编码形式的总纲长度为140,(b )编码形式的总纲长度为220。

上述带宽及对角元地址的计算在PROFIL 子程序中完成,所不同的是PROFIL 子程序中还自动将约束自由度从总刚度矩阵中剔除,而本例中没考虑约束情况。

§10.5 给定位移约束的处理

经过总刚度矩阵的组集,并将给定载荷加在载荷向量的相应位置,我们即可得到线性代数方程

(a) (b)

图10.6

Ku F

(10.5)

每行对应一个自由度,需要指出的是,与每个自由度对应的有两个量——节点位移u i 和节点载荷f i ,在这两个量中有且仅有一个是已知的。在例1中前四个自由度已知节点位移而节点力未知,后面各自由度则已知节点力而节点位移未知。在求解线性代数方程组时,要求右端项全部已知,而未知量应在左端的u 中,因此需要根据边界约束条件对方程组进行改造。

若第i 个自由度上u i =?i 为已知,而f i 未知,则有三种处理方法:

105.1 加大数法

加大数法是程序处理最简单的方法,但也是一种近似的方法。作法是令对角元k ii =M i ,f i =?i M i ,其中M i 远远大于第i 行中其他各量,则其他变量的影响就可忽略,因此有δi ≈?i 。需要指出的是M i 的选择是技术性比较强的,若M i 选择得不够大则其他变量的影响不能较好地消除,则误差较大,若M i 选择得太大则影响方程的性态,引起较大误差甚至奇异,因此这是一种简单但不能尽如人意的方法。

10.5.2 划行划列法

划行划列法是约束处理的一种精确的方法,由于节点力f i 未知,第i 个方程引入了一个新的未知数,所以对求解a 是无意义的,因此放弃这一方程,换为u i =?i 。作法是将第i 行除对角元外的各元素都充零,对角元充1,令f i = ?i 。即

i i ij ii f i j n j k k ?=≠===)

,1(0

1

完成上述操作就可实现给定约束的处理,处理之后为

11

1211112122222212001

0i n i

n i i n n ni

nn n n k k k k u f k k k k u f u f k k k k u f ????????????

????????????=??????????

??

?????????????

?????

(10.6)

很明显上述操作破坏了方程的对称性,是不能接受的,因此需将第i 列也都充零(除对角元),以恢复其对称性。但这一操作不能影响其他各方程,因此需将?i 对其它各行的影响移至等号右边,即做如下操作

11

121111121222222212

0000100n i i n i i i i i

n n nn n n n ni i k k k u f f k k k k u f f k u f k k k u f f k ????-????????????????

-????????

??????

??==????????????????

????????????????-?

???????

(10.7)

完成上述操作则完成了约束处理,既保持了方程的对称性,又精确引入了位移边界条件,若还有其他约束则再做同样处理即可。

10.5.3 消自由度法

划行划列的方法无疑是一个切实可行的方法,但若受约束自由度很多则方程中多储存了许多的零,而且在解方程中也增加了许多计算量,因此就有了对上述算法的改进算法。FEMED 就采用了

这一算法。

对于需要划掉的行列一开始就不再组集入总刚度阵,若有l 个自由度被约束,则方程组就从n 元线性代数方程组降为n-l 元,存储空间和计算量都会有比较大的降低。具体作法是,开辟一个数组ID(2,*),记录约束情况,即活化自由度与方程组的对应关系。对于例1,形成的ID 数组应为:

-1,-2;-3,-4;1,2;3,4;5,6;7,8;9,10。

负值代表该自由度被约束,如-3代表2节点的x 方向是第3个位移约束,而6则代表第5个节点的y 方向的自由度是与方程组的第6个未知数对应。在组集总刚度矩阵时,应先查ID 数组,例如原来应组集在k 85中的数现在则应该组集在k 41中,原来应组集在k 41中的数则无需组集,原来应组集在k 72中的数k ,对应的是第3自由度和第2个约束,现在则在总刚度矩阵中不体现,而应将f 3修改为f 3-k ?2。

在FEMED 中,ID 数组是在BOUN 和PROFIL 中自动形成的。

§10.6 单元应力及节点力的计算

在所有节点位移求出后,要进一步求出应力及反力等,这时需再次对单元循环,对各单元的积分点计算应力

e =σDBa

其中,a e 为单元的节点位移向量,B 、D 矩阵与前面各章相同分别为应变矩阵和弹性矩阵,在B 矩阵中代入不同的局部坐标(ξ, η)可以得到不同点的应力。但一般认为,积分点的应力精度最好,而节点应力是最差的,所以一般有限元程序得到的都是积分点应力而非节点应力,当需要节点应力时,往往再由积分点应力通过插值格式得到节点应力。

计算单元节点力可由下式得到

e e e T T e V

V

dV dV ==??σF B B DB a

将所有各单元上计算出的节点力迭加起来,即可得到各节点的节点力。对于给定位移的自由度,如此求得的即是节点反力,而对于自由节点则会求出一个很小的节点力,这是计算误差使然。在线性分析中,上述计算主要是为了得到节点反力,而在非线性分析中,上述节点力则作为残差参与迭代计算。

§10.7 FEMED 程序

FEMED 程序是一个可用来求解线弹性静力学问题的有限元教学程序,它不包含复杂的前后处理,目的是让读者通过程序的学习初步掌握有限元主体程序的初步编程技巧,能够读懂和编制有限元程序,FEMED 程序采用动态数组存储技术。在线性代数方程组的存储上,采用了变带宽存储技术。在约束节点的处理上,采取了预先消除自由度的方法。所有这些都使其内存能够得到最有效的利用。本程序由于篇幅所限只包含了平面单元和平面桁架单元两种单元,但它采用了单元库的设计方式,使得使用者可以方便地添加自己所需要的单元,丰富自己的单元库。为了方便读者阅读,在10.9节给出了本程序的子程序索引、输入索引及主要变量索引。在程序之后给出了3个简单的算例。

PROGRAM FEMED !MAIN 1 C.....FINITE ELEMENT PROGRAM FOR EDUCATIONAL PURPOSE ONLY. !MAIN 2 C.....AVAILABLE FOR SOLVING MULTI-TYPES ELEMENTS STATIC PROBLEMS !MAIN 3 COMMON/PSIZE/MAXM,MAXA !MAIN 4

OPEN(5,FILE='EXP.DAT') !MAIN 5 OPEN(6,FILE='EXP.OUT') !MAIN 6 MAXM=6000 !MAIN 7 MAXA=90000 !MAIN 8 CALL PCONTR !MAIN 9 STOP !MAIN 10 END !MAIN 11 C............................................

SUBROUTINE PCONTR !PCON 1 C.....MAIN CONTROLING SUBROUTINE !PCON 2

COMMON/M/M(6000) !PCON 3 COMMON/A/A(90000) !PCON 4 COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PCON 5 COMMON/PDATA/MIX,MID,MLD,MJP,NXC,ND,NXL,NUL,NS,NP,NF,NU,NSK !PCON 6 COMMON/PSIZE/MAXM,MAXA !PCON 7 COMMON/MDATA/NUMMAT !PCON 8 C.....READ IN CONTROL DATA !PCON 9

READ(5,*) NJ,NE,NUMMAT,NDF,NDM,NEN !PCON 10 WRITE(6,2000) NJ,NE,NUMMAT,NDF,NDM,NEN !PCON 11 C.....CALCULATE THE ADDRESS OF ARRAY !PCON 12

CALL POINT(NNEQ) !PCON 13 C.....READ MOST OF THE PROBLEM DATA !PCON 14

CALL INPUT(A(NXC),M(MIX),M(MID),A(ND),A(NF)) !PCON 15 C.....DEFINE THE PROFIL FOR GLOBAL STIFFNESS MATRIX STORAGE !PCON 16

CALL PROFIL(M(MJP),M(MID),M(MIX),NK) !PCON 17 NEND=NSK+NK !PCON 18 IF(NEND.GT.MAXA) STOP '{A} ERR' !PCON 19 CALL PZERO(A(NSK),NK) !PCON 20 C.....FORM LOAD ARRAY !PCON 21

CALL PLOAD (M(MID),A(NF),A(NU),NNEQ) !PCON 22 C.....FORM ELEMENT STIFFNESS !PCON 23

CALL PFORM(A(NF),M(MLD),A(NUL),A(NXL),A(NS),A(NP),A(NU),A(NXC), !PCON 24

1 M(MIX),A(ND),M(MID),M(MJP),A(NSK),2) !PCON 25

C.....LDLT TRIANGULAR DECOMPOSITION !PCON 26

CALL LDLT(A(NSK),M(MJP),NEQ) !PCON 27 C.....FORWARD ELIMINATION AND BACKWARD SUBSTITUTION !PCON 28

CALL FORBACK(A(NSK),A(NU),M(MJP),NEQ) !PCON 29 C.....PRINT DISPLACEMENTS !PCON 30

CALL PRTDIS(M(MID),A(NU),A(NF)) !PCON 31 CALL PZERO(A(NF),NNEQ) !PCON 32 https://www.doczj.com/doc/876813557.html,PUTE ELEMENT STRESSES AND THE REACTION !PCON 33

CALL PFORM(A(NF),M(MLD),A(NUL),A(NXL),A(NS),A(NP),A(NU),A(NXC), !PCON 34

1 M(MIX),A(ND),M(MID),M(MJP),A(NSK),3) !PCON 35

C.....PRINT THE REACTION !PCON 36

CALL PRTREAC(A(NF)) !PCON 37 2000 FORMAT(' NUMBER OF NODAL POINTS =',I5/ !PCON 38

1 ' NUMBER OF ELEMENTS =',I5/ !PCON 39

1 ' NUMBER OF MATERIALS =',I5/ !PCON 40

1 ' NUMBER OF DEGREES/NODE =',I5/ !PCON 41

1 ' NUMBER OF DIMENSION SPACE =',I5/ !PCON 42

1 ' NUMBER OF NODES/ELEMENT =',I5/) !PCON 43

RETURN !PCON 44

END !PCON 45 C............................................

SUBROUTINE POINT(NNEQ) !POIN 1 C.....CALCULATE THE ADDRESS OF ARRAY !POIN 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !POIN 3 COMMON/PSIZE/MAXM,MAXA !POIN 4 COMMON/PDATA/MIX,MID,MLD,MJP,NXC,ND,NXL,NUL,NS,NP,NF,NU,NSK !POIN 5 COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL,NST !POIN 6 COMMON/MDATA/NUMMAT !POIN 7 COMMON/A/A(90000) !POIN 8 NEN1=NEN+1 !POIN 9 NNEQ=NJ*NDF !POIN 10 NST=NDF*NEN !POIN 11 C.....CALCULATE THE ADDRESS OF INTEGRAL ARRAY {M} !POIN 12

MIX=1 !POIN 13 MID=MIX+NEN1*NE !POIN 14 MLD=MID+NNEQ !POIN 15 MJP=MLD+NST !POIN 16 MEND=MJP+NNEQ !POIN 17 IF(MEND.GT.MAXM) STOP '{M} ERR' !POIN 18 C.....CALCULATE THE ADDRESS OF REAL ARRAY {A} !POIN 19

NXC=1 !POIN 20 ND =NXC+NJ*NDM !POIN 21 NXL=ND +20*NUMMAT !POIN 22 NUL=NXL+NEN*NDM !POIN 23 NS =NUL+NST*2 !POIN 24 NP =NS +NST*NST !POIN 25 NF =NP +NST !POIN 26 NU =NF +NNEQ !POIN 27 NSK=NU +NNEQ !POIN 28 CALL PZERO(A(1),NSK) !POIN 29 RETURN !POIN 30 END !POIN 31 C............................................

SUBROUTINE INPUT(XC,IX,ID,D,F) !INPU 1 C.....READ MOST OF THE PROBLEM DATA !INPU 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !INPU 3 COMMON/PDATA/MIX,MID,MLD,MJP,NXC,ND,NXL,NUL,NS,NP,NF,NU,NSK !INPU 4 COMMON/PSIZE/MAXM,MAXA !INPU 5 COMMON/MDATA/NUMMAT !INPU 6 COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL,NST !INPU 7 DIMENSION XC(NDM,*),IX(NEN1,*),ID(NDF,*),D(20,*),F(NDM,*) !INPU 8 CHARACTER COOR*6,XX*18 !INPU 9 DATA XX/'NODAL COORDINATES'/COOR/'COOR'/ !INPU 10 CHARACTER FD*18,FORC*6 !INPU 11 DATA FORC/' FORC '/FD/' NODAL FORCE/DISPL'/ !INPU 12 DIMENSION P(1),S(1) !INPU 13 C.....INPUT NODAL COORDINATES DATA !INPU 14

CALL GENVEC(XC,XX,COOR) !INPU 15 C.....INPUT ELEMENT DATA !INPU 16

CALL ELDAT(IX) !INPU 17 C.....INPUT MATERIAL DATA !INPU 18

DO N=1,NUMMAT !INPU 19

READ(5,*) MA,IEL !INPU 20

WRITE(6,2000) MA,IEL !INPU 21

IF(MA.EQ.0) EXIT !INPU 22

D(20,MA)=IEL !INPU 23

CALL ELEMLIB(D(1,MA),P,XC,IX,S,P,NDF,NDM,1) !INPU 24 END DO !INPU 25 C.....INPUT BOUNDARY CONDITIONS !INPU 26

CALL BOUN(ID) !INPU 27 C.....INPUT FORCES AND DISPLACEMENTS !INPU 28

CALL GENVEC(F,FD,FORC) !INPU 29 2000 FORMAT(/' MATERIAL TYPE ',I4,' ELEMENT TYPE ',I4/) !INPU 30 RETURN !INPU 31 END !INPU 32 C.................................................

SUBROUTINE GENVEC(X,CD,FC) !GENV 1 C.....GENERATE REAL DATA ARRAYS BY LINEAR INTERPOLATION !GENV 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !GENV 3 CHARACTER CD*18,FC*6 !GENV 4 DIMENSION X(NDM,*),XL(6) !GENV 5 DO K=1,NJ+1 !GENV 6

READ(5,*) N,NG,(XL(I),I=1,NDM) !GENV 7

IF(N.EQ.0) EXIT !GENV 8

DO I=1,NDM !GENV 9

X(I,N)=XL(I) !GENV 10 END DO !GENV 11 IF(NG.NE.0) THEN !GENV 12

LI=(N-L)/NG !GENV 13

IF(LI.LT.1) STOP 'GENVEC ERR' !GENV 14

DO I=1,NDM !GENV 15

XL(I)=(X(I,N)-X(I,L))/LI !GENV 16 END DO !GENV 17

L1=L+NG !GENV 18

L2=N-NG !GENV 19 DO L=L1,L2 !GENV 20 DO I=1,NDM !GENV 21

LMLG=L-NG !GENV 22

X(I,L)=X(I,LMLG)+XL(I) !GENV 23 END DO !GENV 24 END DO !GENV 25 ENDIF !GENV 26 L=N !GENV 27 END DO !GENV 28 CALL OUTPUT(X,CD,FC) !GENV 29 RETURN !GENV 30 END !GENV 31 C.................................................

SUBROUTINE OUTPUT(X,CD,FC) !OUTP 1 COMMON /CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !OUTP 2 CHARACTER CD*25,FC*6 !OUTP 3 DIMENSION X(NDM,*) !OUTP 4 WRITE(6,2000) CD,(K,FC,K=1,NDM) !OUTP 5

DO J=1,NJ !OUTP 6

WRITE(6,3000)J,(X(L,J),L=1,NDM) !OUTP 7 END DO !OUTP 8 2000 FORMAT(/5X,A18/2X,5HNODE ,6(I6,A6)) !OUTP 9 3000 FORMAT(I6,6F12.4) !OUTP 10 RETURN !OUTP 11 END !OUTP 12 C...............................................................

SUBROUTINE ELDAT(IX) !ELDA 1 C.....ELEMENT DATA INPUT !ELDA 2

COMMON /CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !ELDA 3 DIMENSION IX(NEN1,*),IXL(9) !ELDA 4 WRITE(6,2000) (K,K=1,NEN) !ELDA 5 DO J=1,NE+1 !ELDA 6

READ(5,*) N,NN,NG,LK,(IXL(K),K=1,NEN) !ELDA 7

IF(N.EQ.0) EXIT !ELDA 8

N1=N+NN-1 !ELDA 9

IX(NEN1,N) =LK !ELDA 10 DO K=1,NEN !ELDA 11

IX(K,N) =IXL(K) !ELDA 12 END DO !ELDA 13 IF(NN.EQ.1) CYCLE !ELDA 14 DO I=N+1,N1 !ELDA 15

IX(NEN1,I) =LK !ELDA 16

I1=I-1 !ELDA 17

DO K=1,NEN !ELDA 18

IX(K,I)=IX(K,I1)+NG !ELDA 19

IF(IX(K,I1).EQ.0) IX(K,I)=0 !ELDA 20

END DO !ELDA 21 END DO !ELDA 22 END DO !ELDA 23 DO I=1,NE !ELDA 24 WRITE(6,3000) I,IX(NEN1,I),(IX(K,I),K=1,NEN) !ELDA 25 END DO !ELDA 26 2000 FORMAT(//5X,'ELEMENTS CONNECTION'//1X,'ELEM MATE', !ELDA 27

1 10(I2,' NODE')/(10X,10(I2,' NODE'))) !ELDA 28 3000 FORMAT(2I5,10I7/(10X,10I7)) !ELDA 29 RETURN !ELDA 30 END !ELDA 31 C........................................................

SUBROUTINE BOUN(ID) !BOUN 1 C.....RESTRAINT CONDITIONS INPUT !BOUN 2

COMMON /CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !BOUN 3 DIMENSION ID(NDF,*),IDL(6) !BOUN 4 CHARACTER*6 BC !BOUN 5 DATA BC/' B.C. '/ !BOUN 6 WRITE(6,2000) (I,BC,I=1,NDF) !BOUN 7 DO I=1,NJ !BOUN 8

DO J=1,NDF !BOUN 9

ID(J,I)=0 !BOUN 10 END DO !BOUN 11 END DO !BOUN 12

DO K=1,NJ !BOUN 13 READ(5,*) N,NN,NG,(IDL(I),I=1,NDF) !BOUN 14 IF(N.EQ.0) EXIT !BOUN 15 N1=N+(NN-1)*NG !BOUN 16 DO J=N,N1,NG !BOUN 17

DO I=1,NDF !BOUN 18

ID(I,J) =IDL(I) !BOUN 19

END DO !BOUN 20

WRITE(6,3000) J,(ID(I,J),I=1,NDF) !BOUN 21 END DO !BOUN 22 END DO !BOUN 23 2000 FORMAT(/5X,'NODAL BOUNDARY CONDITIONS'//2X,'NODE',7(I4,A6)/1X) !BOUN 24 3000 FORMAT(I6,7(I7,3X)) !BOUN 25 RETURN !BOUN 26 END !BOUN 27 C.........................................................

SUBROUTINE PROFIL(JP,ID,IX,NK) !PROF 1 C.....SET UP EQUATION NUMBERS !PROF 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PROF 3 DIMENSION JP(1),ID(NDF,1),IX(NEN1,*) !PROF 4 https://www.doczj.com/doc/876813557.html,PUTE THE NUMBER OF EQUATION !PROF 5

NEQ=0 !PROF 6 NAD=0 !PROF 7 DO N=1,NJ !PROF 8

DO I=1,NDF !PROF 9

IF(ID(I,N).EQ.0) THEN !PROF 10

NEQ=NEQ+1 !PROF 11

ID(I,N)=NEQ !PROF 12

JP(NEQ)=0 !PROF 13

ELSE !PROF 14

NAD=NAD-1 !PROF 15

ID(I,N)=NAD !PROF 16

ENDIF !PROF 17 END DO !PROF 18 END DO !PROF 19 https://www.doczj.com/doc/876813557.html,PUTE ROW BANDWIDTH !PROF 20

DO N=1,NE !PROF 21 LL=NEQ !PROF 22 DO I=1,NEN !PROF 23

II=IX(I,N) !PROF 24

IF(II.EQ.0) CYCLE !PROF 25

DO K=1,NDF !PROF 26

KK=ID(K,II) !PROF 27

IF(KK.LE.0) CYCLE !PROF 28

LL=MIN0(LL,KK) !PROF 29

END DO !PROF 30 END DO !PROF 31 DO J=1,NEN !PROF 32

JJ=IX(J,N) !PROF 33

IF(JJ.EQ.0) CYCLE !PROF 34

DO L=1,NDF !PROF 35

KK=ID(L,JJ) !PROF 36

IF(KK.LE.0) CYCLE !PROF 37

M=KK-LL !PROF 38

JP(KK)=MAX0(JP(KK),M) !PROF 39

END DO !PROF 40 END DO !PROF 41 END DO !PROF 42 https://www.doczj.com/doc/876813557.html,PUTE DIAGONAL POINTERS FOR PROFIL !PROF 43

NK=1 !PROF 44 JP(1)=1 !PROF 45 DO N=2,NEQ !PROF 46 JP(N)=JP(N)+JP(N-1)+1 !PROF 47 END DO !PROF 48 NK=JP(NEQ) !PROF 49 WRITE(6,2000) NK !PROF 50 2000 FORMAT(' THE MAXIMUM STORAGE OF GLOBAL STIFFNESS =',I6/) !PROF 51 RETURN !PROF 52 END !PROF 53 C...........................................................

SUBROUTINE PLOAD(ID,F,U,NN) !PLOA 1 C.....FORM LOAD ARRAY !PLOA 2

DIMENSION F(1),ID(1),U(1) !PLOA 3 DO N=1,NN !PLOA 4 J=ID(N) !PLOA 5 IF(J.GT.0) U(J)=U(J)+F(N) !PLOA 6 END DO !PLOA 7 RETURN !PLOA 8 END !PLOA 9 C............................................................

SUBROUTINE PFORM(F,LD,UL,XL,S,P,U,X,IX,D,ID,JP,AA,ISW) !PFOR 1 C.....CONTROLING SUBROUTINE OF LOOPING OVER ELEMENT !PFOR 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PFOR 3 COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL,NST !PFOR 4 DIMENSION UL(NDF,*),U(NDF,*),F(NDF,*),S(NST,*),P(*),XL(NDM,*), !PFOR 5

1 X(NDM,*),D(20,*),IX(NEN1,*),ID(NDF,*),LD(NDF,*),JP(*), !PFOR 6

1 AA(*) !PFOR 7

IEL=0 !PFOR 8 MCT=0 !PFOR 9 DO N=1,NE !PFOR 10 CALL PZERO(UL,NST*(NST+3)) !PFOR 11 MA=IX(NEN1,N) !PFOR 12 C.....FORM ELEMENT DATA !PFOR 13

DO I=1,NEN !PFOR 14

II=IX(I,N) !PFOR 15

IF(II.EQ.0) THEN !PFOR 16

DO J=1,NDF !PFOR 17

LD(J,I)=0 !PFOR 18

END DO !PFOR 19

ELSE !PFOR 20

IID=II*NDF-NDF !PFOR 21

NEL=I !PFOR 22

DO J1=1,NDM !PFOR 23

XL(J1,I)=X(J1,II) !PFOR 24

END DO !PFOR 25

DO J=1,NDF !PFOR 26

IF(ISW.EQ.2) THEN !PFOR 27

K=ID(J,II) !PFOR 28

IF(K.LT.0) UL(J,I)=F(J,II) !PFOR 29

LD(J,I)=K !PFOR 30

ELSE !PFOR 31

UL(J,I)=U(J,II) !PFOR 32

LD(J,I)=IID+J !PFOR 33

END IF !PFOR 34

END DO !PFOR 35

END IF !PFOR 36 END DO !PFOR 37 IE20=D(20,MA) !PFOR 38 IF(IE20.NE.IEL) MCT=0 !PFOR 39 IEL=IE20 !PFOR 40 C.....CALL SUBROUTINE OF DIFFERENT ELEMENT TYPE !PFOR 41

CALL ELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,ISW) !PFOR 42 IF(ISW.EQ.2) THEN !PFOR 43 C.....STIFFNESS MATRIX ASSEMBLAGE !PFOR 44

CALL ADDSTF(LD,JP,S,AA,UL,U) !PFOR 45 ELSE !PFOR 46 C.....FORM FORCE OF NODE !PFOR 47

CALL BASBLY(F,P,LD,NST) !PFOR 48 ENDIF !PFOR 49 END DO !PFOR 50 RETURN !PFOR 51 END !PFOR 52 C.............................................................

SUBROUTINE ADDSTF(LD,JP,S,A,UL,U) !ADDS 1 C.....STIFFNESS MATRIX ASSEMBLAGE !ADDS 2

COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL,NST !ADDS 3 DIMENSION LD(*),JP(*),S(NST,NST),A(*),UL(*),U(*) !ADDS 4 DO J=1,NST !ADDS 5

K=LD(J) !ADDS 6

IF(K.GT.0) THEN !ADDS 7

L=JP(K)-K !ADDS 8

DO I=1,NST !ADDS 9

M=LD(I) !ADDS 10

IF(M.GT.K.OR.M.LE.0) CYCLE !ADDS 11

M=L+M !ADDS 12

A(M)=A(M)+S(J,I) !ADDS 13

END DO !ADDS 14 ELSE !ADDS 15

UU=UL(J) !ADDS 16

DO I=1,NST !ADDS 17

M=LD(I) !ADDS 18

IF(M.LE.0) CYCLE !ADDS 19

U(M)=U(M)-S(J,I)*UU !ADDS 20

END DO !ADDS 21 ENDIF !ADDS 22 END DO !ADDS 23

RETURN !ADDS 24 END !ADDS 25 C...........................................................

SUBROUTINE PZERO(P,N) !PZER 1 DIMENSION P(N) !PZER 2 DO I=1,N !PZER 3

P(I)=0.0 !PZER 4 END DO !PZER 5 RETURN !PZER 6 END !PZER 7 C...........................................................

SUBROUTINE BASBLY(B,P,LD,NS) !BASB 1 C.....ASSEMBLE THE LOAD VECTOR !BASB 2

DIMENSION B(*),P(*),LD(*) !BASB 3 DO I=1,NS !BASB 4

II=LD(I) !BASB 5

IF(II.EQ.0) CYCLE !BASB 6

B(II)=B(II)-P(I) !BASB 7 END DO !BASB 8 RETURN !BASB 9 END !BASB 10 C..........................................................

SUBROUTINE LDLT(A,JP,N) !LDLT 1 C.....LDLT TRIANGULAR DECOMPOSITION !LDLT 2

DIMENSION A(*),JP(*) !LDLT 3 NK=JP(N) !LDLT 4 A(1)=1.0/A(1) !LDLT 5 IJ=1 !LDLT 6 II=2 !LDLT 7 ID1=JP(1) !LDLT 8 DO I=2,N !LDLT 9

ID2=JP(I) !LDLT 10 MI=I-ID2+ID1+1 !LDLT 11 MJ=1 !LDLT 12 DO J=MI,I-1 !LDLT 13

IJ=IJ+1 !LDLT 14

IF(J.NE.1) MJ=J-JP(J)+JP(J-1)+1 !LDLT 15

AIJ=A(IJ) !LDLT 16

MIJ=MAX0(MI,MJ) !LDLT 17

IK=JP(I)-I+MIJ !LDLT 18

JK=JP(J)-J+MIJ !LDLT 19

DO K=MIJ,J-1 !LDLT 20

AIJ=AIJ-A(IK)*A(JK) !LDLT 21

IK=IK+1 !LDLT 22

JK=JK+1 !LDLT 23

END DO !LDLT 24

A(IJ)=AIJ !LDLT 25 END DO !LDLT 26 II=IJ+1 !LDLT 27 AII=A(II) !LDLT 28 IJ=JP(I-1)+1 !LDLT 29 DO J=MI,I-1 !LDLT 30

JJ=JP(J) !LDLT 31

AIJ=A(IJ) !LDLT 32

A(IJ)=A(IJ)*A(JJ) !LDLT 33

AII=AII-AIJ*A(IJ) !LDLT 34

IJ=IJ+1 !LDLT 35 END DO !LDLT 36 A(II)=1.0/AII !LDLT 37 ID1=ID2 !LDLT 38 END DO !LDLT 39 RETURN !LDLT 40 END !LDLT 41 C.................................................................

SUBROUTINE FORBACK(A,B,JP,N) !FORB 1 C.....FORWARD REDUCTION !FORB 2

DIMENSION A(*),B(*),JP(*) !FORB 3 IJ=2 !FORB 4 DO I=2,N !FORB 5

MI=I-JP(I)+JP(I-1)+1 !FORB 6

DO J=MI,I-1 !FORB 7

B(I)=B(I)-A(IJ)*B(J) !FORB 8

IJ=IJ+1 !FORB 9

END DO !FORB 10 IJ=IJ+1 !FORB 11 END DO !FORB 12 C.....REFORM THE RIGHT HAND SIDE !FORB 13

DO I=1,N !FORB 14 JPI=JP(I) !FORB 15 B(I)=B(I)*A(JPI) !FORB 16 END DO !FORB 17 C.....BACK SUBSTITUTION !FORB 18

IJ=IJ-1 !FORB 19 DO I=N,2,-1 !FORB 20 MI=I-JP(I)+JP(I-1)+1 !FORB 21 DO J=I-1,MI,-1 !FORB 22

IJ=IJ-1 !FORB 23

B(J)=B(J)-A(IJ)*B(I) !FORB 24 END DO !FORB 25 IJ=IJ-1 !FORB 26 END DO !FORB 27 RETURN !FORB 28 END FORB 29 C..........................................................

SUBROUTINE PRTDIS(ID,U,F) !PRTD 1 C.....PRINT DISPLACEMENTS !PRTD 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PRTD 3 DIMENSION ID(NDF,*),U(*),F(NDF,*) !PRTD 4 WRITE(6,2000)(K,K=1,NDF) !PRTD 5 NN=NJ*NDF !PRTD 6 NM=NEQ !PRTD 7 C.....FORM COMPLETE DISPLEMENT ARRAY !PRTD 8

DO N=NJ,1,-1 !PRTD 9

DO I=NDF,1,-1 !PRTD 10

K=ID(I,N) !PRTD 11

IF(K.LT.0) THEN !PRTD 12

U(NN)=F(I,N) !PRTD 13

ELSE !PRTD 14

U(NN)=U(NM) !PRTD 15

NM=NM-1 !PRTD 16

ENDIF !PRTD 17

NN=NN-1 !PRTD 18 END DO !PRTD 19 END DO !PRTD 20 C.....PRINT DISPLACEMENTS !PRTD 21

K=0 !PRTD 22 DO N=1,NJ !PRTD 23 WRITE(6,3000) N,(U(K+I),I=1,NDF) !PRTD 24 K=K+NDF !PRTD 25 END DO !PRTD 26 RETURN !PRTD 27 2000 FORMAT(//1X,'NODAL DISPLACEMENTS'/2X,'NODE',6(I8,'DISP')) !PRTD 28 3000 FORMAT(I6,6E12.4) !PRTD 29 END !PRTD 30 C.........................................................

SUBROUTINE PRTREAC(R) !PRTR 1 C.....PRINT NODAL REACTIONS !PRTR 2

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !PRTR 3 REAL R(NDF,*),RSUM(6) !PRTR 4 NNEQ=NDF*NJ !PRTR 5 DO K=1,NDF !PRTR 6

RSUM(K)=0.0 !PRTR 7 END DO !PRTR 8 WRITE(6,2000) !PRTR 9 DO I=1,NJ !PRTR 10 DO K=1,NDF !PRTR 11

RSUM(K)=RSUM(K)+R(K,I) !PRTR 12

IF(ABS(R(K,I)).GT.1.0E-3) WRITE(6,3000) I,K,R(K,I) !PRTR 13 END DO !PRTR 14 END DO !PRTR 15 WRITE(6,4000)(RSUM(K),K=1,NDF) !PRTR 16 2000 FORMAT(/5X,'NODAL REACTIONS'/) !PRTR 17 3000 FORMAT(2X,'NODE',I3,I20,'DOF',F12.4) !PRTR 18 4000 FORMAT(3X,'SUM',6F12.4) !PRTR 19 RETURN !PRTR 20 END !PRTR 21 C..........................................................

SUBROUTINE ELEMLIB(D,U,X,IX,S,P,I,J,ISW) !ELEM 1 C.....ELEMENT LIBRARY !ELEM 2

REAL P(NST),S(NST,NST),U(1) !ELEM 3 DIMENSION IX(1),D(1),X(1) !ELEM 4 COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ,NEN1 !ELEM 5 COMMON /ELDATA/ LM,N,MA,MCT,IEL,NEL,NST !ELEM 6 SELECT CASE(IEL) !ELEM 7 CASE(1) !ELEM 8 CALL TRUSS(D,U,X,IX,S,P,I,J,ISW) !ELEM 9

CASE(2) !ELEM 10 CALL PLANE(D,U,X,IX,S,P,I,J,ISW) !ELEM 11 CASE(3:) !ELEM 12 WRITE(6,4000) IEL !ELEM 13 STOP !ELEM 14 END SELECT !ELEM 15 RETURN !ELEM 16 4000 FORMAT(5X,' **FATAL ERROR** ELEMENT CLASS NUMBER',I5,' INPUT') !ELEM 17 END !ELEM 18 C............................................................

SUBROUTINE TRUSS(D,UL,XL,IX,S,P,NDF,NDM,ISW) !TRUS 1 C.... TRUSS LINEAR ELASTIC ELEMENT ROUTINE !TRUS 2

COMMON /ELDATA/ LM,N,MA,MCT,IEL,NEL,NST !TRUS 3 DIMENSION D(10),UL(NDF,1),XL(NDM,1),IX(1),S(NST,1),P(1),DU(2) !TRUS 4 C.... GO TO CORRECT ARRAY PROCESSOR !TRUS 5

SELECT CASE (ISW) !TRUS 6 C.... INPUT MATERIAL PROPERTIES !TRUS 7

CASE(1) !TRUS 8

READ(5,*) E,A !TRUS 9

D(1)=E*A !TRUS 10 WRITE(6,2000) E,A !TRUS 11 C.....CALCULATE ELEMENTAT STIFFNESS MATRIX !TRUS 12

CASE(2,3) !TRUS 13 DX=XL(1,2)-XL(1,1) !TRUS 14 DY=XL(2,2)-XL(2,1) !TRUS 15 DL=SQRT(DX*DX+DY*DY) !TRUS 16 EAL=D(1)/DL !TRUS 17 CC=DX/DL !TRUS 18 SS=DY/DL !TRUS 19 S(1,1)=CC*CC*EAL !TRUS 20 S(2,2)=SS*SS*EAL !TRUS 21 S(1,2)=CC*SS*EAL !TRUS 22 S(2,1)=S(1,2) !TRUS 23 IF(ISW.EQ.2) THEN !TRUS 24

DO I=1,NDF !TRUS 25

DO J=1,NDF !TRUS 26

S(I+2,J) = -S(I,J) !TRUS 27

S(I,J+2) = -S(I,J) !TRUS 28

S(I+2,J+2) = S(I,J) !TRUS 29

END DO !TRUS 30

END DO !TRUS 31 C.....CALCULATE ELEMENTAT NODAL FORCE !TRUS 32

ELSE !TRUS 33

DU(1)=UL(1,1)-UL(1,2) !TRUS 34

DU(2)=UL(2,1)-UL(2,2) !TRUS 35

DO I=1,NDF !TRUS 36

DO J=1,NDF !TRUS 37

P(I)=P(I)-S(I,J)*DU(J) !TRUS 38

END DO !TRUS 39

P(I+2)=-P(I) !TRUS 40

END DO !TRUS 41 C.....CALCULATE ELEMENTAT AXIAL FORCE !TRUS 42

FN=P(1)*CC+P(2)*SS !TRUS 43

WRITE(6,3000) N,IX(1),IX(2),FN !TRUS 44 END IF !TRUS 45 END SELECT !TRUS 46 2000 FORMAT(/5X,'PLANE TRUSS LINEAR ELASTIC ELEMENT'// !TRUS 47

1 10X,7HMODULUS,E18.5/10X,4HAREA,E21.5/) !TRUS 48 3000 FORMAT(/5X,'ELEMENT',I8,5X,'NODE1',I5,5X,'NODE2',I5,5X, !TRUS 49

1 'AXIAL FORCE',E15.6) !TRUS 50

RETURN !TRUS 51 END !TRUS 52 C............................................................

SUBROUTINE PLANE(D,UL,XL,IX,S,P,NDF,NDM,ISW) !PLAN 1 C.... PLANE LINEAR ELASTIC ELEMENT ROUTINE !PLAN 2

COMMON /ELDATA/ LM,N,MA,MCT,IEL,NEL,NST !PLAN 3 DIMENSION D(10),UL(NDF,1),XL(NDM,1),IX(1),S(NST,1),P(1), !PLAN 4

1 SHP(3,9),SG(16),TG(16),WG(16),SIG(6),EPS(3),WD(2) !PLAN 5

CHARACTER WD*4 !PLAN 6 DATA WD/'RESS','RAIN'/ !PLAN 7 C.... GO TO CORRECT ARRAY PROCESSOR !PLAN 8

SELECT CASE (ISW) !PLAN 9 C.... INPUT MATERIAL PROPERTIES !PLAN 10

CASE(1) !PLAN 11 READ(5,*)E,XNU,L,I !PLAN 12 D(4)=L !PLAN 13 IF(I.NE.0) I = 1 !PLAN 14 IF(I.EQ.0) I = 2 !PLAN 15 D(1) = E*(1.+(1-I)*XNU)/(1.+XNU)/(1.-I*XNU) !PLAN 16 D(2) = XNU*D(1)/(1.0E0+(1-I)*XNU) !PLAN 17 D(3) = E/2.0E0/(1.0E0+XNU) !PLAN 18 WRITE(6,2000) WD(I),E,XNU,L !PLAN 19 C.... STIFFNESS COMPUTATION !PLAN 20

CASE(2) !PLAN 21 L = D(4) !PLAN 22 IF(L*L.NE.LINT) CALL PGAUSS(L,LINT,SG,TG,WG) !PLAN 23 C.... FAST STIFFNESS COMPUTATION, COMPUTE INTEGRALS OF SHAPE FUNCTIONS !PLAN 24

DO L = 1,LINT !PLAN 25

CALL SHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX) !PLAN 26

XSJ = XSJ*WG(L) !PLAN 27 C.... LOOP OVER ROWS !PLAN 28

J1 = 1 !PLAN 29

DO J=1,NEL !PLAN 30

W11 = SHP(1,J)*XSJ !PLAN 31

W12 = SHP(2,J)*XSJ !PLAN 32 C.... LOOP OVER COLUMNS (SYMMETRY NOTED) !PLAN 33

K1 = J1 !PLAN 34

DO K = J,NEL !PLAN 35

S(J1 ,K1 ) = S(J1 ,K1 ) + W11*SHP(1,K) !PLAN 36

S(J1 ,K1+1) = S(J1 ,K1+1) + W11*SHP(2,K) !PLAN 37

S(J1+1,K1 ) = S(J1+1,K1 ) + W12*SHP(1,K) !PLAN 38

S(J1+1,K1+1) = S(J1+1,K1+1) + W12*SHP(2,K) !PLAN 39

K1 = K1 + NDF !PLAN 40

END DO !PLAN 41

平面三角形单元有限元程序设计

. 一、题目 如图1所示,一个厚度均匀的三角形薄板,在顶点作用沿板厚方向均匀分布的竖向载荷。已知:P=150N/m ,E=200GPa ,=0.25,t=0.1m ,忽略自重。试计算薄板的位移及应力分布。 要求: 1. 编写有限元计算机程序,计算节点位移及单元应力。(划分三角形 单元,单元数不得少于30个); 2. 采用有限元软件分析该问题(有限元软件网格与程序设计网格必 须一致),详细给出有限元软件每一步的操作过程,并将结果与程序计算结果进行对比(任选取三个点,对比位移值); 3. 提交程序编写过程的详细报告及计算机程序; 4. 所有同学参加答辩,并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下: 1)整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2)单元分析,建立单元刚度矩阵。 3)整体分析,建立总刚矩阵。 4)建立整体结构的等效节点荷载和总荷载矩阵 5)边界条件处理。 6)解方程,求出节点位移。 7)求出各单元的单元应力。 8)计算结果整理。 一、程序设计 网格划分 如图,将薄板如图划分为6行,并建立坐标系,则

刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵,已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律,继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算:(1)每个单元对应2种单元刚度矩阵中的哪一种; (2)该单元对应总刚度矩阵的那几行哪几列 (3)将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环:(1)最外层:逐行计算 (2)中间层:该行逐个计算 (3)最里层:区分为第 奇/偶 数个计算 单元刚度的集成:[ ][][][][][]' '''''215656665656266256561661e Z e e e Z e Z e e e e k k k K k k k k k k +?++=? =?==?==?=?????? 边界约束的处理:划0置1法 X Y P X Y P

有限元分析程序设计

结构有限元分析程序设计 绪论 §0.1 开设“有限元程序设计”课程的意义和目的 §0.2 课程特点 §0.3 课程安排 §0.4 课程要求 §0.5 基本方法复习 $0.1 意义和目的 1.有限元数值分析技术本身要求工程设计研究人员掌握 1). 有限元数值分析技术的完善标志着现代计算力学的真正成熟和实用化,已在各种 力学中得到了广泛的应用。比如:,已杨为工程结构分析中最得以收敛的技术手段,现代功用大致有: a). 现代结构论证。对结构设计从内力,位移等方面进行优劣评定,从而进 行结构优化设计。 b)可取代部份实验,局部实验+有限元分析,是现代工程设计研究方法的一大 特点。 c)结构的各种功能分析(疲劳断裂,可靠性分析等)都以有限元分析工具作为 核心的计算工具。 2). 有限元数值分析本身包括着理论+技术实现(本身功用所绝定的) 有限元数值分析本身包括着泛函理论+分片插值函数+程序设计 2. 有限元分析的技术实现(近十佘年的事)更依赖于计算机程序设计 有限元分析的技术取得的巨大的成就,从某种意义上说,得益于计算机硬件技术的发展和程序设计技术的发展,这两者的依赖性在当代表现得更加突出。(如可视化技术) 3.从学习的角度,不仅要学习理论,而且要从程序设计设计角度对这些理论的技术实现有 一个深入的了解,应当致力于掌握这些技术实现能力,从而开发它,发展它。(理论本身还有待于进一步完美相应的程序设计必须去开发) 4.程序设计不仅是实现有限元数值分析的工具和桥梁,而且在以下诸方面也有意义: 1). 精通基本概念,深化理论认识; 2). 锻炼实际工程分析,实际动手的能力; 3). 获得以后工作中必备的工具。(作业+老师给元素库) 目的:通过讲述有限元程序设计的技术与技巧,便能达到自编自读的能力。 §0.2 课程特点 总描述:理论+算法+数据结构(程序设计的意义) 理论:有限元算法,构造,步骤,解的等外性,收敛性,稳定性,误差分析 算法;指求解过程的技术方法,含两方面的含义;a. 有限元数值分析算法,b, 与数据结构有关的算法(总刚稀疏存贮,提取,节点优化编号等) 数据结构:指各向量矩阵存贮管理与实现,辅助管理结构(指针,数据记录等) 具体特点: 理论性强:能量泛函理论+有限元构造算法+数据结构构造算法 内容繁杂:理论方法+技术方法+技术技巧 技巧性强:排序,管理结构(指针生成,整型运算等)

《有限元》教学大纲

《有限元分析》课程教学大纲 【课程编号】XXXXX 【课程名称】有限元分析/ Finite Element Analysis 【课程性质】专业核心课 【学时】144学时【实验/上机学时】144学时 【考核方式】试卷考【开课单位】XX学院 【授课对象】本科、机械设计制造及其自动化学生 一、课程的性质、目的和任务 有限元法作为边值问题的近似计算方法,随着计算机和计算技术的迅猛发展,其应用已从固体力学发展到流体力学、热力学、电磁学、声学、光学、生物学等多耦合场问题。《有限元分析基础》是材料成型类专业的一门专业基础课,主要介绍固体力学有限单元法的基本理论和应用。在对有限单元法的原理、方法进行讲授的同时配以相应的计算算例及大型工程软件的使用示例,加深学生的理解和消化。 课程教学所要达到的目的是:1、有限单元法的基本理论和实施方法;2、掌握工程结构和设备的受力及变形分析技能并最终提高他们的工程设计能力和解决实际问题的能力; 3、利用ANSYS软件上机实践完成两个上机练习:刚架结构有限元分析和三维固体有限元分析; 4、掌握利用有限元的加权残值法求解场问题的概念,重点介绍1维和2维热传导问。 题有限元分析。 二、教学内容、基本要求和学、课时分配 第一章:ANSYS概论(13学时) (一)基本要求:了解有限元法的分析过程,ANSYS 15.0的安装与启动,前处理、加载并求解、后处理。 (二)教学内容和课时分配: 1、有限元法的分析过程,ANSYS 15.0的安装与启动(2学时) 2系统要求、设置运行参数(1学时)

3、A NSY分析的基本过程(1学时) 4、实验内容(9学时) 实验1梁的有限元建模与变形分析(1学时) 实验目的和要求: 1)要求选择不同形状的截面分别进行计算; 2)梁截面分别采用以下三种截面; 3)设置计算类型; 重点:有限元法的分析过程,ANSYS 15.0的安装与启动; 难点:ANSYS^析的基本过程; 第二章:图形用户界面(13学时) (一)基本要求:了解ANSY S^件界面下各窗口的功能,具体包括应用命令 菜单、主菜单、工具栏、输入窗口、图形窗口和输出窗口。ANSYS架构及命令,具体 包括简单模型的建立、材料属性输入、单元的选择和划分、求解处理和后置处理。 (二)教学内容和课时分配: 1、A NSYS 15.0图形用户界面的组成(1学时) 具体包括应用命令菜单、主菜单、工具栏、输入窗口、图形窗口和输出窗口。ANSYS 架构及命令,具体包括简单模型的建立、材料属性输入、单元的选择和划分、求解处理和后置处理2、对话框及其组件、通用菜单,输入窗口 2、主菜单,输出窗口,图形窗口的功能(1学时) 3、个性化界面(1学时) 4、实验内容(10学时) 实验1超静定桁架的有限元建模与分析 实验目的和要求:上机熟悉ANSY歎件的命令,并对简单的例题进行有限元静、动态分析。 重点(黑体,小四号字):ANSYS 15.0图形用户界面的组成; 难点(黑体,小四号字):主菜单,输出窗口,图形窗口的功能;

有限元法的基本思想及计算 步骤

有限元法的基本思想及计算步骤 有限元法是把要分析的连续体假想地分割成有限个单元所组成的组合体,简称离散化。这些单元仅在顶角处相互联接,称这些联接点为结点。离散化的组合体与真实弹性体的区别在于:组合体中单元与单元之间的联接除了结点之外再无任何关联。但是这种联接要满足变形协调条件,即不能出现裂缝,也不允许发生重叠。显然,单元之间只能通过结点来传递内力。通过结点来传递的内力称为结点力,作用在结点上的荷载称为结点荷载。当连续体受到外力作用发生变形时,组成它的各个单元也将发生变形,因而各个结点要产生不同程度的位移,这种位移称为结点位移。在有限元中,常以结点位移作为基本未知量。并对每个单元根据分块近似的思想,假设一个简单的函数近似地表示单元内位移的分布规律,再利用力学理论中的变分原理或其他方法,建立结点力与位移之间的力学特性关系,得到一组以结点位移为未知量的代数方程,从而求解结点的位移分量。然后利用插值函数确定单元集合体上的场函数。显然,如果单元满足问题的收敛性要求,那么随着缩小单元的尺寸,增加求解区域内单元的数目,解的近似程度将不断改进,近似解最终将收敛于精确解。 用有限元法求解问题的计算步骤比较繁多,其中最主要的计算步骤为: 1)连续体离散化。首先,应根据连续体的形状选择最能完满地描述连续体形状的单元。常见的单元有:杆单元,梁单元,三角形单元,矩形单元,四边形单元,曲边四边形单元,四面体单元,六面体单元以及曲面六面体单元等等。其次,进行单元划分,单元划分完毕后,要将全部单元和结点按一定顺序编号,每个单元所受的荷载均按静力等效原理移植到结点上,并在位移受约束的结点上根据实际情况设置约束条件。 2)单元分析。所谓单元分析,就是建立各个单元的结点位移和结点力之间的关系式。现以三角形单元为例说明单元分析的过程。如图1所示,三角形有三个结点i,j,m。在平面问题中每个结点有两个位移分量u,v和两个结点力分量F x,F y。三个结点共六个结点位移分量可用列

《有限单元法》编程作业

湖南大学 《有限单元法》编程大作业 专业:土木工程 姓名: 学号: 2013年12月

目录 程序作业题目: (3) 1、程序编制总说明 (3) 2、Matlab程序编制流程图 (3) 3、程序主要标示符及变量说明 (4) 4、理论基础和求解过程 (5) 4.1、构造插值函数 (5) 4.2位移插值函数及应变应力求解 (5) 5.程序的验证 (6) 附录:程序代码 (15)

程序作业题目: 完成一个包含以下所列部分的完整的有限元程序( Project) 须提供如下内容的文字材料(1500字以上): ①程序编制说明; ②方法的基本理论和基本公式; ③程序功能说明; ④程序所用主要标识符说明及主要流程框图; ⑤ 1~3 个考题:考题来源、输出结果、与他人成果的对比结果(误差百分比); ⑥对程序的评价和结论(包括正确性、适用范围、优缺点及其他心得等)。 须提供源程序、可执行程序和算例的电子文档或文字材料。选题可根据各自的论文选题等决定。 1、程序编制总说明 a.该程序采用平面三角形等参单元,能解决弹性力学的平面应力、平面应变问题。 b.能计算单元受集中力的作用。 c.能计算结点的位移和单元应力。 d.考题计算结果与理论计算结果比较,并给出误差分析。 e.程序采用MATLAB R2008a编制而成。 2、Matlab程序编制流程图

图1 整个程序流程图 3、程序主要标示符及变量说明 1、变量说明: Node ------- 节点定义 gElement ---- 单元定义 gMaterial --- 材料定义,包括弹性模量,泊松比和厚度 gBC1 -------- 约束条件 gNF --------- 集中力 gk------------总刚 gDelta-------结点位移 输入结构控制参数 输入其它数据 形成整体刚度阵 引入支承条件 解方程,输出位移 求应力,输出应力 形成节点荷载向量 开始 结束 1 单元面积 求弹性矩阵 单元刚度矩阵 位移-应变矩阵 6 7 8 9 10 2 3 4 5

有限单元法与有限元分析

有限单元法与有限元分析 1.有限单元法 在数学中,有限元法(FEM,Finite Element Method)是一种为求解偏微分方程边值问题近似解的数值技术。求解时对整个问题区域进行分解,每个子区域都成为简单的部分,这种简单部分就称作有限元。它通过变分方法,使得误差函数达到最小值并产生稳定解。类比于连接多段微小直线逼近圆的思想,有限元法包含了一切可能的方法,这些方法将许多被称为有限元的小区域上的简单方程联系起来,并用其去估计更大区域上的复杂方程。它将求解域看成是由许多称为有限元的小的互连子域组成,对每一单元假定一个合适的(较简单的)近似解,然后推导求解这个域总的满足条件(如结构的平衡条件),从而得到问题的解。这个解不是准确解,而是近似解,因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解,而有限元不仅计算精度高,而且能适应各种复杂形状,因而成为行之有效的工程分析手段。 随着电子计算机的发展,有限单元法是迅速发展成一种现代计算方法。它是50年代首先在连续体力学领域--飞机结构静、动态特性分析中应用的一种有效的数值分析方法,随后很快广泛的应用于求解热传导、电磁场、流体力学等连续性问题。 1.1.有限元法分析本质 有限元法分析计算的本质是将物体离散化。即将某个工程结构离散为由各种单元组成的计算模型,这一步称作单元剖分。离散后单元与单元之间利用单元的节点相互连接起来;单元节点的设置、性质、数目等应视问题的性质,描述变形形态的需要和计算精度而定(一般情况单元划分越细则描述变形情况越精确,即越接近实际变形,但计算量越大)。所以有限元中分析的结构已不是原有的物体或结构物,而是同新材料的由众多单元以一定方式连接成的离散物体。这样,用有限元分析计算所获得的结果只是近似的。如果划分单元数目非常多而又合理,则所获得的结果就与实际情况相符合。 1.2.特性分析 1)选择位移模式: 在有限单元法中,选择节点位移作为基本未知量时称为位移法;选择节点力作为基本未知量时称为力法;取一部分节点力和一部分节点位移作为基本未知量时称为混合法。位移法易于实现计算自动化,所以,在有限单元法中位移法应用范围最广。 当采用位移法时,物体或结构物离散化之后,就可把单元总的一些物理量如

《有限元分析与应用课程标准》

《有限元分析及应用》课程标准 课程代码:汽车学分:3 建议课时数:64 英文名称: 适用专业:计算机辅助设计与分析 先修课程:《计算机辅助设计》 课程团队负责人及成员:陈良萍、刘宏强、王云、赵静、李蕾、黄艺、史俊玲、 毛新 1.课程定位和设计思路 1.1课程定位 本课程是为计算机辅助设计与分析专业本科生开设的一门专业核心课程,重点介绍有限元法的基本原理和方法、一些成熟的有限元软件功能和简单的分析步骤,同时结合工程实际,为他们进一步学习或实际应用及参加科研工作开辟道路。其任务是通过先修课程中所学知识的综合运用和新知识的获取,使学生初步掌握现代设计中的一种重要方法,开阔视野,提高能力,以适应科学技术发展的要求。 1.2设计思路 在教学中,首先通过力学中的矩阵位移法思想的对比教学,引出连续介质力学有限单元法的学习重点在于单元的插值函数如何构造。这因为,虽说矩阵位移法是对杆系结构而言的,但其结构的离散化和组建整体刚度方程的思想完全可以借鉴到连续介质力学,它们的不同点只是在单元刚度矩阵的建立;而不同单元类型的单元刚度矩阵的建立,又取决于对应单元插值函数的构造。这样处理,不但使学生抓住了本课程的教学重点,而且对有限单元法的整体思想有了宏观上掌握;起到主动学习而非被动接受的作用。在单元构造的教学中,理论学习的重点在于常规单元的介绍;通过常规单元介绍插值函数的完备性与收敛性等。接之,介绍高次单元、等参单元等教学内容。在理论教学中,强调数学论证的严谨性和工程应用的适应性。

结合工程实例教学,拓宽学生数值分析方面的应用能力在课内对不同的单元类 型进行介绍时,及时抓住不同单元在应用中的对比教学与其适用性,并结合工程实例介绍单元类型的合理选取和单元网格的合理划分等。为学生在实际问题的数值分析中如何选定单元和剖分单元奠定了一定的基础和经验。 2.工作任务和课程目标 2.1工作任务 由于采用有限单元法的分析计算软件大多已商业化,而熟悉应用这些中的常规软件也应是本门课程的主要教学内容。在课内学生学会使用软件建立分析模型的基本步骤,其中包括分析模型抽象、几何模型绘制、单元网格划分、材料定义、边界条件定义、方程求解方法等。因课内教学时数的不足,学生应利用课余时间学习,以提高对实际问题的数值分析能力。 2.2课程目标 从教学思想和方法上对原课程进行改革,使学生从较高层次上理解有限元方法的实质,掌握有限元分析的工具,并具备初步处理工程问题的能力;使该课程成为具有较宽口径和较大覆盖面的、面向计算机辅助设计方面的专业基础课;注意课程体 系的整体优化,强调课程的深度、广度与应用。 3.教学方针落实情况

有限元分析软件比较分析

有限元分析软件 有限元分析是对于结构力学分析迅速发展起来的一种现代计算方法。它是50 年代首先在连续体力学领域--飞机结构静、动态特性分析中应用的一种有效的数值分析方法,随后很快广泛的应用于求解热传导、电磁场、流体力学等连续性问题。 有限元分析软件目前最流行的有:ANSYS、ADINA、ABAQUS、MSC 四个比较知名比较大的公司,其中ADINA、ABAQUS 在非线性分析方面有较强的能力目前是业内最认可的两款有限元分析软件,ANSYS、MSC 进入中国比较早所以在国内知名度高应用广泛。目前在多物理场耦合方面几大公司都可以做到结构、流体、热的耦合分析,但是除ADINA 以外其它三个必须与别的软件搭配进行迭代分析,唯一能做到真正流固耦合的软件只有ADINA。ANSYS是商业化比较早的一个软件,目前公司收购了很多其他软件在旗下。ABAQUS专注结构分析目前没有流体模块。MSC是比较老的一款软件目前更新速度比较慢。ADINA是在同一体系下开发有结构、流体、热分析的一款软件,功能强大但进入中国时间比较晚市场还没有完全铺开。 结构分析能力排名:ABAQUS、ADINA、MSC、ANSYS 流体分析能力排名:ANSYS、ADINA、MSC、ABAQUS 耦合分析能力排名:ADINA、ANSYS、MSC、ABAQUS 性价比排名:最好的是ADINA,其次ABAQUS、再次ANSYS、最后MSC ABAQUS 软件与ANSYS 软件的对比分析: 1.在世界范围内的知名度:两种软件同为国际知名的有限元分析软件,在世界范围内具有各自广泛的用户群。ANSYS 软件在致力于线性分析的用户中具有很好的声誉,它在计算机资源的利用,用户界面开发等方面也做出了较大的贡献。ABAQUS软件则致力于更复杂和深入的工程问题,其强大的非线性分析功能在设计和研究的高端用户群中得到了广泛的认可。由于ANSYS 产品进入中国市场早于ABAQUS,并且在五年前ANSYS 的界面是当时最好的界面之一,所以在中国,ANSYS 软件在用户数量和市场推广度方面要高于ABAQUS。但随着ABAQUS北京办事处的成立,ABAQUS软件的用户数目和市场占有率正在大幅度和稳步提高,并可望在今后的几年内赶上和超过ANSYS。 2.应用领域:ANSYS 软件注重应用领域的拓展,目前已覆盖流体、电磁场和多物理场耦合等十分广泛的研究领域。ABAQUS 则集中于结构力学和相关领域研究,致力于解决该领域的深层次实际问题。 3.性价比:ANSYS 软件由于价格政策灵活,具有多种销售方案,在解决常规的

有限元程序课程设计

重庆大学本科学生课程设计任务书 课程设计题目有限元程序设计 学院资源及环境科学学院专业工程力学年级2010级 已知参数和设计要求: 1.独立完成有限元程序设计。 2.独立选择计算算例,并能通过算例判断程序的正确性。 3.独立完成程序设计报告,报告内容包括理论公式、程序框图、程序本 体、计算算例,算例结果分析、结论等。 学生应完成的工作: 1.复习掌握有限单元法的基本原理。 2.掌握弹性力学平面问题3节点三角形单元或4节点等参单元有限元方法 的计算流程,以及单元刚度矩阵、等效节点载荷、节点应变、节点应力 和高斯积分等的计算公式。 3.用Fortran语言编写弹性力学平面问题3节点三角形单元或4节点等参 单元的有限元程序。 4.在Visual Fortran 程序集成开发环境中完成有限元程序的编辑和调试 工作。 5.利用编写的有限元程序,计算算例,分析计算结果。 6.撰写课程设计报告。 目前资料收集情况(含指定参考资料): 1.王勖成,有限单元法,北京:高等教育出版社,2002。 2.O.C. Zienkiewicz, R. L. Taylor, Finite Element Method, 5th Eition, McGraw-Hall Book Company Limited, 2000。 3.张汝清,董明,结构计算程序设计,重庆:重庆大学出版社,1988。 课程设计的工作计划: 1.第1周星期一上午:教师讲解程序设计方法,程序设计要求和任务安 排。 2.第1周星期一至星期二完成程序框图设计。 3.第1周星期三至第2周星期四完成程序设计。 4.第2周星期五完成课程设计报告。 任务下达日期 2013 年 6 月 6 日完成日期 2013 年 07 月 03 日 指导教师(签名) 学生(签名)

有限元法基础试题

有限元法基础试题(A ) 一、填空题(5×2分) 1.1单元刚度矩阵e T k B DBd Ω = Ω? 中,矩阵B 为__________,矩阵D 为___________。 1.2边界条件通常有两类。通常发生在位置完全固定不能转动的情况为_______边界,具体指定有限的非零值位移的情况,如支撑的下沉,称为_______边界。 1.3内部微元体上外力总虚功: ()(),,,,e x x xy y bx xy x y y by d W F u F v dxdy δστδτσδ??=+++++??+(),,,,x x y y xy y x u v u u dxdy σδσδτδδ??+++??的表达式中,第一项为____________________的虚功,第二项为____________________的虚功。 1.4弹簧单元的位移函数1N +2N =_________。 1.5 ij k 数学表达式:令j d =_____,k d =_____,k j ≠,则力i ij F k =。 二、判断题(5×2分) 2.1位移函数的假设合理与否将直接影响到有限元分析的计算精度、效率和可靠性。( ) 2.2变形体虚功原理适用于一切结构(一维杆系、二维板、三位块体)、适用于任何力学行为的材料(线性和非线性),是变形体力学的普遍原理。 ( ) 2.3变形体虚功原理要求力系平衡,要求虚位移协调,是在“平衡、协调”前提下功的恒等关系。 ( ) 2.4常应变三角单元中变形矩阵是x 或y 的函数。 ( ) 2.5 对称单元中变形矩阵是x 或y 的函数。 ( ) 三、简答题(26分) 3.1列举有限元法的优点。(8分) 3.2写出有限单元法的分析过程。(8分) 3.3列出3种普通的有限元单元类型。(6分) 3.4简要阐述变形体虚位移原理。(4分) 四、计算题(54分) 4.1对于下图所示的弹簧组合,单元①的弹簧常数为10000N/m ,单元②的弹簧常数为20000N/m ,单元③的弹簧常数为10000N/m ,确定各节点位移、反力以及单元②的单元力。(10分) 4.2对于如图所示的杆组装,弹性模量E 为10GPa ,杆单元长L 均为2m ,横截面面积A 均为2×10-4m 2,弹簧常数为2000kN/m ,所受荷载如图。采用直接刚度法确定节点位移、作用力和单元②的应力。(10分)

(完整word版)有限元分析软件的比较

有限元分析软件的比较(购买必看)-转贴 随着现代科学技术的发展,人们正在不断建造更为快速的交通工具、更大规模的建筑物、更大跨度的桥梁、更大功率的发电机组和更为精密的机械设备。这一切都要求工程师在设计阶段就能精确地预测出产品和工程的技术性能,需要对结构的静、动力强度以及温度场、流场、电磁场和渗流等技术参数进行分析计算。例如分析计算高层建筑和大跨度桥梁在地震时所受到的影响,看看是否会发生破坏性事故;分析计算核反应堆的温度场,确定传热和冷却系统是否合理;分析涡轮机叶片内的流体动力学参数,以提高其运转效率。这些都可归结为求解物理问题的控制偏微分方程式,这些问题的解析计算往往是不现实的。近年来在计算机技术和数值分析方法支持下发展起来的有限元分析(FEA,Finite Element A nalysis)方法则为解决这些复杂的工程分析计算问题提供了有效的途径。在工程实践中,有限元分析软件与CAD系统的集成应用使设计水平发生了质的飞跃,主要表现在以下几个方面: 增加设计功能,减少设计成本; 缩短设计和分析的循环周期; 增加产品和工程的可靠性; 采用优化设计,降低材料的消耗或成本; 在产品制造或工程施工前预先发现潜在的问题; 模拟各种试验方案,减少试验时间和经费; 进行机械事故分析,查找事故原因。 在大力推广CAD技术的今天,从自行车到航天飞机,所有的设计制造都离不开有限元分析计算,FEA在工程设计和分析中将得到越来越广泛的重视。国际上早20世纪在50年代末、60年代初就投入大量的人力和物力开发具有强大功能的有限元分析程序。其中最为著名的是由美国国家宇航局(NASA)在1965年委托美国计算科学公司和贝尔航空系统公司开发的NASTRAN有限元分析系统。该系统发展至今已有几十个版本,是目前世界上规模最大、功能最强的有限元分析系统。从那时到现在,世界各地的研究机构和大学也发展了一批规模较小但使用灵活、价格较低的专用或通用有限元分析软件,主要有德国的ASKA、英国的PA FEC、法国的SYSTUS、美国的ABQUS、ADINA、ANSYS、BERSAFE、BOSOR、COSMOS、ELAS、MARC和STARDYNE等公司的产品。 以下对一些常用的软件进行一些比较分析: 1. LSTC公司的LS-DYNA系列软件

有限元法及应用---教学大纲

《有限元法及应用》教学大纲 课程代码:050142016 课程英文名称:Application of finite element method 课程总学时:24 讲课:24 实验:0 上机:0 适用专业:材料成型及控制工程 大纲编写(修订)时间:2017.07 一、大纲使用说明 (一)课程的地位及教学目标 (1)课程地位:本课程是材料成型及控制工程专业的专业基础课,必修。 (2)教学目标:掌握材料加工中的有限元法及应用的基础理论与实现方法;能够结合本专业知识进行材料加工中的有限元技术应用;具备初步分析和解决材料加工中的有限元法应用问题的能力。 (二)知识、能力及技能方面的基本要求 (1)知识方面的基本要求: 掌握现代材料加工有限元法的基本概念和建立有限元数学模型的过程与方法。了解有限元法在材料加工领域的各种典型应用。 掌握有限元法的基本概念、解题思路、求解步骤,MSC.MARC有限元分析软件在材料加工中的具体应用技术。了解平面问题与轴对称问题的基本理论。 (2)能力方面的基本要求 具备使用有限元法解决材料加工过程中实际问题的基本能力; 初步具备应用有限元技术分析和解决材料加工过程中各种工艺问题的能力; 具有利用本课程基本理论知识进行材料加工工程中的计算机模拟应用研发与进一步深入学习的能力。 (3)技能方面的基本要求 能够使用有限元软件进行材料加工过程模拟。 (三)实施说明 本教学大纲依据专业教学性指导性计划制定,指导教学环节。 理论教学环节: 教学以课堂讲授为主,多媒体辅助教学,加强基础知识、基本技能、创新意识的培养。结合材料加工领域中的实际案例讲解有限元技术及其在材料加工领域的各种典型应用。对课程中的重点、难点问题着重讲解。由于本课程既具有理论性又具有实践性,因此在教学过程中要注意理论联系实际,通过实例锻炼学生分析解决问题的能力。 采用启发式教学,培养学生思考问题、分析问题和解决问题的能力;通过互动教学调动学生学习的主观能动性,培养学生的独立思考解决问题的能力。 引进与本课有关的发展前沿课题成果,让学生了解本学科的最新发展动态,扩大学生知识视野。 (四)对先修课的要求 在讲授本课前,学生应修完材料成型原理、计算机程序设计。本课程为毕业设计、创新创业训练等实践环节打下基础。 (五)对习题课、实验环节的要求

有限元法中的几个基本概念

诚信·公平·开放·共赢 Loyalty Fair Opening Win-win 有限元法中的几个基本概念 有限元法是把要分析的连续体假想地分割成有限个单元所组成的组合体,简称离散化。 这些单元仅在顶角处相互联接,称这些联接点为结点。 离散化的组合体与真实弹性体的区别在于:组合体中单元与单元之间的联接除了结点之外再无任何关联。但是这种联接要满足变形协调条件,即不能出现裂缝,也不允许发生重叠。显然,单元之间只能通过结点来传递内力。 通过结点来传递的内力称为结点力,作用在结点上的荷载称为结点荷载。当连续体受到外力作用发生变形时,组成它的各个单元也将发生变形,因而各个结点要产生不同程度的位移,这种位移称为结点位移。 在有限元中,常以结点位移作为基本未知量。并对每个单元根据分块近似的思想,假设一个简单的函数近似地表示单元内位移的分布规律,再利用力学理论中的变分原理或其他方法,建立结点力与位移之间的力学特性关系,得到一组以结点位移为未知量的代数方程,从而求解结点的位移分量。然后利用插值函数确定单元集合体上的场函数。显然,如果单元满足问题的收敛性要求,那么随着缩小单元的尺寸,增加求解区域内单元的数目,解的近似程度将不断改进,近似解最终将收敛于精确解。 附:FELAC 2.0软件简介 FELAC 2.0采用自定义的有限元语言作为脚本代码语言,它可以使用户以一种类似于数学公式书写和推导的方式,非常自然和简单的表达待解问题的微分方程表达式和算法表达式,并由生成器解释产生完整的并行有限元计算C程序。 FELAC 2.0的目标是通过输入微分方程表达式和算法之后,就可以得到所有有限元计算的程序代码,包含串行程序和并行程序。该系统采用一种语言(有限元语言)和四种技术(对象技术、组件技术、公式库技术生成器技术)开发而成。并且基于FELAC 1.0的用户界面,新版本扩充了工作目录中右键编译功能、命令终端输入功能,并且丰富了文本编辑功能,改善了用户的视觉体验,方便用户快速便捷的对脚本或程序进行编辑、编译与调试。其中并行版在前后处理上进行了相应的改进。

《有限元》教学大纲

《有限元分析》课程教学大纲 【课程编号】××××× 【课程名称】有限元分析/ Finite Element Analysis 【课程性质】专业核心课 【学时】144学时【实验/上机学时】144学时 【考核方式】试卷考【开课单位】XX学院 【授课对象】本科、机械设计制造及其自动化学生 一、课程的性质、目的和任务 有限元法作为边值问题的近似计算方法,随着计算机和计算技术的迅猛发展,其应用已从固体力学发展到流体力学、热力学、电磁学、声学、光学、生物学等多耦合场问题。《有限元分析基础》是材料成型类专业的一门专业基础课,主要介绍固体力学有限单元法的基本理论和应用。在对有限单元法的原理、方法进行讲授的同时配以相应的计算算例及大型工程软件的使用示例,加深学生的理解和消化。 课程教学所要达到的目的是:1、有限单元法的基本理论和实施方法;2、掌握工程结构和设备的受力及变形分析技能并最终提高他们的工程设计能力和解决实际问题的能力;3、利用ANSYS软件上机实践完成两个上机练习:刚架结构有限元分析和三维固体有限元分析;4、掌握利用有限元的加权残值法求解场问题的概念,重点介绍1维和2维热传导问。 题有限元分析。 二、教学内容、基本要求和学、课时分配 第一章:ANSYS概论(13学时) (一)基本要求:了解有限元法的分析过程,ANSYS 15.0的安装与启动,前处理、加载并求解、后处理。 (二)教学内容和课时分配: 1、有限元法的分析过程,ANSYS 15.0的安装与启动(2学时)

2、系统要求、设置运行参数(1学时) 3、ANSYS分析的基本过程(1学时) 4、实验内容(9学时) 实验1 梁的有限元建模与变形分析(1学时) 实验目的和要求: 1)要求选择不同形状的截面分别进行计算; 2) 梁截面分别采用以下三种截面; 3) 设置计算类型; 重点:有限元法的分析过程,ANSYS 15.0的安装与启动; 难点:ANSYS分析的基本过程; 第二章:图形用户界面(13学时) (一)基本要求:了解ANSYS软件界面下各窗口的功能,具体包括应用命令菜单、主菜单、工具栏、输入窗口、图形窗口和输出窗口。ANSYS架构及命令,具体包括简单模型的建立、材料属性输入、单元的选择和划分、求解处理和后置处理。 (二)教学内容和课时分配: 1、ANSYS 15.0图形用户界面的组成(1学时) 具体包括应用命令菜单、主菜单、工具栏、输入窗口、图形窗口和输出窗口。ANSYS 架构及命令,具体包括简单模型的建立、材料属性输入、单元的选择和划分、求解处理和后置处理2、对话框及其组件、通用菜单,输入窗口 2、主菜单,输出窗口,图形窗口的功能(1学时) 3、个性化界面(1学时) 4、实验内容(10学时) 实验1 超静定桁架的有限元建模与分析 实验目的和要求:上机熟悉ANSYS软件的命令,并对简单的例题进行有限元静、动态分析。 重点(黑体,小四号字):ANSYS 15.0图形用户界面的组成;

有限元法基础重点归纳(精)

1、有限元这种数值计算方法起源于20世纪50年代中期航空工程中飞机结构的矩阵分析。 2、有限单元法的基本思想:在力学模型上将一个原来连续的物体离散成为有限个具有一定 大小的单元,这些单元仅在有限个节点上相连接,并在节点上引进等效力以代替实际作用于单元上的外力。 3、节点:网格间相互连接的点。 4、边界:网格与网格的交界线。 5、有限元的优点:①理论基础简明,物理概念清晰,且可在不同的水平上建立起对该法的 理解②具有灵活性和适用性,应用范围极为广泛③该法在具体推导运算中,广泛采用了矩阵方法。 6、有限单元法分类(从选择基本未知量的角度:位移法(以节点位移为基本未知量,通用 性广、力法(以节点力、混合法(一部分以节点位移,另一部分以节点力 7、有限元法分析计算的基本步骤:①结构的离散化②单元分析(选择位移模式,建立单元 刚度方程,计算等效节点力③整体分析④求解方程,得出节点位移⑤由节点位移计算单元的应变与应力。 8、单元划分:将某个机械结构划分为由各种单元组成的计算模型。 9、有限元法基本近似性------几何近似。

10、弹性力学的任务:分析弹性体在受外力作用并处于平衡状态下产生的应力、应变和位移状态及其相互关系等。 11、弹性力学假设所研究的物体是连续的、完全弹性的、均匀的、各向同性的、微小变形的和无初应力的 12、外力:体力(分布在物体体积内的力---重力、惯性力、电磁力面力(分布在物体表面上的力---流体压力、接触力、风力 13、应力:物体受外力作用,或由于温度有所改变,其内部发生的内力。σ={ σx σy σz τx τy τz } = [σx σy σz τx τy τz ]T 14、应变:物体受到外力作用时,其形状发生改变时的形变。---长度和角度。 ε={ εx εy εz γx γy γz } = [εx εy εz γx γy γz ]T 15、位移:弹性体在载荷作用下,不仅会发生形变,还将产生位移,即弹性体位置 的移动。 δ={u v w }=[u v w ]T 16:、变形协调条件:设想在变形前,把弹性体分为许多微小立方单元体。变形后,每个单元体都产生任意变形而变成一些六面体。可能发生这样的情况,这些六面体

有限元编程的c++实现算例

有限元编程的c++实现算例 1. #include<> 2. #include<> 3. 4. 5. #define ne 3 #define nj 4 #define nz 6 #define npj 0 #define npf 1 #define nj3 1 2 #define dd 6 #define e0 #define a0 #define i0 #define pi 16. 17. 18. int jm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}}; /*gghjghg*/ 19. double gc[ne+1]={,,,}; 20. double gj[ne+1]={,,,}; 21. double mj[ne+1]={,a0,a0,a0}; 22. double gx[ne+1]={,i0,i0,i0}; 23. int zc[nz+1]={0,1,2,3,10,11,12}; 24. double pj[npj+1][3]={{,,}}; 25. double pf[npf+1][5]={{0,0,0,0,0},{0,-20,,,}}; 26. double kz[nj3+1][dd+1],p[nj3+1]; 27. double pe[7],f[7],f0[7],t[7][7]; 28. double ke[7][7],kd[7][7]; 29. 30. 31. 36. void jdugd(int); 38. void zb(int); 39. void gdnl(int); 40. void dugd(int);

有限元分析 教学大纲

《有限元分析》课程教学大纲 一、课程的地位、目的和任务 本课程地位: 《有限元分析》课程是机械设计制造及其自动化专业的一门重要专业选修课。有限元分析方法是一种数值分析方法,在大型数值运算中得到广泛的应用。 本课程目的: 《有限元分析》课程在教学内容方面着重机械分析的基本知识、基本理论和基本方法的传授。在培养学生的设计能力方面着重设计构思和设计技能的基本训练。 本课程任务: 1.树立正确的设计思想和创新意识,了解本课程基本理论的创立、运用和发展; 2.了解国家当前的有关技术、经济政策,具有正确运用标准、规范、手册、图册和查阅有关技术资料的能力; 二、本课程与其它课程的联系 本课程应在学完《画法几何与机械制图》、《理论力学》、《材料力学》课程等课程以后进行,可与《互换性与技术测量》课程同时开设。本课程学习结束后,为学生顺利进入后续专业课学习打下基础,本课程在机械类专业教学计划中起到承前启后的作用,是一门设计性的主干技术课程。在整个人才培养中有不可或缺的总要作用。 三、教学内容及要求 第一篇总论 第一章绪论 教学要求: (1)了解有限元研究的内容与方法; (2)初步理解其在解决固体力学与结构分析方面的问题,而且应用与传热学、流体力学、电磁学等领域的重要地位。 教学内容: 第一节机械结构设计与有限元分析的关系 (一)有限元方法的提出 (二)有限元方法的重要性 第二节用有限元分析方法解决一些工程上的问题 (一)有限元法在工程中的应用

第二章弹性力学的基本理论 教学要求: (1)重点掌握真实解释一个函数,基函数是一组函数,试探函数是某一类函数。教学内容: 第一节有限元相关的数学与力学的知识 (一)有限元数学方程 (二)有限元力学方程 第二节弹性力学变分原理 (一)弹性力学原理 (二)弹性力学的表达式 第三章连续体弹性问题的有限元分析原理 教学要求: (1)掌握该原理; (2)熟知几种常用的单元的节点参数、表达形式和使用范围。 教学内容: 第一节二维、三维建模的有限元分析技术 (一)二维建模有限元技术 (一)三维建模有限元技术 第二节连续体的离散过程 (一)连续体的离散过程 (二) 2D单元的构造 (三) 3D单元的构造 第四章软件使用及结构分析实例与应用教学要求: (1)掌握软件的使用方法,结构问问题的分析与过程; (2)能够应用软件进行一般的结构分析。 教学内容: 第一节分析方法 (一)掌握该种分析方法 (二)解决处理实际工程问题 第二节实践练习 (一)上机练习,尽快掌握分析的原理 第五章接触问题的有限元分析 教学要求: (1)掌握边界接触问题法人解决方法和分析思路。 教学内容: 第一节接触问题的分析方法

有限元单元法复习资料

1.1有限单元法中“离散”的含义是什么?有限单元法是如何将具 有无限自由度的连续介质的问题转变为有限自由度问题的?位移有限单元法的标准化程式是怎样的?(1)离散的含义即将结构离散化,即用假想的线或面将连续体分割成数目有限的单元,并在其上设定有限个节点;用这些单元组成的单元集合体代替原来的连续体,而场函数的节点值将成为问题的基本未知量。(2)给每个单元选择合适的位移函数或称位移模式来近似地表示单元内位移分布规律,即通过插值以单元节点位移表示单元内任意点的位移。因节点位移个数是有限的,故无限自由度问题被转变成了有限自由度问题。(3)有限元法的标准化程式:结构或区域离散,单元分析,整体分析,数值求解。 1.2单元刚度矩阵和整体刚度矩阵各有哪些性质?各自的物理意义是什么?两者有何区别? 单元刚度矩阵的性质:对称性、奇异性(单元刚度矩阵的行列式为零)。整体刚度矩阵的性质:对称性、奇异性、稀疏性。 单元刚度矩阵Kij物理意义Kij即单元节点位移向量中第j个自由度发生单位位移而其他位移分量为零时,在第i个自由度方向引起的节点力。 整体刚度矩阵K中每一列元素的物理意义是:要迫使结构的某节点位移自由度发生单位位移,而其他节点位移都保持为零的变形状态,在所有个节点上需要施加的节点荷载。 2.1 为了使计算结果能够收敛于精确解,位移函数需要满足什么条件?为什么? 满足完备性和协调性。 原因:完备性包括两个条件:即刚体位移条件与常应变条件。首先,位移函数必须包含单元的刚体位移。结构中的单元不仅产生与该单元本身变形相应的位移,还可能因其他单元变形而通过节点位移产生单元刚体位移。为了正确反映单元的实际位移形态,位移函数必须具有反映刚体位移的能力。 其次,由于单元位移函数采用多项式,故在单元内部协调条件总能满足,要求反映在相邻单元之间。实质上来说,要求相邻单元间协调是为了保证单元交界面上应变有限。 3.1构造单元形函数有那些基本原则?试采用构造单元几何方法,构造T10单元的形函数,并对其收敛性进行讨论。 答:形函数是定义于单元内坐标的连续函数。通常单元位移函数采用多项式,其中的待定常数由节点位移参数确定,因此其个数应与单元节点自由度数相等。根据实体结构的几何方程,单元的应变是位移的一次导数。为了反映单元刚体位移和常应变即满足完备性要求,位移函数中必须包含常数项和一次项,即完全一次多项式。 3.3 何谓面积坐标?其特点是什么? 答:三角形单元中,任一点P(x,y)与其3个角点相连形成3个子三角形,其位置可以用下述称为面积坐标的三个比值来确定: L1=A1/A L2=A2/A L3=A3/A 其中A1,A2,A3分别为P23,P31,P12的面积。 各三角形面积为:Ai=1/2* =(ai+bi+ci)/2 由于A1+A2+A3=A,所以有L1+L2+L3=1,Li=(ai+bi+ci)/(2A) 特点:①T3单元的形函数Ni就是面积坐标Li ②面积坐标与三角形在整体坐标系中的位置无关,故称为局部坐标。 ③三个节点的面积坐标分别为节点1(1,0,0),节点2(0,1,0),节点3(0,0,1),形心的面积坐标(1/3,1/3 ,1/3)。④单元边界方程为Li=0 (i=1,2,3); ⑤在平行于2,3边的一条直线上,所有点都要相同的面积坐标。⑥面积坐标与直角坐标互为线性关系。 体积坐标:P点与四面体四个面围成的四个子四面体的体积与原来四面体体积的比值。即 剪切闭锁现象:当梁的高度与梁的长度之比t/l趋于零时,这种单元将出现这种现象,算得的挠度趋于零。 为克服剪切闭锁,使C0型单元适用于各种高度的梁。采用减缩积分方案与假设减应变法。 零能模式:对应于某种非刚体位移模式,减缩积分时高斯点上的应变正好等于零,此时的应变能当然也为零,这种非刚体位移模式称为零能模式。采用减缩积分时会发生零能模式。 5.1、等参单元:将整体坐标系中xy中形状中较复杂的真实单元变换成局部坐标系xy中规则的标准单元,然后在标准单元中构造形函数。由于坐标变换式与单元位移函数中用了相同的形函数N i(ξ,η),故称这种变换为等参变换,相应的单元称为等参单元。 2、等参单元的优越性:①有些工程较复杂,用直边单元离散这些结构需要大量的单元才能得到较好的近似,而曲边的等参单元可方便地离散复杂结构。②如在单元内多取些节点,单元便具有较多的位移自由度,从而就能够插值表示较复杂的单元内部位移场,这样就提高了单元本身的精度。③等参单元刚度矩阵、荷载矩阵的计算是在规则单元域内进行的因此不管被积函数多么复杂,都可以方便地采用标准化数值积分。 3、数值积分的阶次:对于N点积分,当被积函数为m次多项式且m<=2N-1时,可得精确积分值。反之,对于m次多项式的被积函数,精确积分要求的积分点数N>=(m+1)/2。 6.1、工程梁和剪切梁的基本假设?有两种梁弯曲理论①工程梁理论基本假设:平截面假设与横向纤维无挤压假设。前者认为梁横截面变形后仍为平面,且垂直于变形后的中性轴。该假设意味着横向剪切应变γxy =0,后者认为梁的横向纤维无挤压,即εy=0。②剪切梁理论基本假设:横向纤维元无挤压与另一假设认为法平面变形后仍为平面,但不再垂直于变形后的中性轴。 6.2. 剪切梁怎么考虑剪切影响:在结构单元分析中,可在工程梁单元的基础上考虑剪切变形的影响,也可通过挠度与转角各自独立插值直接构造剪切梁单元。 6.3对于杆系结构单元,为什么要在局部坐标系内建立单元刚度矩阵?为什么还要坐标变换?(1)在局部坐标系内可以更方便的建立单 元刚度矩阵。(2)在整体分析中,对所有单元都应采用同一个坐标系即整体坐标系X Y,否则围绕同一节点的不同单元对节点施加的节点力不能直接相加。因此,在进行整体分析之前,还需要进行坐标转换工作,把局部坐标系中得出的单元刚度方程转换成整体坐标系中的单元刚度方程,从而得出整体坐标系中的单元刚度矩阵。 7.1. 薄板弯曲理论基本假定:第一条:板厚方向的挤压变形可忽略不计,即εz=0,。这项假设类似于梁的横向纤维间无挤压假设。第二条:在板弯曲变形中,中面法线保持为直线且仍为弹性曲面(挠度曲面)的法线。第三条:薄板中面只发生弯曲变形,没有面内的伸缩变形,即中面水平位移(u)z=0=(v) z=0=0. 7.2. 厚板理论基本假设:板的中面法线变形后基本保持为直线,但因横向变形的缘故,该直线不在垂直于变形后的中面。因此,法线绕坐标轴的转角θx、θy不再是挠度的导数,而是独立变量;中面内的线位移和板厚方向的挤压变形也可忽略。 7.3. 薄板、厚板基本假定的不同:薄板:板弯曲变形中,中面法线保持为直线且仍为弹性曲面法线。厚板:板中面法线变形后仍基本保持为直线,但该直线不再垂直于变形后的中面。 7.4. DKT单元:离散Kirchhoff理论的基本思想是在若干离散点上满足Kirchhoff直法线假设。基于这种理论构造薄板单元时,w,θx,,θy 也各自独立插值;然后在若干离散点上引入直法线假设。这样构造的单元叫做DKT单元 8.1. 薄壳单元基本假设:薄壳理论假设:薄壳发生微笑变形时,忽略沿壳体厚度方向的挤压变形;且认为直法线假设成立,即变形后中面法线保持为直线且仍为中面的法线;壳体变形时中面不但发生弯曲,而且面内也将产生面内伸缩变形;折板假设;非耦合假设。 薄壳与薄板理论的假设的异同点:相同点:直法线假设和法向(板厚度方向)的纤维无挤压假设均成立。不同点:薄板中面只发生弯曲变形,没有面内的伸缩变形,即中面水平位移为零,而壳体变形时中面不但发生弯曲,而且也将产生面内伸缩变形。 厚壳分析的假设:变形前后的中曲面法线变形后仍基本保持为直线,但因横向剪切变形的缘故,该直线不再垂直于变形后的中曲面,此外,壳体厚度方向的挤压变形可以忽略。 与厚板理论的假设的 相同点:中面法线变形后仍基本保持为直线,但因横向剪切变形的缘故,该直线不在垂直于变形后的中面。厚度方向的挤压变形忽略不计。不同点:厚板理论的假设中,中面内的线位移可以忽略,而厚壳理论的假设中,中面内的位移不可忽略,并且厚壳的位移场可用中面位移表示。 8.2. 平板型单元:组成的折板系统去代替原来的壳体,由平面应力状态与平板弯曲应力状态加以组合而得壳体的应力单元。 分析这种单元时所提出的假设:理论假设:薄壳发生微笑变形时,忽略沿壳体厚度方向的挤压变形,且认为直法线假设成立,即变形后中面法线保持为直线且仍为中面的法线。,折板假设,非耦合假设。 应用平板型壳单元可能会出现的问题,如何解决:1.单元共面问题,解法:引入唯一边界条件可解方程Ka=P 。2.虚拟旋转刚度,解法:在特殊节点上给以任意的虚拟刚度系数。Kθzθzθzi=0,经坐标变换,整体坐标系中该节点平衡方程将满足唯一解条件。赋予Kθzθz任何值。3.新型平面膜元:在平面膜元角点上增加旋转自由度θz,使其有对应的刚度。 8.3. 面内变形与弯曲变形之间非耦合的假设是针对什么提出来的?试说明单元组装时,面内效应与弯曲效应的耦合将会出现。 答:面内变形与弯曲变形之间非耦合的假设是针对局部坐标系下的单元提出的。 9.1. 减少自由度的措施有哪些?各自基本概念如何? 答:1.恰当利用结构对称性。基本思想:利用结构的对称性,取结构一部分建立有限元模型。根据荷载对称性,分析对称面上的位移状态,以确定对称面上节点的位移边界条件。2.采用子结构技术。基本思想:在大型复杂结构的有限元分析中,可将原结构分成若干区域,每个区域作为一个子结构,这些子结构在其公共边界上互相连接起来。 2. 为什么说位移法中应力解的精度低于位移解? 答:在位移有限单元法中,沿单元边界是连续的,而位移的导数通常不连续,因此,在单元边界上应力是不连续的;基本未知量是位移,而单元应变和应力是由位移求导得到的,因此应力精度低于位移精度。 3. 在无法获得精确解的条件下,如何进行误差估计? 答:有限元解法的误差估计有:残值法,后处理法。后处理法:由于无法获得精确解,一般以修匀后的改进值σ*作为“精确解”进行误差估计,通过与精确值误差范数对比,这样做非常有效。

相关主题
相关文档 最新文档