非线性有限元作业--刘晓蓬--21206181

大连理工大学

2013年“非线性有限元”课程考查

姓 名：刘 晓 蓬

学 号：21206181

指导教师：白 瑞 祥

DALIAN UNIVERSITY OF TECHNOLOGY

一、材料非线性部分

1、基本理论与概念

根据相关有限元列式阐述弹塑性物理非线性分析的特点和程序设计要点

塑性是一种在某种给定载荷下,材料产生永久变形的材料特性,对大多的工程材料来说,当其应力低于比例极限时,应力应变关系是线性的,另外,大多数材料在其应力低于屈服点时表现为弹性行为也就是说当移走载荷时其应变也完全消失。由于屈服点和比例极限相差很小,因此在程序中假定它们相同在应力应变的曲线中低于屈服点的叫作弹性部分,超过屈服点的叫作塑性部分,也叫作应变强化部分。塑性分析中考虑了塑性区域的材料特性。路径相关性,既然塑性是不可恢复的,那么这种问题的就与加载历史有关,这类非线性问题叫作与路径相关的或非保守的。非线性路径相关性是指对一种给定的边界条件可能有多个正确的解,内部的应力-应变分布存在,为了得到真正正确的结果必须按照系统真正经历的加载过程,加载率相关性塑性应变的大小可能是加载速度快慢的函数。如果塑性应变的大小与时间有关,这种塑性叫作率无关性塑性,相反,与应变率有关的性叫作率相关的塑性。大多的材料都有某种程度上的率相关性,但在大多数静力分析所经历的应变率范围两者的应力应变曲线差别不大,所以在一般的分析中可变为是与率无关的。工程应力应变与真实的应力应变,塑性材料的数据一般以拉伸的应力应变曲线形式给出,材料数据可能是工程应力与工程应变,也可能是真实应力与真实应变。大应变的塑性分析一般采用真实的应力应变数据,而小应变分析一般采用工程的应力应变数据。其又服从固有的屈服准则、流动准则和强化准则。

程序设计要遵循一些基本原则,下面的这些原则有助于程序进行精确的塑性分析。1.塑性材料常数必须能够足以描述所经历的应力或应变范围内的材料特性。2.缓慢加载应该保证在一个时间步内最大的塑性应变增量小于5%,一般来说如果Fy是系统刚开始屈服时的载荷,那么在塑性范围内的载荷增量应近似为0.05*Fy(对用面力或集中力加载的情况)或Fy (对用位移加载的情况)。3.当模拟类似梁或壳的几何体时必须有足够的网格密度,为了能够足够的模拟弯曲反应在厚度方向必须至少有二个单元。4.除非那个区域的单元足够大,应该避免应力奇异,由于建模而导致的应力奇异有单点加载或单点约束,凹角,模型之间采用单点连接,单点耦合或接触条件。5.加强收敛性的方法是要使用小的时间步长。6.要满足非线性稳定性的要求,要注意迭代的选取。

2、考题计算

由于该厚壁筒模型是轴对称模型,所以在求解过程中,我们选取1/4模型进行了进行建模分析,具体如下图:

建模时取了柱坐标系下厚壁筒从0°~90°范围内的部分,高度取为100mm,模型完成后进行网格的划分,这里利用了Patran的Mesh Seed功能,通过在径向、周向,高度方向撒种生成Mesh网格,网格划分如上图。

考虑到实体的变形情况,关于模型的边界条件,定义如下:

(1)模型的上、下表面为两个平面,在该两平面上限制z方向的位移为0；

(2)对于模型的内外两圆弧面,为了方便定义边界条件,建立了柱坐标,该两平是延径向变形的,所以ρ坐标是放开的,为了限制模型的刚体移动,这里限制角坐标θ为0。

(3)对于模型两个侧平面,是属于模型的对称面,所以该两平面的单元在垂直于平面的方向上位移为零,这里利用柱坐标,即沿周向的位移为零,所以同样要限制角坐标θ为0。

由于厚壁筒受到均匀内压,所以在施加载荷时选择均布载荷Pressure,大小为

2

p=12.5N/mm,作用在内圆弧表面上。对于材料塑性的定义,首先定义样式模量和泊松比,然后在弹塑性对话框里定义屈服载荷和硬化系数或通过在Stress/Strain Curve栏中添加事先定义的材料属性场来表征弹塑性比例系数m。对于求解分析,求解器选择Nastran进行计算分析,单元属性选择3D Solid属性,分析类型定义为非线性并设置大变形和跟随力及载荷增量步等,以此进行弹塑性非线性求解。 结果分析:

(a)对于理性塑性材料,即硬化系数为0,求解结果如下:

该图为100%载荷作用下模型的应力云图及变形情况。观察可知,筒内壁应力较高且首先达到屈服应力发生塑性变形,沿径向方向向外,各层应力逐渐递减,且外层部分属于弹性变形的范畴,模型某一层为弹塑性变形的分界面。

上图为径向方向应力分布图,其中横坐标为相应节点沿径向距圆心的距离,纵坐标为应力值,坐标范围为0~18Mpa,红色曲线中水平直线部分为应力达到屈服应力,在转折点即为弹塑性变形分界面,斜线部分为弹性变形区域,从曲线分析,结果与上图中的应力云图一致。 应力σｒ和σθ沿径向r的分布曲线分别如下:

σ

ｒ-r

图

σ

θ-r

图

(b) 对于线性强化材料,此处选择E1=0.6E,分析方法同上,结果如下: 应力云图及变形情况如下:

应力σｒ和σθ沿径向r的分布曲线分别如下:

σ

ｒ-r

图

σ

θ-r

图

二、几何非线性部分

1.基本理论与概念

详细阐述等弧长法的特点和求解过程,并根据文献检索,介绍几种有关弧长法的改进计算方法。 等弧长法

等弧长法是跟踪结构非线性平衡路径的一种方法,最初由Riks和Wempner提出,继而由Crisfield和Ramm等人加以改进和发展,目前已成为结构非线性稳定分析中的主要方法, Crisfields提出的柱面等弧长法是目前最流行、最有效的等弧长法,其迭代将沿着半径为ΔL(弧长)的空间柱面进行的,因此被称为柱面等弧长法。

等弧长法,结构非线性静力分析中增量形式的平衡方程为:

?KT??q=???P0????Z? (1)

式中:?KT?为平衡路径上某一点的切线刚度矩阵;?q为位移增量向量;??为荷载增量控制参数;?P0?为参考(基准)荷载向量;??Z?为不平衡力向量。

式(1)中?q及??共N+1个未知数,但只有N个平衡方程,需补充一个约束方程:

C??q，???=0 (2)

不同的弧长法均具有相同的约束方程形式:

((qa)(ij))T?qa?i??2?i?j???i?1(j)2

(j)??2?P0?T?P0??li2 (3)

其中:?在等弧长法中取0,?qa?i的括号中下表a表示累积量,即第i增量步内(j)次迭代前的累积位移。

?qa?i(j)??q?i??qa?i?1 (4)

(j)

迭代求解过程

(1)(i=1,j=1)从初始状态(?0?0,q0?0)开始,用一个初始荷载?P1????1jj??1?j?由经验确定)作一线性解,由方程?KT??q0?1?q1?????1???P0?(此时

(j)?j??P0?(其中

?j??j?。 ?1?j????1?j?,?Zi?j?1??0)解得:q1??q1(2)(i=1,j≥2)求?q1,??1,按以下公式计算:

?2??2??KT?1?1??q1?2?=??1?2??P0????Z?1?1? (5)

式中:?KT?1??KT?1(在同一个增量步内采用等刚度迭代);??Z?1是第1次迭代结束时结构的不平衡力。

迭代结束时判断是否满足收敛准则,如果不满足收敛精度要求,则在原增量步内计算

?1??0??1???1?3?,??1?3?,?直至满足收敛条件,再进入到下一个增量步的迭代计算。

第一增量步内迭代k次结束时的总荷载及总位移分别为:

kk?1????1 q1???q1?j? (6)

j?1j?1?j?(3)(i≥2,j=1)按照前述方法计算??i。然后重复上一步中的计算,只是此时切线刚度阵为

?1??KT?2??KT?(q1)??KT??21???KT??2j? (7)

(4)(i≥2,j≥1)由下式解出增量位移??q?i:

(j)?KT?i?qi?j?=??i?j??P0????Z?i?j?1? (8)

式中:??Z?i式中(?F?int)i平;(?F?int)i?j?1??j?1?为上次迭代结束时的不平衡力,即??Z?i?j?1??(?F?int)i?j?1??(?F?ext)i?j?1???v?B?i??j?1??T????j?1?dv为结构的反力,????j?1?为结构的当前应力水

?j?1???i?j?1??P0?为上次迭代结束时的外荷载。

?j?由于在同一增量步中采用等刚度迭代,所以??q?i可写成:

??q?i?????i?j??qT?i???qZ?i?? (9)

jj式中:??qZ?i?j????K?T?i???Z?i?1?j?1?;?qT?i为切线位移也称参考位移向量,即在单位参

考荷载作用下位移的方向,而其大小是没有意义的;?qZ为由当前不平衡力引起的位移。这样,经过j次迭代后,总位移和总荷载参数变为:

qi???qi?jj?1?j??qi?? (10) ?i?j???i?j?1????i?j? (11) 如此重复迭代,直到满足收敛准则。 约束方程及荷载增量控制参数??i??j??的求解 一般工程结构中,往往荷载比较大而位移却很小,在式(3)中,位移的改变对荷载的影响非常小,使得式中左端第二项处于支配地位,载荷近于不变接近与通常的改进的N-R迭代。为此,Crisfield将式(3)改为下式: ??qn?i?qn?i????qn?i?j?j??T?qn?i???li2 (12) j由此得到的弧长法实际上是控制本增量步内N维位移向量的模保持不变,因此称为等弧长法,由式(9)得: j?1????i?j??qT?i???qZ?i (13) Tj?1??j?这样式(12)可表示为: ??qn?i?j?1????qZ?i???i?j??j??qT?i???qn?i????qZ?i???i??j?j??qT?i??li2 (14) 得到一个以??ij为未知数的一元二次方程: A??i式中: ??j??2?B??i?j??C?0 (15) ??A??qT?i?qT?iB?2?qn?iC?T??j?1????qZ?i?j???qn?i?j?1????qZ?i?j???q????q?TTniTi (16) ?j?1????qZ?iji2?j???l2i若方程(15)的两个根??i??j??和??????,而所要选取的根应使得荷载—位移曲线继续向1前”延伸”,而不能”按原路返回”。为此要使本增量步内相邻两次积累位移向量??qn?i和?j???qn?i?j?1?夹角?的余弦为正值。 cos?1,2???qn?i?j?1????i??j???qT?i???qZ?i1,2?j??T?qn?i?j?1? (17) 上式中取cos?>0时的??i作为(j)次迭代的荷载增量因子。 弧长增量控制 在非线性跟踪分析中,适当选取(调整)每增量步中弧长增量非常重要,要做到这一点,就必须对上一步收敛后的一些信息进行分析,一般采用的弧长增量控制方式是由Bellini提出的方法: ?j?li?li?1?Jd/Ji?1? (18) 式中:li为第i增量步的弧长;Jd为期望迭代次数,下标d表示期望值;Ji?1为第(i-1)步增量步内满足收敛时所经过的迭代次数。如果在规定的最大迭代次数内不能满足收敛准则,则可将li减半,从而减小荷载增量步,从第i-1步的收敛状态重新计算。初始弧长l1则由经验来定。 初始荷载增量控制参数??i由式(12)得 1/2??1??的确定 ?1?l??qn?i2i??T?qn?i?1????ili??1???q?2TTi?qT?i (19) 所以有 ??i?1???1?qT?i?qT?iT (20) ??i??的符号决定了荷载值是增加或是减小,根据Berganetal增量功方法判断 ?Wi??qT?i?P0?i (21) 式中?Wi表示参考荷载在当前位移上做的功。若?Wi?0,表示结构处于强化阶段,能够继续加载,取??i?j?j?0;反之?Wi?0,结构处于软化阶段,取??i???0,即卸载。 T 弧长法的改进 本文在Crisfield弧长法的基础上,提出了两种修正弧长的策略:一是根据结构刚度参数的变化去加权修正弧长;二是利用已知平衡点的信息去外插修正弧长。上述两种修正弧长策略

的结合,形成了一种适用性强、率高的改进弧长法。同时,利用改进弧长法实现了对板/壳结构在任一指定载荷下的求解。大量算例说明改进弧长法不仅对结构后屈曲路径全过程的踪,而且对获得指定载荷的收敛解,在求解效率上都有明显的提高。 一种改进弧长法的策略

结构非线性静力分析中增量形式的平衡控制方程为:

??????FI???R??KT???a (22)

式中 : [ KT ]为切线刚度矩阵,

????a为位移增量向量,Δλ为载荷增量控制参数 , FI为参

考载荷向量, R为残余力向量.运用弧长法求解结构非线性有限元方程时,第i步内第j次迭代后

的累加位移为:

????i???a????i??ajj?1??I?i???a??R?i???ij?aj (23)

jj??R?i?????a?a??i式中 : 为第i步第j次迭代的载荷因子, Ii为参考载荷下的切线位移, 为

由残余力引起的位移增量。

荷载因子Δλ的变化受第i步弧长i控制,一般采用Crisfield提出的迭代控制方程:

l????i??a????i??a而弧长i一般由Bellini提出的方法计算:

jTj?li2 (24)

lli?li?1( Jd/Ji?1)1/2 (25)

式中:i为第i增量步的弧长,Jd为期望迭代次数,Ji-1为第i-1增量步满足收敛时的迭代次数。 加权修正弧长

为了提高计算效率,需要针对具体问题(如结构形状、荷作用方式、大小等)合理地选择步长 (Δλ和Δλd) ,根据公式 (4) ,首先要合理地控制弧长i和次数

llld。即仅仅考虑上一增量步的迭代

Ji?1和人为主观上的期望迭代应次数Jd。而可以反应结构当前非线性变形状态的参数

还有其他,如Bergan提出的用结构刚度参数的变化值来描述结构当前所处的非线性状态。于是对公式(4)作用一个权系数α,

li0??li?1Jd/Ji?1 (26)

?是以结构刚度参数变化为主的权系数。

令

???1?2,其中?=?S/S?S,为结构刚度参数的变化,S为第i增量步的无量纲刚度

i2i?1i?2?1为可反映其他因素影响的权系数,一般可根据经验确定。

参数,?S为预定的铡度参数增量;通过大量算例可知,?S和外插修正弧长

?1的建议数据范围要根据经验选取。

采用弧长法求解时能够实现跟踪结构非线性变形的平衡路径,一般情况下 (除分支屈曲外)结构变形的平衡路径随载荷步的增加是逐渐地缓慢地向前延伸的,因此有可能通过对平衡路径上最近一些点的插值,外推出下一个平衡点的大致的位置,然后在新的位置上开始迭代,往往能很快收敛到真正的平衡点,根据这一思想提出了外插弧长法。

如图1,O点为上一步(第i-1步)收敛点,我们利用已知的前m+1个点的位置进行插值,外推到O′点作为下一步的初始位置,很显然O′要比O点能更靠近平衡点A。它的插值表达式是:

?G??,???0i?i?1x?xl???????G??,???kk?i?m?1?l?i?m?1,l?kxk?xl? (27)

i?1?G??,???,?G??,???i分别是已知第k个

式中:

k0增量步的收敛值向量和第i步预测点的初值向量( N+1维,包括节点位移及载荷因子) ( k =

?G?k,?G?l。xk,x为分

i-m-1,…,i-1) ,以下简写作

别表示路径上已知平衡点 k和预测点i的N+1维空间位置坐标,将第i-m-1步平衡作为计算始点(

0xi?m?1?0.0),则后面m个点(i-m,…,i-1)以及预测

点的位置坐标可由已知平衡点的

?G?k?k?i?m?1,?i?1?算得,外插弧长修正

中通常采用三点抛物线插值。

xk?xk?1??G?k??G?k?1

0?k?i?m,i?1?

0x?xi?1?li (28)

?G?i计算预测1'(图1)的切线刚度矩阵,并重新计算弧长, l式中:i为给定的弧长值,然后利用

(29)

作为以后迭代过程中的弧长约束,迭代步骤和一般的弧长法相同。

li??G?i??G?k?12、考题计算

对下图梁结构进行几何非线性分析，载荷，材料参数，几何参数和有限元软件均自选，给出详细分析报告。

M

P

本文来源：https://www.bwwdw.com/article/d3or.html