搜档网
当前位置:搜档网 › 有限元动力学分析方程及解法

有限元动力学分析方程及解法

有限元动力学分析方程及解法
有限元动力学分析方程及解法

动力分析中平衡方程组的解法

1前言

描述结构动力学特征的基本力学变量和方程与静力问题类似,但所有的变量都是时间的函数。

基本变量

三大类变量(,)i u t ξ、(,)ij t εξ和(,)ij t σξ是坐标位置(,,)x y z ξ和时间t 的函数,一般将其记为()()()i ij ij u t t t εσ。

基本方程

(1) 平衡方程

利用达朗贝尔原理将惯性力和阻尼力等效到静力平衡方程中,有

,()()()()0ij j i i i t b t u

t u t σρν+--= (1) 其中ρ为密度,ν为阻尼系数。

(2) 几何方程

,,1

()(()())2ij i j j i t u t u t ε=+ (2)

(3) 物理方程

()()ij ijkl kl t D t σε= (3)

其中ijkl D 为弹性系数矩阵。

(4) 边界条件

位移边界条件()BC u 为,

()()i i u t u t = 在u S 上 (4)

力的边界条件()BC p 为,

()()ij j i t n p t σ= 在p S 上 (5)

初始条件

0(,0)()i i u t u ξξ== (6)

0(,0)()i i u

t u ξξ== (7)

虚功原理

基于上述基本方程,可以写出平衡方程及力边界条件下的等效积分形式,

,()()0p

ij j i i i ij j i S u

u b u d n p dA δσρνδσΩ∏=---+Ω+-=?? (8) 对该方程右端第一项进行分部积分,并应用高斯-格林公式,整理得,

()()0p

ijkl ij kl i i i i i i i i S D u u u u d b u d p u dA εδερδνδδδΩΩ-++Ω-Ω+=??? (9) 有限元分析列式

单元的节点位移列阵为,

111222()[(),(),(),(),(),()(),(),()]e t k k k U t u t v t w t u t v t w t u t v t w t = (10)

单元内的插值函数为,

(,)()()e t u t N U t ξξ= (11)

其中()N ξ为单元的形状函数矩阵,与相应的静力问题单元的形状函数矩阵完全相同,ξ为单元中的几何位置坐标。

基于上面的几何方程和物理方程及(11)式,将相关的物理量表达为节点位移的关系,有,

(,)[](,)[]()()()()e e t t t u t N U t B U t εξξξξ=?=?= (12)

(,)()()()()e e t t t D DB U t S U t σξεξξ=== (13)

(,)()()e t

u t N U t ξξ= (14) (,)()()e t

u t N U t ξξ= (15) 将(12)-(15)供稿到虚功方程(9)中,有,

[()()()()]()0e e e e e e e T e t t t t

t M U t C U t K U t R t U t δδ∏=++-= (16) 由于()e t U t δ具有任意性,消去该项并简写有,

e e e e e t t t t

U C U KU R ++= (17) 其中,

e e T M N Nd ρΩ=

Ω? (18)

e

e T C N Nd νΩ=Ω? (19)

e e T K B DBd Ω=

Ω? (20)

e M 为单元质量矩阵,e C 为单元阻尼矩阵,e K 为单元刚度矩阵。同样,将单元的各个矩阵进行组装,可形成系统的整体有限元方程,即,

MU

CU KU R ++= (21) 其中M ,C 和K 分别是系统的质量、阻尼和刚度矩阵,R 是外荷载向量,U ,U

和U 分别是有限元分割体的加速度、速度和位移向量。方程(21)是通过考虑在时刻t 的静力平衡而推导出来的。

对静力或动力分析的选择(即在分析中是考虑或忽略与速度及加速度有关的力),一般取决于工程上的判断,其目的在于减少所需要的分析工作量。但是,应该认识到,一个静力分析的假定,应该有理由说明它是正确的,否则,分析的结果就是无意义的。确实,在非线性分析中,采用忽略惯性力和阻尼力的假定,可能严重到难以求得甚至无法求得解答。

在数学上,方程(21)是一个二阶线性微分方程组,原则上可用求解常系数微分方程组的标准过程来求得方程组的解。但是,如果矩阵的阶数很高,则采用求解一般微分方程组的过程可能要付出很高的费用,除非特别利用系数矩阵K ,C 和M 的特殊性质。因此,在实用的有限元分析中,主要对几种有效的方法感兴趣,下面将集中介绍这几种方法。我们所考虑的基本过程,可分为两种求解方法:直接积分法和振型叠加法。初看起来,这两种方法似乎完全不同,但事实上它们有着密切的关系,至于选择这种或那种方法,只取决于它们的数值效果。 2直接积分法

在直接积分中对方程(21)是逐步地进行数值积分的,“直接”的意思是,进行数值积分前没有进行把方程变为另一种形式的变换。实质上,直接积分是基于下面的两个想法,第一个想法是只在相隔t ?的一些离散的时间区间上而不是试图在任一时刻t 上满足方程(21)即包含有惯性力和阻尼力作用的(静力)平衡是在求解区间上的一些离散时刻点上获得的。因此,似乎在静力分析中使用过的所有求解方法,在直接积分法中或许也能有效地使用;第二个想法是假定位移、速度和加速度在每一时间区间t ?内变化。

下面假设分别用0

00U ,U ,U 来表示初始时刻)t (0=的位移、速度和加速度向量为已知,要求出方程(21)从0=t 到T t =的解。在求解时,把时间全程T 划分为几个相等的时间区间t ?(即n /T t =?),所用的积分格式是在时刻t ,?0,

T ,,t t ,t ,,t ?+?2上确定方程的近似解。由于计算下一个时刻的解的算法要考虑到前面各个时刻的解,因此假定在时刻t ,,t , ?0的解为已知,来推导出求时刻t t ?+的解的算法。计算时刻t t ?+的解对于计算自此以后t ?的时刻上的解是有代表意义的,这样就可建立用来计算在所有离散时间点上解的一般算法。

(a ) 中心差分法

若把式(21)的平衡关系看作是一个常系数常微分方程组,便可以用任一有限差分表达式通过位移来近似表示加速度和速度。因此,在理论上,许多不同的有限差分表达式均可使用。但是,我们要求求解格式必须是有效的,这样便只需考虑少数几种计算格式。对某些问题求解是非常有效的一个过程是中心差分法,这个方法假定

{}{}

t t t t t t t t t t t U U t U U U U t U ?+?-?+?-+-?=+-?=21212 (22)

将式(22)代入t 时刻的式(21),可得

t t t t t t U C t M t

U M t K R U C t M t ?-?+??? ???-?-??? ???--=??? ???+?2112211222 (23)

从式(23)我们可以求出t t U ?+。应该注意,t t U ?+的解是基于利用在时刻t 的平衡条件。因此,该积分过程称为显式积分方法,且这样的积分格式在逐步解法中不需要对(有效)刚度矩阵进行分解。另一方面,以后所考虑的Houbolt ,Wilson θ及Newmark 方法,要利用在t t ?+上的平衡条件,因而称为隐式积分方法。

另外还应注意到,应用中心差分法时,t t U ?+的计算包含有t U 和t t U ?-,因此,

计算在时刻t ?的解,必需用一个具体的起始过程。由于0

00U ,U ,U 都是已知的,由关系式(22)可求

02002

U t U t U U t ?+?-=?- (24) 具体计算步骤为

A . 初始计算

1.

形成刚度矩阵K 、质量矩阵M 和阻尼矩阵C 。 2. 计算初始值0

00U ,U ,U 。

3.

选取时间步长t ?,要求cr t t ?≤?(临界值)。 4. 计算系数201t a ?=,t a ?=211,022a a =,231a a =。 5.

计算0300U a U t U U t +?-=?-。 6.

形成有效质量矩阵C a M a M ?10+=。 7. 对M ?作三角分解:T LDL M ?=

B . 每一时间步长内的计算

1. 计算在时刻t 的有效荷载:

()()t

t t t t U C a M a U M a K R R ??-----=102。 2.

计算时刻t t ?+的位移:t t t T R ?U LDL =?+。 3.

必要时,按照式(11.3)计算时刻t 速度和加速度。 假设所考虑的系统没有物理阻尼,即C 是零矩阵,在这种情形下式(23)可简

化为

t t t R ?MU t =??+21 (25) 其中

()()t

t t t t U C a M a U M a K R R ??-----=102 因此,如果质量矩阵是对角形的,则解方程组(11.1)时就不需要进行矩阵的分

解,即只需进行矩阵相乘便可求得右端项的有效荷载向量t

R ?,从而利用 ???

? ???=?+ii )

i (t )

i (t t m t R ?U 2 (26) 可得出位移向量的各个分量,其中)i (t t U ?+和)i (t R ?分别表示向量t t U ?+和t

R ?的第i 个分量,而ii m 是质量矩阵的第i 个对角线元素,并且假定0>ii m 。

如果对总刚度矩阵和质量矩阵都不需进行三角分解,也就不必形成总体的K 和M 。此时,求解式(23)可以在单元一级来解决,然后将每个单元的结果累加即可,即

)U U (M t

U K R R ?t i t t i t i i t t 212-?--=∑∑?- (27)

使用式(26)和(27)形式的中心差分法的优点是很明显的,因为它不需要计算总刚度矩阵和总质量矩阵,求解过程基本上是在单元一级上进行,所需要的内存比较少。如果所有相继的单元刚度矩阵和质量矩阵均相同,则该方法就显得更有效,因为这时只需计算或从后备存贮器上连续读出对应于系统中第一个单元的矩阵。

至于中心差分法的缺点,必需承认,该过程的效果与对角形质量矩阵中采用和忽略通常依赖于速度的阻尼力有关,若只包合一个对角形阻尼矩阵,则仍然可保持在单元一级上进行求解的优点。从实用上看,只能用于对角形质量矩阵的这个缺点通常是不很严重的,因为可以采用足够精细的有限元离散化来使解有良好的精度。

使用中心差分格式的另一个十分重要的考虑是,该积分方法要求时间步长t ?小于一个临界值cr t ?,可由整个单元分割体的刚度和质量的性质来算出cr t ?。更准确地说,要得到一个有效的解必须

πn

cr T t t =?≤? (28)

其中n T 是分析物体的最小周期,n 是单元系统的阶。

要求使用的时间步长t ?小于临界时间步长cr t ?的差分格式,例如中心差分法,称为条件稳定的。如果使用一个大于cr t ?的时间步长,则积分是不稳定的,这意味着由数值积分或在计算机上的舍入所导致的误差都会增大,并且在许多情形下会使响应的计算失去意义,积分稳定性的概念是非常重要的。

上面我们讨论了中心差分方法的主要缺点:其格式是条件稳定的,许多其它积分方法也会是条件稳定的。由于在结构分析中条件稳定方法的有效使用会受到限制,因此在下面几节我们来考虑一些通常使用的无条件稳定的积分格式。无条件稳定的积分格式的有效性可从积分能取得较好精度这一事实显示出来,时间步长t ?可以不按式(28)的要求束选取,在许多情形下t ?可以大于式(28)所允许的t ?几个量级。所有下面讨论的积分方法都是隐式的,即在求解时要求对刚度矩阵K ,或者说得更确切些是对有效刚度矩阵进行三角分解。

(b ) Houbolt 法

Houbolt 积分格式与前面讨论的中心差分方法是有联系的,因为它也是根据标准的有限差分表示式用位移的各个分量来近似地表示加速度和速度的各个分量。Houbolt 积分格式采用了下面的有限差分展开式

{}{}t t t t t t t t t t t t t t t t t t U U U U t

U U U U U t U ?-?-?+?+?-?-?+?+-+-?=-+-?=222291811614521 (29) 为了得到在时刻t t ?+的解,现在我们考虑在时刻t t ?+的式(21)(而不是象中心差分法那样考虑在时刻t 的情形),它给出

t t t t U C t M t

R U K C t M t ??? ???+?+=??? ??+?+??+?+35611222 t t t t U C t M t

U C t M t ?-?-??? ???+?+??? ???+?-222311234 (30) 利用上式求解t t U ?+,需要先知道t U ,t t U ?-和t t U ?-2,虽然知道0

00U ,U ,U 对开始Houbolt 积分格式是有用的,但用其它方法计算t t U ,U ??2会更精确,即需要利用特殊的起始方法。

用来积分式(21)以得出t t U ,U ??2,是用一个不同的积分格式并尽可能是用以t ?的几分之一为时间步长,例如用中心差分格式。

Houbolt 法与中心差分法的一个基本差别在于,刚度矩阵K 是作为所求位移t t U ?+的因子出现的。t t KU ?+项的出现是因为在式(21)中取时刻t t ?+的平衡而不像中心差分法取时刻t 的平衡,因此Houbolt 方法是一个隐式积分格式,而中心差分方法则是一个显式过程。关于积分所用的时间步长t ?不存在临界的时间步长限制,一般t ?可以选得比式(28)所给的中心差分法的步长大一些。

一个值得注意的特点是,如果忽略质量和阻尼的影响,基于Houbolt 方法的逐步积分格式可直接简化为静力分析,但中心差分解法则不能这样做。换言之,若0=C 和0=M ,Houbolt 求解方法就得出与时间有关的荷裁作用下的静力解。

(c ) Wilson θ法

Wilson θ法实质上是线性加速度法的推广,线性加速度法假定加速度从时刻t 到时刻t t ?+为线性变化。Wilson θ方程则假定加速度从时刻t 到时刻t t ?+θ为线性变化,其中01.≥θ,当01.=θ时,这个方法就简化为线性加速度格式。要达到无条件稳定,则必须用371.≥θ,通常用401.=θ。

令τ表示时间的增量,其中t ?≤≤θτ0,于是对从t 到t t ?+θ的时间区间,假定

()t t t t t U U t

U U -?+=?++θτθτ (31)

对上式积分得

()()t t t t t t t t t t t t t U U t

U U U U U U t U U U -?+++=-?++=?++?++θτθτ

θτττθττ6212322 (32) 当t ?=θτ时

()()t t t t t t t t t t t t t U U t U t U U U U t U U 26

222+?+?+=+?+=?+?+?+?+θθθθθθθ (33) 这样,就可以求出

()()

t t t t t t t t t t t t t t U t U U U t U U U t U U t U 22326622?---?=-?--?=

?+?+?+?+θθθθθθθθ (34) 这样,要得到在时刻t t ?+的位移、速度和加速度的解,就只需考虑在时刻t t ?+θ的平衡方程(21)。然而,因为假定加速度为线性变化,故所用的投影荷载向量是线性变化的,即

()t t t t t t R R R R -+=?+?+θθ (35)

具体计算步骤为

A . 初始计算

1.

形成刚度矩阵K 、质量矩阵M 和阻尼矩阵C 。 2.

计算初始值000U ,U ,U 。 3. 选取时间步长t ?。计算系数()2

06

t a ?=θ,t a ?=θ31,122a a =,23t

a ?=θ,θ0

4a a =,θ2

5a a -=,θ316-=a ,27t a ?=,6

28t a ?= 4.

形成有效刚度矩阵C a M a K K ?10++=。 5. 对K ?作三角分解:T LDL K ?=

B . 每一时间步长内的计算

1. 计算在时刻t t ?+θ的有效荷载:

()()()t

t t t t t t t t t t t U a U U a C U U a U a M R R R R ? 312022++++++-+=?+?+θθ。 2.

计算时刻t t ?+θ的位移:t t t t T R ?U LDL ?+?+=θθ。 3. 计算时刻t t ?+的位移、速度和加速度:

()t t t t t t t U a U a U U a U 654++-=?+?+θ

()t t t t t t U U a U U ++=?+?+7

()t t t t t t t U U a U t U U 28++?+=?+?+

正如前面所指出的,Wilson θ法也是一个隐式积分方法,因为刚度矩阵K 是未知位移向量的系数矩阵。还可以注意到,该方法不需要特别的初始过程,因为在时刻t t ?+的位移、速度和加速度只是利用在时刻t 的相同的量来表示。

(d ) Newmark 法

Newmark 积分格式也可以认为是线性加速度法的推广,它所用的假定如下

()[]

2211t U U t U U U t U U U U t t t t t t t t t t t t t ???

????+??? ??-+?+=?+-+=?+?+?+?+ ααδδ (36) 其中α和δ是参数,根据积分的精度和稳定性的要求来确定这两个参数。当2

1

=δ和6

1=α时,式(11.17)相应于线性加速度法(它也可由Wilson θ法取1=θ得到)。Newmark 最初提出以‘恒定—平均—加速度法’作为无条件稳定的格式,在这种情形下,2

1=δ,41=α。 要得到在时刻t t ?+的位移、速度和加速度的解,就只需考虑在时刻t t ?+的平衡方程(21),具体步骤如下:

A . 初始计算

1.

形成刚度矩阵K 、质量矩阵M 和阻尼矩阵C 。 2.

计算初始值000U ,U ,U 。 3. 选取时间步长t ?,参数δ和α。计算积分系数:

50.≥δ,()250250δα+≥..,201t a ?=

α,t a ?=αδ1,t a ?=α12,

1213-=αa ,14-=αδa ,??

? ??-?=225αδt a ,()δ-?=16t a ,t a ?=δ7 4.

形成有效刚度矩阵C a M a K K ?10++=。 5. 对K ?作三角分解:T LDL K ?=

B . 每一时间步长内的计算

1. 计算在时刻t t ?+的有效荷载:

()()t t t t t t t t t t U a U a U a C U a U a U a M R R ? 541320++++++=?+?+。

2.

计算时刻t t ?+的位移:t t t t T R ?U LDL ?+?+=。 3. 计算时刻t t ?+的速度和加速度:

()t t t t t t t U a U a U U a U 320---=?+?+

t t t t t t U a U a U U ?+?+++= 76

应该注意Newmark 法和Wilson θ法在计算机上执行时的密切关系,利用这个关系就有可能在一个简单的计算机程序中方便地使用这两个积分格式。 11.1 振型叠加法

在直接积分法中,质量矩阵是对角形的且无阻尼,则对一个时间步长的运算次数稍多于k nm 2,其中n 和k m 分别是所考虑的刚度矩阵的阶和半带宽。在中心差分法中,刚度矩阵与位移向量相乘需要k nm 2次运算;而在Houbolt 法、Wilson θ法和Newmark 法中,则在每个时间步长求解方程组时约需要k nm 2次运算。有效刚度矩阵的初始三角分解还要求一些附加的运算,并且,如果在分析中使用一致质量矩阵或考虑阻尼矩阵,对其中任一情形,每一时间步长所需要的附加运算次数正比于k nm 。因此,若略去初始计算的运算次数,则整个积分所要求的总运算次数约为s nm k α,其中α与所用的矩阵的性质有关,2≥α,而s 是时间步长的步数。

上面的考虑说明,直接积分所需要的运算次数直接正比于分析中的时间步数。因此,一般说来,当要求较短时间(即n 个时间步效)的响应时,可以预料,使用直接积分法是有效的。但是,如果积分必须对许多时间步数进行,则先把平

衡方程(21)变换,使之能以较少的费用进行逐步求解就可能会更有效。具体地说:由于所需要的运算次数直接正比于刚度矩阵的半带宽k m ,因而k m 的减小会按比例地降低逐步解法的费用。

平衡方程组(21)是将有限元插值函数应用在虚功方程的计算中得到的,而最终的矩阵K ,M 和C 的带宽由有限元节点的编号束决定。因此,有限元网格的决定了系统矩阵的阶数和带宽。虽然我们可以通过重新编排节点号达到减小系统矩阵带宽的目的,但是,用这种方式能够得到的最小带宽是有限度的,因而有必要探索另外的途径,这就是振型叠加法。

利用下面的变换,对有限元分析的位移进行变换,即

)t (PX )t (U = (37)

其中P 是个待确定的方阵,)t (X 是与时间有关的n 阶向量,)t (X 的各个分量称为广义位移。将上式代入式(21),并左乘T P 得到

)t (R ~)t (X K ~)t (X C ~)t (X M ~=++ (38)

其中

MP P M ~T =,CP P C ~T =,KP P K ~T =

实际上,这个变换就是将位移的表达式变为

)t (PX N )t ,z ,y ,x (U )m ()m (= (39)

变换的目的是要得到新的系统刚度、质量和阻尼矩阵M ~,K ~和C ~,使它们的

带宽比原来系统矩阵的小,同时也应适当地选取变换矩阵P 。此外应注意,为了使式(37)所表示的任何向量U 和X 之间有唯一的关系,P 必须是非奇异的。

在理论上可以有许多不同的变换矩阵P 都可以减小原来系统矩阵的带宽,但

在实用上,一个有效的变换矩阵是由无阻尼自由振动平衡方程0=+KU U

M 的解来建立的,该方程的可以假定解为

()0t t sin U -=ωφ (40)

其中φ是n 阶向量,t 是时间变量,0t 是时间常数,而ω是表示向量φ的振动频率的常量。把式(40)代入无阻尼自由振动平衡方程,得到广义特征值问题

φωφM K 2= (41)

从中可以确定φ和ω,特征问题(41)式有n 个特征解()121φω,,()

222φω,,……,

()

n n ,φω2

,其中特征向量与M 是正交的,即 ??

?≠====)

j i (,)j i (,M j T i 01φφ (42) 而 222210n

ωωω≤≤≤≤ (43) 向量i φ称为第i 阶振型向量,i ω是相应的振动频率。应该强调指出,n 个位移解公)n ,,,i )(t t (sin i i 210=-ωφ中的任何一个均满足振动平衡方程。 定义一个矩阵Φ,它的各列为特征向量i φ,并定义一个矩阵2Ω,其对角线上的元素为特征值2i ω,即

[]n φφφ 21=Φ,????????????

??=Ω222

212n ωωω (44) 则式(41)可以写成

2ΦΩ=ΦM K (45)

由于向量M 是正交的,于是

I M ,K T T =ΦΦΩ=ΦΦ2 (46)

很显然,矩阵Φ就是一个合适的P 。利用

)t (X )t (U Φ= (47)

可得到对应于振型广义位移的平衡方程

)t (R )t (X )t (X C )t (X

T Φ=Ω+ΦΦ+2 (48) 根据正交性,可得初始条件

00MU X T Φ=,0

0U M X T Φ= (49) 方程(48)表明,若在分析中不包含阻尼矩阵,当利用有限元系统的自由振动的振型向量组成变换矩阵P 时,则有限元平衡方程组是不耦合的,由于在许多情形下阻尼矩阵不能明显确定而只能近似地考虑阻尼的影响,采用下述阻尼矩阵是合理的:它包含所有要求的影响同时又能有效地求出平衡方程的解。在许多分析中是完全忽略阻尼影响的。

(a ) 忽略阻尼的分析

若在分析中不考虑与速度有关的阻尼影响,方程(11.29)就简化为

)t (R )t (X )t (X

T Φ=Ω+2 (50) 亦即下面的n 个独立的方程

n ,,,i )t (R )t (r )t (r x )t (x T i i i i i i 212=?

??==+φω (51) 其中应注意到,在(51)中的第i 个典型方程是一个具有单位质量和刚度2i ω的单自由度系统的平衡方程,这个系统运动时的初始条件可从式(49)得出

000U M )(x ,MU )(x T i i T i i φφ== (52) (51)中的每一项积分,都可以用上节的直接积分法求出,也可以利用杜哈密(Duhamel )积分求得

?++-=t i i i i i i i i t cos t sin d )t (sin )(r )t (x 01

ωβωαττωτω (53)

式中i i ,βα是由初始条件确定的。一般上式的杜哈密积分要采用数值积分来计算。

对于完整的响应,必须算出式(51)中的所有n 个方程的解,然后将每个振型的响应叠加就得到有限元的节点位移,即

∑=i

i i )t (x )t (U φ (54)

因此,概括地讲,用振型叠加的响应分析要求:首先求解方程(41)的特征值和特征向量;然后求解不耦合的平衡方程组(51);最后,如式(54)所示,将每一个特征向量的响应进行叠加。在分析中,特征向量是有限元分割体自由振动的振型。如前面讨论过的,选择振型叠加分析或用上节所述的直接积分法,仅仅取决于其数值效果,而这两种过程所得的解在所用的时间积分格式的数值误差和计算机的舍入误差两方面部是相同的。

对于许多类型的实际荷载条件来说,为了得到系统真实响应的一个好的近似解,仅需要考虑全部不耦合方程中的一部分就够了。通常仅需用到头p 个平衡方程,即为了得到一个满意的近似解,我们在分析中需要考虑式(51)中有关p ,,,i 21=的方程,其中n p <<。这意味着仅需对(41)式解出最小的p 个特征值及相应的特征向量,同时,在式(54)中仅对头p 个振型进行叠加。

事实上,在振型叠加分析中可以仅需考虑少数几个振型,振型叠加过程比直接积分有效得多。但是,由此也可知,振型叠加的效果与分析中必须考虑的振型

数有关,通常所考虑的结构和荷载的空间分布以及频率成份决定着所用振型的个数。对于地震荷载,有些情形只需考虑十个最低的振型,虽然系统的阶可能大于1000。另一方面,对于冲击或振动荷载,一般要包含更多的振型,可能要n p 3

2=。 当考虑振型叠加分析中所应包含的振型的个数的选取问题时,应该记住的是,找出的是动力方程(21)的一个近似解。因此,如果所考虑的振型的个数不足,方程(21)的解就不可能足够精确。实际上这意味着计算出的近似响应,并不满足包含惯性力的平衡。以p U 表示考虑p 个振型时由振型叠加而算得的响应,则可以用下面的误差估计p ε,来估计任意时刻t 的分析精度。 []22)t (R )t (KU )t (U

M )t (R )t (p p p +-= ε (55)

如果己得到系统平衡方程(21)一个满意的近似解,则)t (p ε在任何时刻t 都是很小的。但必需注意,)t (U p 的求得,一定要基于所考虑的p 个振型中每一个振型的响应的精确计算,因为该方法的唯一误差是由于在分析中没有包含足够的振型所致。

应该注意,由上式算出的误差度量决定着包括惯性力的平衡被满足到什么程度,同时它又是惯性和弹性节点力与节点荷载三者平衡的度量.另一方面,因为不是所有振型向量都被用到,所以我们可以说,p ε是没有包括在振型叠加分析中的那些外荷载向量部分的度量,由此可以看出,在直接积分分析中,p ε在时刻 ,t ,t ,??20总是等于零的。总起来说,假定已精确地解出(51)中互不耦合的方程,用n p <的振型叠加分析中的误差是由于没有使用足够的振型所致,而在直接积分分析时,误差的产生是由于使用了过大的时间步长所致。

还要指出一个更为重要的情况,到目前为止仅讨论了有限元系统平衡方程

(21)的精确解。但是,我们真正的目的是要求得所考虑的结构的真实的精确响应的一个良好的近似解。事实上,有限元分析能够求得真实结构的“精确”频率的上界。通常,有限元分析非常好地逼近于最低的精确频率,而可以预料,在逼近高频和相应振型方面就几乎没有精度可言。振型叠加过程在分析中不必包括对应于有限元系统中可能不精确的高频响应,与直接积分相比,这是该方法的特有优点。因而,假定在有限元分析中准确地算出了所有重要的频率,则就不必计算系统中高振型的响应,也不会严重影响解的精度。

(b ) 考虑阻尼的分析

以特征向量)n ,,,i (i 21=φ为基的有限元系统平衡方程的一般形式由式(48)给出。该式表明,假如略去阻尼的影响,则平衡方程就互不耦合,即可以对每一方程逐个地进行时间积分。对不能略去阻尼影响的系统的分析,我们仍然希望基本上能使用相同的计算过程去处理互不耦合的平衡方程(48)。通常,阻尼矩阵C 不能象形成单元分割体的质量矩阵和刚度矩阵那样,由单无阻尼矩阵形成,而它的用处在于近似表示出在系统响应期间的全部能量的消耗。如果能够假定阻尼是成比例的,则振型叠加分析特别有效,该假定为

ij i i j T i C δξωφφ2= (56)

其中i ξ是振型阻尼参数,ij δ是Kronecker δ符号。因此,利用上式,假定特征向量)n ,,,i (i 21=φ也是C 正交的,则方程(48)简化为下面形式的n 个方程:

)t (r x )t (x )t (x

i i i i i i i =++22ωξω (57) 同样,这个方程可以用前节的直接积分计算,也可以用Duhamel 积分 []?++-=---t i i i i t i )t (i i i t cos t sin e d )t (sin e )(r )t (x i i i i 01

ωβωαττωτωωξτξω (58) 其中21i i i ξωω-=,其它同无阻尼一样。

若阻尼满足关系(56)式,在振型叠加分析中就易于考虑阻尼的影响。然而,在真实的阻尼比)p ,,,i (i 21=ξ为已知,则需要用显式算出矩阵C ,把它代人式

(56)中即可得出确定的阻尼比。如果2=p ,可假定Rayleigh 阻尼为如下形式

K M C βα+= (59)

式中α和β是常数。

例如:假设一个多自由度系统,21=ω,32=ω且在该两个振型中我们分别要求有2%和10%的临界阻尼,即要求0201.=ξ和1002.=ξ,试确定Rayleigh 阻尼的常数。

根据Rayleigh 阻尼,有

K M C βα+=

而利用关系(11.37),可得

i i i i T i )K M (ξωβωαφβαφ22=+=+

将两队阻尼比和相应的频率代入

60

090804..=+=+βαβα 解得3360.-=α,1040.=β。这样阻尼矩阵为

K .M .C 10403360+-=

由于已给出阻尼矩阵,如果采用Rayleigh 阻尼矩阵,我们就可以确定任何值i ω的阻尼比

i

i i ωβωαξ22+=。 在实际的分析中,很可能已知两个以上频率的阻尼比,在这种情形下,可用两个平均值1ξ和2ξ来计算α和β。

例如:假设对一个多自由度系统的近似阻尼指定如下:

019

0150070030021401000400300205432154321==========ωωωωωξξξξξ,,,,.,.,.,.,. 试选择合适的Rayleigh 阻尼常数。

由于只要两对阻尼比与频率,根据已知数据,选择

1712040302211====ωξωξ,.;,.

代入可得

08

428924016..=+=+βαβα 解得014980.=α,014050.=β。这样,阻尼矩阵为K .M .C 014050014980+=,计算时的阻尼比由下式求出

i

i i ..ωωξ20140500149802+= 对于用两个以上的阻尼比,可用柯西(Caughey )级数的方式,来建立阻尼矩阵,具体的步骤可参考相关有限元书籍。但是,大多数采用直接积分的实际分析仍采用Rayleigh 阻尼的假定。Rayleigh 阻尼的一个缺点是高振型衰减大于低振型,因为Rayleigh 常数是按低振型选定的。

实际上,对一个具体结构的分析,其合适的Rayleigh 系数常常可以利用一个典型的同类结构的阻尼特性的已有数据来选取,即,把同样的α和β近似地用在同类的结构分析中。Rayleigh 系数的大小,在很大程度上是由结构材料的能量消

耗特性来决定的。

在以上的讨论中,不论在振型叠加分析或者在直接积分过程中,我们都假定可以利用比例阻尼来适当地表示结构的阻尼特性。在多数分析中,比例阻尼的假定(即式(56)满足)是满足要求的。然而,在那些材料性质有很大变化的结构的分机中,就可能需要使用不成比例的阻尼,例如,在基础与结构相互影响的问题中,可以观察到基础的阻尼显著大于地面上结构的阻尼。在这种情形下,在形成结构的阻尼矩阵时,按不同的结构部分规定不同的Rayleigh系数是合理的,这就导致一个不满足关系式(56)的阻尼矩阵。当把阻尼集中在特定的自由度上时(如在结构的支点上)就遇到另一种阻尼不成比例的情形。

具有不成比例阻尼的有限元平衡方程组,可不加修改地利用上节的直接积分算法进行求解,因为阻尼矩阵的性质并不影响到求解过程的推导。另一方面,对用无阻尼自由振动振型作为基向量的振型叠加分析而言,我们看到,式(48)中T在阻尼不成比例的情形下是一个满矩阵,换言之,以振型向量为基的平衡ΦC

Φ

方程组不再是互不耦合的。

有限元动力学分析知识点汇总

复习目录 一、模型输入、建模 A 输入几何模型 1、两种方法:No defeaturing 和 defeaturing (Merge合并选项、Solid实体选项、Small选项) 2、产品接口。输入IGES 文件的方法虽然很好,但是双重转换过程CAD > IGES > ANSYS 在很多情况下并不能实现100%的转换.ANSYS 的产品接口直接读入“原始”的CAD 文件,解决了上面提到的问题. 3、输入有限元模型。除了实体几何模型外, ANSYS 也可输入由某些软件包生成的有限元单元模型数据(节点和单元)。 B 实体建模 1、定义实体建模:建立实体模型的过程。(两种途径) 1)自上而下建模:首先建立体(或面),对这些体或面按一定规则组合得到最终需要的形状. ?开始建立的体或面称为图元. ?工作平面用来定位并帮助生成图元. ?对原始体组合形成最终形状的过程称为布尔运算 ?总体直角坐标系 [csys,0] 总体柱坐标系[csys,1] 总体球坐标系[csys,2] 工作平面 [csys,4] 2)自下而上建模:按照从点到线,从线到面,从面到体的顺序建立模型。 B 网格划分 1、网格划分三步骤: 定义单元属性、指定网格的控制参数、生成网格 2、单元属性(单元类型 (TYPE)、实常数 (REAL)、材料特性 (MAT)) 3、单元类型 单元类型是一个重要选项,它决定如下单元特性: 自由度(DOF)设置、单元形状、维数、假设的位移形函数。 1)线单元(梁单元、杆单元、弹簧单元) 2)壳用来模拟平面或曲面。 3)二维实体用于模拟实体截面 4)三维实体 ?用于几何属性,材料属性,荷载或分析要求考虑细节,而无法采用更简单的单元进行建模的结构。

有限元动力学分析知识点

有限元动力学分析知识 点

复习目录 一、模型输入、建模 A 输入几何模型 1、两种方法:No defeaturing 和 defeaturing (Merge合并选项、Solid实体选项、Small选项) 2、产品接口。输入IGES 文件的方法虽然很好,但是双重转换过程CAD > IGES > ANSYS 在很多情况下并不能实现100%的转 换.ANSYS 的产品接口直接读入“原始”的CAD 文件,解决了上面提到的问题. 3、输入有限元模型。除了实体几何模型外, ANSYS 也可输入由某些软件包生成的有限元单元模型数据(节点和单元)。 B 实体建模 1、定义实体建模:建立实体模型的过程。(两种途径) 1)自上而下建模:首先建立体(或面),对这些体或面按一定规则组合得到最终需要的形状. ?开始建立的体或面称为图元. ?工作平面用来定位并帮助生成图元. ?对原始体组合形成最终形状的过程称为布尔运算 ?总体直角坐标系 [csys,0] 总体柱坐标系[csys,1] 总体球坐标系[csys,2] 工作平面 [csys,4] 2)自下而上建模:按照从点到线,从线到面,从面到体的顺序建立模型。

B 网格划分 1、网格划分三步骤: 定义单元属性、指定网格的控制参数、生成网格 2、单元属性(单元类型 (TYPE)、实常数 (REAL)、材料特性 (MAT)) 3、单元类型 单元类型是一个重要选项,它决定如下单元特性: 自由度(DOF)设置、单元形状、维数、假设的位移形函数。 1)线单元(梁单元、杆单元、弹簧单元) 2)壳用来模拟平面或曲面。 3)二维实体用于模拟实体截面 4)三维实体 ?用于几何属性,材料属性,荷载或分析要求考虑细节,而无法采用更简单的单元进行建模的结构。 ?也用于从三维CAD系统转化而来的几何模型,而这些几何模型转化成二维模型或壳体会花费大量的时间和精力 4、单元阶次与形函数 ?单元阶次是指单元形函数的多项式阶次。 ?什么是形函数? –形函数是指给出单元内结果形态的数值函数。因为FEA 的解答只是节点自由度值,需要通过形函数用节点自由 度的值来描述单元内任一点的值。 –形函数根据给定的单元特性给出。

多体动力学软件和有限元软件的区别(优.选)

有限元软件与多体动力学软件 数值分析技术与传统力学的结合在结构力学领域取得了辉煌的成就,出现了以ANSYS 、NASTRAN 等为代表的应用极为广泛的结构有限元分析软件。计算机技术在机构的静力学分析、运动学分析、动力学分析以及控制系统分析上的应用,则在二十世纪八十年代形成了计算多体系统动力学,并产生了以ADAMS 和DADS 为代表的动力学分析软件。两者共同构成计算机辅助工程(CAE )技术的重要内容。 商业通用软件的广泛应用给我们工程师带来了极大的便利,很多时候我们不需要精通工程问题中的力学原理,依然可以通过商业软件来解决问题,不过理论基础的缺失还是会给我们带来不少的困扰。随着动力有限元与柔性多体系统分析方法的成熟,有时候正确区分两者并不是很容易。 机械领域应用比较广泛的有两类软件,一类是有限元软件,代表的有:ANSYS, NASTRAN, ABAQUS, LS-DYNA, Dytran 等;另一类是多体动力学软件,代表的有ADAMS, Recurdyn , Simpack 等。在使用时,如何选用这两类软件并不难,但是如果深究这两类软件根本区别并不容易。例如,有限元软件可以分析静力学问题,也可以分析“动力学”问题,这里的“动力学”与多体动力学软件里面的动力学一样吗?有限元软件在分析动力学问题时,可以模拟物体的运动,它与多体动力学软件中模拟物体运动相同吗?多体动力学软件也可以分析柔性体的应力、应变等,这与有限元软件分析等价吗? 1 有限元软件 有限单元法是一种数学方法,不仅可以计算力学问题,还可以计算声学,热,磁等多种问题,我们这里只探讨有限元法在机械领域的应用。 计算结构应力、应变等的力学基础是弹性力学,弹性力学亦称为弹性理论,主要研究弹性体在外力作用或温度变化等外界因素下所产生的应力、应变和位移,从而为工程结构或构件的强度、刚度设计提供理论依据和计算方法。也就是说用有限元软件分析力学问题时,是用有限元法计算依据弹性力学列出的方程。 考虑下面这个问题,在()0t , 时间内给一个结构施加一个随时间变化的载荷()P t ,我们希望得到结构的应力分布,在刚刚施加载荷的时候,结构中的应力会有波动,应力场是变化的,但很久以后,应力场趋于稳定。 如果我们想得到载荷施加很久以后,稳定的应力场分布,那么应该用静力学分析方法分析

有限元动力学分析方程及解法

动力分析中平衡方程组的解法 1前言 描述结构动力学特征的基本力学变量和方程与静力问题类似,但所有的变量都是时间的函数。 基本变量 三大类变量(,)i u t ξ、(,)ij t εξ和(,)ij t σξ是坐标位置(,,)x y z ξ和时间t 的函数,一般将其记为()()()i ij ij u t t t εσ。 基本方程 (1) 平衡方程 利用达朗贝尔原理将惯性力和阻尼力等效到静力平衡方程中,有 ,()()()()0ij j i i i t b t u t u t σρν+--=&&& (1) 其中ρ为密度,ν为阻尼系数。 (2) 几何方程 ,,1 ()(()())2ij i j j i t u t u t ε=+ (2) (3) 物理方程 ()()ij ijkl kl t D t σε= (3) 其中ijkl D 为弹性系数矩阵。 (4) 边界条件 位移边界条件()BC u 为, ()()i i u t u t = 在u S 上 (4) 力的边界条件()BC p 为, ()()ij j i t n p t σ= 在p S 上 (5) 初始条件 0(,0)()i i u t u ξξ== (6) 0(,0)()i i u t u ξξ==&& (7)

虚功原理 基于上述基本方程,可以写出平衡方程及力边界条件下的等效积分形式, ,() ()0p ij j i i i ij j i S u u b u d n p dA δσρνδσΩ∏=---+Ω+-=??&&& (8) 对该方程右端第一项进行分部积分,并应用高斯-格林公式,整理得, ()()0p ijkl ij kl i i i i i i i i S D u u u u d b u d p u dA εδερδνδδδΩΩ-++Ω-Ω+=???&&& (9) 有限元分析列式 单元的节点位移列阵为, 111222()[(),(),(),(),(),()(),(),()]e t k k k U t u t v t w t u t v t w t u t v t w t =L (10) 单元内的插值函数为, (,)()()e t u t N U t ξξ= (11) 其中()N ξ为单元的形状函数矩阵,与相应的静力问题单元的形状函数矩阵完全相同,ξ为单元中的几何位置坐标。 基于上面的几何方程和物理方程及(11)式,将相关的物理量表达为节点位移的关系,有, (,)[](,)[]()()()()e e t t t u t N U t B U t εξξξξ=?=?= (12) (,)()()()()e e t t t D DB U t S U t σξεξξ=== (13) (,)()()e t u t N U t ξξ=&& (14) (,)()()e t u t N U t ξξ=&&&& (15) 将(12)-(15)供稿到虚功方程(9)中,有, [()()()()]()0e e e e e e e T e t t t t t M U t C U t K U t R t U t δδ∏=++-=&&&g (16) 由于()e t U t δ具有任意性,消去该项并简写有, e e e e e t t t t U C U KU R ++=&&& (17) 其中, e e T M N Nd ρΩ= Ω? (18) e e T C N Nd νΩ=Ω? (19)

UG有限元分析教程

第1章高级仿真入门 在本章中,将学习: ?高级仿真的功能。 ?由高级仿真使用的文件。 ?使用高级仿真的基本工作流程。 ?创建FEM和仿真文件。 ?用在仿真导航器中的文件。 ?在高级仿真中有限元分析工作的流程。 1.1综述 UG NX4高级仿真是一个综合性的有限元建模和结果可视化的产品,旨在满足设计工程师与分析师的需要。高级仿真包括一整套前处理和后处理工具,并支持广泛的产品性能评估解法。图1-1所示为一连杆分析实例。 图1-1连杆分析实例 高级仿真提供对许多业界标准解算器的无缝、透明支持,这样的解算器包括NX Nastran、MSC Nastran、ANSYS和ABAQUS。例如,如果结构仿真中创建网格或解法,则指定将要用于解算模型的解算器和要执行的分析类型。本软件使用该解算器的术语或“语言”及分析类型来展示所有网格划分、边界条件和解法选项。另外,还可以求解模型并直接在高级仿真中查看结果,不必首先导出解算器文件或导入结果。 高级仿真提供基本设计仿真中需要的所有功能,并支持高级分析流程的众多其他功能。 ?高级仿真的数据结构很有特色,例如具有独立的仿真文件和FEM文件,这有利于在分布式工作环境中开发有限元(FE)模型。这些数据结构还允许分析师轻松 地共享FE数据去执行多种类型分析。

UG NX4高级仿真培训教程 2 ?高级仿真提供世界级的网格划分功能。本软件旨在使用经济的单元计数来产生高质量网格。结构仿真支持完整的单元类型(1D、2D和3D)。另外,结构级仿真 使分析师能够控制特定网格公差。例如,这些公差控制着软件如何对复杂几何体 (例如圆角)划分网格。 ?高级仿真包括许多几何体简化工具,使分析师能够根据其分析需要来量身定制CAD几何体。例如,分析师可以使用这些工具提高其网格的整体质量,方法是消 除有问题的几何体(例如微小的边)。 ?高级仿真中专门包含有新的NX传热解算器和NX流体解算器。 NX传热解算器是一种完全集成的有限差分解算器。它允许热工程师预测承受热载荷系统中的热流和温度。 NX流体解算器是一种计算流体动力学(CFD)解算器。它允许分析师执行稳态、不可压缩的流分析,并对系统中的流体运动预测流率和压力梯度,也可 以使用NX传热和NX流体一起执行耦合传热/流体分析。 1.2仿真文件结构 当向前通过高级仿真工作流时,将利用4个分离并关联的文件去存储信息。要在高级仿真中高效地工作,需要了解哪些数据存储在哪个文件中,以及在创建那些数据时哪个文件必须是激活的工作部件。这4个文件平行于仿真过程,如图1-2所示。 图1-2仿真文件结构 设计部件文件的理想化复制 当一个理想化部件文件被建立时,默认有一.prt扩展名,fem#_i是对部件名的附加。例如,如果原部件是plate.prt,一个理想化部件被命名为plate_fem1_i.prt。 一个理想化部件是原设计部件的一个相关复制,可以修改它。 理想化工具让用户利用理想化部件对主模型的设计特征做改变。不修改主模型部件,

有限元分析基础教程(ANSYS算例)(曾攀)

有限元分析基础教程Fundamentals of Finite Element Analysis (ANSYS算例) 曾攀 清华大学 2008-12

有限元分析基础教程曾攀 有限元分析基础教程 Fundamentals of Finite Element Analysis 曾攀 (清华大学) 内容简介 全教程包括两大部分,共分9章;第一部分为有限元分析基本原理,包括第1章至第5章,内容有:绪论、有限元分析过程的概要、杆梁结构分析的有限元方法、连续体结构分析的有限元方法、有限元分析中的若干问题讨论;第二部分为有限元分析的典型应用领域,包括第6章至第9章,内容有:静力结构的有限元分析、结构振动的有限元分析、传热过程的有限元分析、弹塑性材料的有限元分析。本书以基本变量、基本方程、求解原理、单元构建、典型例题、MATLAB程序及算例、ANSYS算例等一系列规范性方式来描述有限元分析的力学原理、程序编制以及实例应用;给出的典型实例都详细提供有完整的数学推演过程以及ANSYS实现过程。本教程的基本理论阐述简明扼要,重点突出,实例丰富,教程中的二部分内容相互衔接,也可独立使用,适合于具有大学高年级学生程度的人员作为培训教材,也适合于不同程度的读者进行自学;对于希望在MATLAB程序以及ANSYS平台进行建模分析的读者,本教程更值得参考。 本基础教程的读者对象:机械、力学、土木、水利、航空航天等专业的工程技术人员、科研工作者。

目录 [[[[[[\\\\\\ 【ANSYS算例】3.3.7(3) 三梁平面框架结构的有限元分析 1 【ANSYS算例】4.3.2(4) 三角形单元与矩形单元的精细网格的计算比较 3 【ANSYS算例】5.3(8) 平面问题斜支座的处理 6 【ANSYS算例】6.2(2) 受均匀载荷方形板的有限元分析9 【ANSYS算例】6.4.2(1) 8万吨模锻液压机主牌坊的分析(GUI) 15 【ANSYS算例】6.4.2(2) 8万吨模锻液压机主牌坊的参数化建模与分析(命令流) 17 【ANSYS算例】7.2(1) 汽车悬挂系统的振动模态分析(GUI) 20 【ANSYS算例】7.2(2) 汽车悬挂系统的振动模态分析(命令流) 23 【ANSYS算例】7.3(1) 带有张拉的绳索的振动模态分析(GUI) 24 【ANSYS算例】7.3(2) 带有张拉的绳索的振动模态分析(命令流) 27 【ANSYS算例】7.4(1) 机翼模型的振动模态分析(GUI) 28 【ANSYS算例】7.4(2) 机翼模型的振动模态分析(命令流) 30 【ANSYS算例】8.2(1) 2D矩形板的稳态热对流的自适应分析(GUI) 31 【ANSYS算例】8.2(2) 2D矩形板的稳态热对流的自适应分析(命令流) 33 【ANSYS算例】8.3(1) 金属材料凝固过程的瞬态传热分析(GUI) 34 【ANSYS算例】8.3(2) 金属材料凝固过程的瞬态传热分析(命令流) 38 【ANSYS算例】8.4(1) 升温条件下杆件支撑结构的热应力分析(GUI) 39 【ANSYS算例】8.4(2) 升温条件下杆件支撑结构的热应力分析(命令流) 42 【ANSYS算例】9.2(2) 三杆结构塑性卸载后的残余应力计算(命令流) 45 【ANSYS算例】9.3(1) 悬臂梁在循环加载作用下的弹塑性计算(GUI) 46 【ANSYS算例】9.3(2) 悬臂梁在循环加载作用下的弹塑性计算(命令流) 49 附录 B ANSYS软件的基本操作52 B.1 基于图形界面(GUI)的交互式操作(step by step) 53 B.2 log命令流文件的调入操作(可由GUI环境下生成log文件) 56 B.3 完全的直接命令输入方式操作56 B.4 APDL参数化编程的初步操作57

有限元动力学分析方程及解法

动力分析中平衡方程组的解法 1前言 描述结构动力学特征的基本力学变量和方程与静力问题类似,但所有的变量都是时间的函数。 基本变量 三大类变量(,)i u t ξ、(,)ij t εξ和(,)ij t σξ是坐标位置(,,)x y z ξ和时间t 的函数,一般将其记为()()()i ij ij u t t t εσ。 基本方程 (1) 平衡方程 利用达朗贝尔原理将惯性力和阻尼力等效到静力平衡方程中,有 ,()()()()0ij j i i i t b t u t u t σρν+--= (1) 其中ρ为密度,ν为阻尼系数。 (2) 几何方程 ) ,,1()(()())2 ij i j j i t u t u t ε=+ (2) (3) 物理方程 ()()ij ijkl kl t D t σε= (3) 其中ijkl D 为弹性系数矩阵。 (4) 边界条件 位移边界条件()BC u 为, ()()i i u t u t = 在u S 上 (4) 力的边界条件()BC p 为, ()()ij j i t n p t σ= 在p S 上 (5) 初始条件 0(,0)()i i u t u ξξ== (6)

{ 0(,0)()i i u t u ξξ== (7) 虚功原理 基于上述基本方程,可以写出平衡方程及力边界条件下的等效积分形式, ,()()0p ij j i i i ij j i S u u b u d n p dA δσρνδσΩ∏=---+Ω+-=?? (8) 对该方程右端第一项进行分部积分,并应用高斯-格林公式,整理得, ()()0p ijkl ij kl i i i i i i i i S D u u u u d b u d p u dA εδερδνδδδΩΩ-++Ω-Ω+=??? (9) 有限元分析列式 单元的节点位移列阵为, 111222()[(),(),(),(),(),()(),(),()]e t k k k U t u t v t w t u t v t w t u t v t w t = (10) 单元内的插值函数为, (,)()()e t u t N U t ξξ= (11) % 其中()N ξ为单元的形状函数矩阵,与相应的静力问题单元的形状函数矩阵完全相同,ξ为单元中的几何位置坐标。 基于上面的几何方程和物理方程及(11)式,将相关的物理量表达为节点位移的关系,有, (,)[](,)[]()()()()e e t t t u t N U t B U t εξξξξ=?=?= (12) (,)()()()()e e t t t D DB U t S U t σξεξξ=== (13) (,)()()e t u t N U t ξξ= (14) (,)()()e t u t N U t ξξ= (15) 将(12)-(15)供稿到虚功方程(9)中,有, [()()()()]()0e e e e e e e T e t t t t t M U t C U t K U t R t U t δδ∏=++-= (16) 由于()e t U t δ具有任意性,消去该项并简写有, e e e e e t t t t U C U KU R ++= (17)

多体动力学和非线性有限元联合仿真

A New Solution For Coupled Simulation Of Multi-Body Systems And Nonlinear Finite Element Models Giancarlo CONTI, Tanguy MERTENS, Tariq SINOKROT (LMS, A Siemens Business) Hiromichi AKAMATSU, Hitoshi KYOGOKU, Koji HATTORI (NISSAN Motor Co., Ltd.) 1 Introduction One of the most common challenges for flexible multi-body systems is the ability to properly take into account the nonlinear effects that are present in many applications. One particular case where these effects play an important role is the dynamic modeling of twist beam axles in car suspensions: these components, connecting left and right trailing arms and designed in a way that allows for large torsional deformations, cannot be modeled as rigid bodies and represent a critical factor for the correct prediction of the full-vehicle dynamic behavior. The most common methods to represent the flexibility of any part in a multi-body mechanism are based on modal reduction techniques, usually referred to as Component Mode Synthesis (CMS) methods, which predict the deformation of a body starting from a preliminary modal analysis of the corresponding FE mesh. Several different methods have been developed and verified, but most of them can be considered as variations of the same approach based on a limited set of modes of the structure, calculated with the correct boundary conditions at each interface node with the rest of the mechanism, allowing to greatly reduce the size of system’s degrees of freedom from a large number of nodes to a small set of modal participation factors. By properly selecting the number and frequency range of the modes, as well as the boundary conditions at each interface node [1], it is possible to accurately predict the static and dynamic deformation of the flexible body with remarkable improvements in terms of CPU time: this makes these methods the standard approach to reproduce the flexibility of components in a multi-body environment. Still, an important limitation inherently lies in their own foundation: since displacements based on modal representation are by definition linear, any nonlinear phenomena cannot be correctly simulated. For example, large deformations like twist beam torsion during high lateral acceleration cornering maneuvers typically lead to geometric nonlinearities, preventing any linear solution from accurately predicting most of the suspension’s elasto-kinematic characteristics like toe angle variation, wheel center position, vertical stiffness. One possible solution to overcome these limitations while still working with linear modal reduction methods is the sub-structuring technique [2]: the whole flexible body is divided into sub-structures, which are connected by compatibility constraints preventing the relative motion of the nodes that lie between two adjacent sub-structures. Standard component mode synthesis methods are used in formulating the equations of motion, which are written in terms of generalized coordinates and modal participation factors of each sub-structure. The idea behind it is that each sub-portion of the whole flexible structure will undergo smaller deformations, hence remaining in the linear flexibility range. By properly selecting the cutting sections it is usually possible to improve the accuracy of results (at least in terms of nodal displacements: less accuracy can be expected for stress and strain distribution). Another limitation of these methods is the preliminary work needed to re-arrange the FE mesh, although some CAE products already offer automatic processes enabling the user to skip most of the re-meshing tasks and hence reducing the modeling efforts. An alternative approach to simulate the behavior of nonlinear flexible bodies is based on a co-simulation technique that uses a Multi-body System (MBS) solver and an external nonlinear Finite Element Analysis (FEA) solver. Using this technique one can model the flexible body in the external nonlinear FEA code and the rest of the car suspension system in the MBS environment. The loads due to the deformation of the body are calculated externally by the FEA solver and communicated to the MBS solver at designated points where the flexible body connects to the rest of the multi-body system. The MBS solver, on the other hand, calculates displacements and velocities of these points and communicates them to the nonlinear FEA solver to advance the simulation. This approach doesn’t suffer from the limitations that arise from the linear modeling of the flexibility of a body. This leads to more accurate results, albeit at the price of much larger CPU time. In fact, simulation results are strongly affected by the size of the communication time step between the two solvers: a better accuracy (and more stable solver convergence) can be generally obtained by using smaller time steps which require larger calculation times, as shown also in [3].

有限元分析基础教程

有限元分析基础教程

前言 有限元分析已经在教学、科研以及工程应用中成为重要而又普及的数值分析方法和工具;该基础教程力求提供具备现代特色的实用教程。在教材的内容体系上综合考虑有限元方法的力学分析原理、建模技巧、应用领域、软件平台、实例分析这几个方面,按照教科书的方式深入浅出地叙述有限元方法,并体现出有限元原理“在使用中学习,在学习中使用”的交互式特点,在介绍每一种单元的同时,提供完整的典型推导实例、MATLAB实际编程以及ANSYS应用数值算例,并且给出的各种类型的算例都具有较好的前后对应性,使学员在学习分析原理的同时,也进行实际编程和有限元分析软件的操作,经历实例建模、求解、分析和结果评判的全过程,在实践的基础上深刻理解和掌握有限元分析方法。 一本基础教材应该在培养学员掌握坚实的基础理论、系统的专业知识方面发挥作用,因此,教材不但要提供系统的、具有一定深度的基础理论,还要介绍相关的应用领域,以给学员进一步学习提供扩展空间,本教程正是按照这一思路进行设计的;全书的内容包括两个部分,共分9章;第一部分为有限元分析基本原理,包括第1章至第5章,内容有:绪论、有限元分析过程的概要、杆梁结构分析的有限元方法、连续体结构分析的有限元方法、有限元分析中的若干问题讨论;第二部分为有限元分析的典型应用领域,包括第6章至第9章,内容有:静力结构的有限元分析、结构振动的有限元分析、传热过程的有限元分析、弹塑性材料的有限元分析。在基本原理方面,以基本变量、基本方程、求解原理、单元构建等一系列规范的方式进行介绍;在阐述有限元分析与应用方面,采用典型例题、MATLAB程序及算例、ANSYS算例的方式,以体现出分析建模的不同阶段和层次,引导学员领会有限元方法的实质,还提供有大量的练习题。 本教程的重点是强调有限元方法的实质理解和融会贯通,力求精而透,强调学员综合能力(掌握和应用有限元方法)的培养,为学员亲自参与建模、以及使用先进的有限元软件平台提供较好的素材;同时,给学员进一步学习提供新的空间。 本教程力求体现以下特点。 (1)考虑教学适应性:强调对学员在数学原理、分析建模、软件应用几个方面的培养目标要求,注重学员在工程数值方面的基础训练,培养学员“使用先进软件+分析实际问题”的初步能力。 (2)考虑认知规律性:力求按照有限元分析方法的教学规律和认知规律,在教材中设计了“基本变量、基本方程、求解原理、单元构建”这样的模块;并体现出有限元原理“在使用中学习,在学习中使用”的交互式特点,在介绍每一种单元的同时,提供实用的MATLAB实际编程和数值实例;在每一章还进行要点总结,给出典型例题,以引导学员领会有限元方法的实质,体现教材的启发性,有利于激发学员学习兴趣和便于自学。 (3)考虑结构完整性:本教程提供完整的教材结构:绪论、正文、典型例题、基于MATLAB的编程算例与数值算例、具有一定深度的ANSYS算例、各章要点、习题、专业术语的英文标注、关键词中文和英文索引、参考文献,便于学员查阅。 (4)内容上的拓展性:除基本内容外,还介绍了较广泛的应用领域,包括:静力结构分析、结构振动分析、传热过程分析、弹塑性材料分析;提供了有关的典型问题的建模详细分析过程,基本上反映了有限元分析在一些主要领域的应用状况及建模方法。 (5)编排上的逻辑性:本教程力求做到具有分明的层次和清楚的条理,在每一章中重点突出有限元方法的思想、数理逻辑及建模过程,强调相应的工程概念,提供典型例题及详解,许多例题可作为读者进行编程校验的标准考题(Benchmark),还提供了对应的MATLAB编程算例与ANSYS算例,特别是介绍了基于APDL参数化的ANSYS建模方法,并给出具体的实例,力求反映有限元分析的内在联系及特有思维方式。

ProE Mechanica有限元分析入门教程

Pro/E Mechanica有限元分析入门教程 一、进行Mechanica分析的步骤: 1)建立几何模型:在Pro/ENGINEER中创建几何模型。 2)识别模型类型:将几何模型由Pro/ENGINEER导入Pro/MECHANICA中,此步需要用户确定 模型的类型,默认的模型类型是实体模型。我们为了减小模型规模、提高计算速度,一般用面的形式建模。 3)定义模型的材料属性。包括材料、密度、弹性模量、泊松比等。 4)定义模型的约束。 5)定义模型的载荷。 6)有限元网格的划分:由Pro/MECHANICA中的Auto GEM(自动网格划分器)工具完成有限元 网格的自动划分。 7)定义分析任务,运行分析。 8)根据设计变量计算需要的项目。 9)图形显示计算结果。 二、下面将上述每一步进行详解: 1、在Pro/ENGINEER模块中完成结构几何模型后,单击“应用程序”→“Mechanica”,弹出下 图所示窗口, 点击Continue继续。弹出下图,启用Mechanica Structure。一定要记住不要勾选有限元模式前面的复选框,最后确定。

2、添加材料属性单击“材料”,进入下图对话框,选取“More”进入材料库,选取材料 Name---------为材料的名称; References-----参照Part(Components)-----零件/组件/元件 V olumes-------------------体积/容积/容量; Properties-------属性Material-----材料;点选后面的More就可以选择材料的类型 Material Orientation------材料方向,金属材料或许不具有方向性,但是某些复合材料是纤维就具有方向性,可以根据需要进行设置方向及其转角。点选OK,材料分配结束。 3、定义约束 1):位移约束 点击,出现下图所示对话框,

有限元分析-清华大学教程

8.1 进入工程分析模块 8.2施加约束 8.3 施加载荷 8.4 静态有限元计算过程和后处理 8.5动态分析的前处理和显示计算结果8.6有限元分析实例 习题

工程分析指的是有限元分析,包括静态分析(Static Analyses)和动态分析。动态分析又分为限制状态固有频率分析(Frequency Analyses)和自由状态固有频率分析(Free Frequency Analyses),前者在物体上施加一定约束,后者的物体没有任何约束,即完全自由。 8.1 进入工程分析模块 1. 进入工程分析模块前的准备工作 (1)在三维实体建模模块建立形体的三维模型,为三维形体添加材质,见4.7。 (2)将显示模式设置为Shading(着色)和Materials(材料),这样才能看到形体的应力和变形图,详见2.11.6。

2. 进入工程分析模块 选择菜单【Start】→【Analysis & Simulation】→【Generative Structural Analysis】弹出图8-1所示新的分析实例对话框。 在对话框中选择静态分析(Static Analyses)、限制状态固有频率分析(Frequency Analyses)还是自由状态固有频率分析(Free Frequency Analyses),单击OK按钮,将开始一个新的分析实例。 图8-1新的分析实例对话框

3.有限元分析的过程 有限元分析的一般流程为: (1)从三维实体建模模块进入有限元分析模块。(2)在形体上施加约束。 (3)在形体上施加载荷。 (4)计算(包括网格自动划分),解方程和生成应力应变结果。 (5)分析计算结果,单元网格、应力或变形显示。(6)对关心的区域细化网格、重新计算。 上述(1)~(3)过程是有限元分析预(前)处理,(4)是计算过程,(5)、(6)是有限元后处理。 有限元文件的类型为CATAnalysis。

动力学有限元

6.2结构动力有限元法理论与模型 一、基本原理 在实际问题的求解中,应用最广的是基于位移的有限元素法。此法的基本思想是把本来为连续的工程结构分割成在结点上相联的单元组合体。取这些结点的位移为基本未知量,并假定每个单元中的位移用单元位移函数来描述,这实质上是假定了单元的模态。在此基础上,利用能量变分原理进行单元分析的全结构分析,得到全结构的振动平衡方程,从而把连续体的动力学问题化为多自由度系统的振动问题。有限元动力分析的基本过程是首先将工程结构离散化,通过选择合理的单元确定出分析模型,在此基础上选择位移函数,进行单元分析,确定单元的刚度、质量、阻尼、载荷矩阵,再经过坐标变换,通过能量变分原理,进行全结构分析,建立系统的振动平衡方程。最后运用有限元数值方法进行方程的求解。 结构动力有限元法采用的单元位移函数与静力分析相同,基本原理和求解过程也与静力分析相同,不同之处仅在分析模型的确定与运动方程的建立方面。 二、动态分析模型的确定 由于结构动态分析中除考虑弹性力外,还要考虑惯性力和阻尼力,其运动方程是常微分方程组,所以动态分析的复杂程度高,计算工作量大,有限元分析模型要尽量精炼、简单。 1.模型确定的基本原则 ?分析模型应与分析的目的相适应。动力分析的目的各不相同,有的是为了提供固有特性计算动态响应或供控制系统用;有的是为了舱内提供振动环境。不同的目的,通常要求不同的模态数与计算精度。显然,用于估算基本固 有频率的模型应当比计算冲击响应的模型简单。用于设计计算的模型应当比用于校核计算的模型简单。 ?分析模型要与选用的计算工具与计算条件相适应。计算机软件种类日益丰富,选择分析模型要与所用程序、所用计算机容量相适应。如对于容量大的计算机,可选用较为复杂的有限元模型,而对于容量小的计算机则在能反映 结构动态性能的前提下尽量简化模型,使求解规模尽量小。对于大模型,可选用子结构模型,采用模态综合方法 求解。应注意, 不一定模型愈精细精度就愈高。模型愈复杂,往往带来了更繁杂的运算。虽然模型误差小了,但 计算误差加大。不恰当的精细模型反而得到不佳的结果。 ?模型应正确反映结构的实际特性。一个具体结构的动态特性,主要取决于质量、刚度的大小与分布,取决于结构边界条件与阻尼特性。因此,模型应尽量保持整个飞行器、甚至各部件的质量、质心位置不变,保证构造的刚度 特性和传力路线基本不变。真实反映分析对象的边界条件与阻尼特性。因此,选取的单元应保持其几何形状、受 力特性、变形特性与实际结构的几何形状特点、受力传力特点、变形特性相一致。例如,对于薄壁结构翼面或 舱段,不应选用板、壳单元来模拟蒙皮。

动力学问题有限元分析

【问题描述】如图I所示的固定梁模型,梁的材料为结构钢,尺寸为3m×0.5m×0.025m,结构两端固定。用集中力250N代表旋转的机器,作用点位于梁的三分之一处,机器旋转速率为0~3000rpm。 图I 固定梁 【要求】在ANSYS Workbench软件平台上,对该梁结构进行谐响应分析。 1.分析系统选择 (1)运行ANSYS Workbench,进入工作界面,首先设置模型单位。在菜单栏中找到Units下拉菜单,依次选择Units>Metric(kg,m,s,℃,A,N,V)命令。 (2)在左侧工具箱【Toolbox】下方“分析系统”【Analysis Systems】中双击“模态分析”【Model】系统,此时在右侧的“项目流程”【Project Schematic】中会出现该分析系统共7个单元格。相关界面如图1所示。 图1 分析系统选择

(3)拖动左侧工具箱中“分析系统”【Analysis Systems】中的“谐响应”【Harmonic Response】系统进到模态分析系统的【Solution】单元格中,为之后热应力分析做准备。完成后的相关界面如图2所示。 图2 谐响应分析流程图 2.输入材料属性 (1)在右侧窗口的分析系统A中双击工程材料【Engineering Data】单元格,进入工程数据窗口。 (2)可以看见,系统本身默认结构钢【Structural Steel】已在备选材料窗口中,在此不必再另行选择,直接单击【Project】选项卡回到项目流程界面即可。 3.建立几何模型 (1)双击分析系统A中的“几何”【Geometry】单元格。 (2)在菜单栏中依次选择Units>Meter,确认以“米”作为建模单位。之后,单击树形目录中的【XYPlane】,再单击工具栏中的“创建草图”选项即选项以创建草图,此时【XYPlane】分支下出现了名为“Sketch1”的草绘平面。如图3所示。

相关主题