《结构分析中的有限元法》2015-有限元习题-参考答案

本科有限元习题参考答案

2015年3月10日作业

1、简述力学课程中介绍的各种力学模型的简化条件、基本假设和适用范围（包括有拉压杆模型、弯曲梁模型、平面应力和平面应变模型、轴对称模型、板模型、壳模型等）

2、给出弹性力学问题中平衡方程、几何方程、物理方程的表达式及其意义。

（1）平衡方程：

zy

yz xz zx yx xy z yz

xz z y xy

zy y x zx

yx x f y

x z f x z y f z y x ττττττττσττσττσ====+??+??+??=+??+??+??=+??+??+??,000， 物理意义：应力分量与体力分量之间的关系。 （2）几何方程：

z

u x w y w z v x v y u z

w y v x u zx yz xy

z y x ??+??=??+??=??+??=??=??=??=γγγεεε,,,,

物理意义：应变分量与位移分量之间的关系。 （3）物理方程：

[]

[]

[]

zx

zx yz yz xy xy y x z z z x y y z y x x G

G G E E

E τγτγτγσσμσεσσμσεσσμσε1,1,1)

(1

)

(1

)(1

===+-=+-=+-=

物理意义：应变分量与应力分量之间的关系。

3、简述最小势能原理的主要内容和主要公式。

根据虚功原理得到：??=-Γ

T Ω

T T 0Td Γδu d Ω）F δu -σδε（，由

)(21εδσεδδεU T T =??

?

??=则0)21((=Γ-Ω-=∏??ΩΓ）Td u d F u T T T p σεδδ

其中，??ΩΓ

Γ-Ω-=∏Td u d F u T T T p )21

(σε即为系统的总势能，它是弹性体变

形势能和外力势能之和。上面变分为零式表明：在所有区域内满足几何关系，在边界上满足给定位移条件的可能位移中，真实位移使系统的总势能取驻值(可证

明此驻值为最小值)。

4、“有限元法”都有哪些名称(包括中文和英文)？

有限元法，也叫有限单元法、有限元素法、有限元分析；FEM（Finite Element Method）、FEA（Finite Element Analysis）

5、简述有限元法的发展和现状。

近几十年，伴随着计算机科学和技术的快速发展，有限元法作为工程分析的有效方法在理论、方法的研究、计算机程序的开发以及应用领域的开拓者方面均取得了根本性的发展。

（1）单元的类型和形式

为了扩大有限元法的应用领域，新的单元类型和形式不断涌现(等参元,梁板壳,复合材料)

（2）有限元法的理论基础和离散格式

将Hellinger-Reissner、Hu—Washizu(多场变量变分原理)应用于有限元分析，发展了混合模型、杂交型的有限元表达格式，应研究了各自的收敛条件；将加权余量法用于建立有限元的表达格式；进一步研究发展有限元解的后验误差估计和应力磨平方法。

（3）有限元方程的解法(大型复杂工程结构问题——静态, 特征值, 瞬态等) （4）有限元法的计算机软件(专用软件, 通用软件)

6、结构力学分析问题中的三种“非线性”都包含哪些，并解释其含义。

非线性可以是由材料性质、变形状态和边界接触条件引起的，分别称为材料、几何、边界非线性。材料非线性就是材料的本构关系不是线性的。几何非线性时结构在载荷作用过程中产生大的位移和转动，如板壳结构的大挠度。边界非线性是指高挠度部件或由多个部件组成的结构组合件，渐进位移将会增大部件自身或是部件之间产生接触的可靠性，以此特征的特定类型几何非线性为边界条件或者接触非线性。

7、简述有限元法的未来。

以有限元法为代表的计算力学提出一系列新的课题：

（1）为了真实地模拟新材料和新结构的行为，需要发展新的材料本构模型和单元型式。例如对于特种合金、复合材料、陶瓷材料、机敏材料、智能材料、生物材料以及纳米材料等。

（2）为了分析和模拟各种类型和形式的结构在复杂载荷工况和环境作用下的全寿命过程的响应，需要发展新的数值分析方案。如：多重非线性（材料、几何、边界等）相耦合的分析方法；多场（结构、流体、热、电、化学）耦合作用的分析方法；跨时间／空间多尺度；非确定性（随机／模糊）的分析方法；自适应的分析方法。

（3）有限元软件和CAD／CAM等软件系统共同集成完整的虚拟产品开发（VPD）系统。这个系统强烈影响着未来工程系统的设计、制造、和运行，主要体现在：

①它能提供对所设计的工程系统从加工制造到运行，直至失效和破坏的全寿命过程的更深入认识，从而能更好地识别它的属性和特征。

②它能够鉴定和评估所设计对象的性能和质量，并允许以最低的费用在设计过程中就对所设计的对象进行修改和优化。

③它能显著地缩短工程对象设计和投产的周期，降低生产成本，提高市场竞争力。

2015年3月17日作业

1、简述有限元法的基本思想，并结合简单结构来说明。

(1)有限元法，也叫有限单元法，它的基本思想是将一个结构或连续体的求解域离散为若干个子域（单元），并通过它们边界上的结点相互联结成为组合体。

(2)有限元法用每一个单元内所假设的近似函数来分片地表示全求解域内待求的未知场变量。而每个单元内的近似函数由未知函数或其导数在单元各个结点上的数值和与其对应的插值函数来表示。由于在联结相邻单元的结点上，场函数应具有相同的数值，因而将它们用作数值求解的基本未知量。这样一来，求解原

来待求场函数的无穷自由度问题转换为求解场函数结点值的有限自由度问题。

(3)有限元法是通过和原问题数学模型（基本方程、边界条件）等效的变分原理或加权余量法，建立求解基本未知量（场函数的结点值）的代数方程组或微分方程组。此方程组称为有限元求解方程，并表示成规范的矩阵形式。接着用数值方法求解此方程，从而得到问题的解答。

2、简述结构离散(或有限元建模)的内容和要求。

有限元建模的内容：

1）网格划分---即把结构按一定规则分割成有限单元

2）边界处理---即把作用于结构边界上约束和载荷处理为结点约束和结点载荷

有限元建模的要求：

1）离散结构必须与原始结构保形---单元的几何特性

2）一个单元内的物理特性必须相同---单元的物理特性

3、简述结点力和结点载荷的差别。

结点力：单元与单元间通过结点的相互作用力。

结点载荷：作用于结点上的外载荷。

4、列表给出有限元几类基本单元的图形、结点数、结点自由度数和单元总自由度数(包括杆单元、梁单元、平面三角形单元、平面四边形单元、轴对称问题三角形单元、四边形壳单元、四面体单元)。

5、写出适用于插值下表的基函数，并给出拉格朗日插值多项式。

基函数：2

92451813

11413141

)31(2-+-=--?--

=x x x x l k ，31636416)41(2+-=x x l k

，6

1

672)1(2+-=x x l k

根据拉格朗日插值多项式：)()()()()()(000k n k n k

j j j

k j

n k n k k n x f x x x x x p x f x l x p ∑∏∑=≠==--==或。

将2=n 带入：3

35

634938-)(22-+=x x x p

2015年3月24日作业

1、为了保证有限单元法解答的收敛性，位移函数应满足哪些条件？完备协调元、非协调元和完备元分别是什么意思？

为了保证有限单元法解答的收敛性，位移函数应满足： 1) 位移函数必须包括单元的刚性位移（即常量项）； 2) 位移函数必须包括常量应变（即线性项）； 3) 位移函数在单元内部必须连续（连续性条件）； 4) 位移函数应使得相邻单元间的位移协调（协调性条件）

注：上述四个条件称为有限元解收敛于真实解的充分条件；前三个条件称为必要条件。满足四个条件的位移函数构成的单元称为完备协调元；满足前三个条件的单元称为非协调元；满足前两个条件的单元称为完备元。

2、如下图所示4结点平面应力单元，结点1~结点4对应的结点坐标分别为(0,0), (0,1), (2,0),(2,1)，结点1~结点4对应的结点位移分别为(u1,v1), (u2,v2), (u3,v3), (u4,v4)，试基于拉格朗日插值基函数构造如下单元的位移函数。

解：根据拉格朗日插值基函数：

44332211),(),(),(),(),(u y x l u y x l u y x l u y x l y x u +++=

xy

y x l y x y x l y

x y x l y x y x y x l 21

),()

1(21

),()2(21

),()

1)(2(21

)10)(20()1)(2(),(4321=--=--=--=----=

()

1

0,0()

20,1()32,0()

42,1

所以：432121

)1(21)2(21)1)(2(21),(xyu u y x yu x u y x y x u +------=

同求：),(y x v 。

3、说明有限元方法解误差的主要来源？

答：影响有限元解的误差：1) 离散误差。边界上以直线代曲线导致离散化模型与实际物体的差异；2) 位移函数误差。一般情况下单元位移函数不可能与实际单元的位移场一致；3) 计算机计算误差。计算机字长的限制、相差悬殊的数值加减运算。

4、说明用有限单元法解题的主要步骤。

答：研究问题的力学建模；结构离散；单元分析；整体分析与求解；结果分析及后处理。

5、推导基于变分原理的总势能泛函极值条件。 解：有积分形式确立的标量泛函有

,,d E ,,d u

u

ΠF u Ωu x x Ω

Γ

??????

=+Γ ? ???????

??L

L 其中u 是未知函数，F 和E 是特定的算子，Ω是求解域，Γ是Ω的边界。Π称为未知函数u 的泛函，随函数u 的变化而变化。连续介质问题的解u 使泛函Π对于微小的变化u δ取驻值，即泛函的“变分”等于零0Πδ=，此为变分法。

将虚功原理用于弹性变形时，总功W 要包括外力功(T)和内力功(U)两部分，即：W = T - U ；内力功(-U)前面有一负号，是由于弹性体在变形过程中，内力是克服变形而产生的，所有内力的方向总是与变形的方向相反，所以内力功取负值。

根据虚功原理，总功等于零得：T- U = 0，即外力虚功T = 内力虚功U 弹性力学中的虚功原理可表达为：在外力作用下处于平衡状态的弹性体，如果发生了虚位移，那么所有的外力在虚位移上的虚功(外力功)等于整个弹性体内应力在虚应变上的虚功(内力功)。

根据虚功原理得到()0T T T d d δδδΩ

Γ

-Ω-Γ=??εu F u T σ

1

()02T T T p d d δδΩΓ

??∏=-Ω-Γ= ??

?

??u F u T εσ

其中的p ∏即为总势能泛函。由上面变分为零式表明：在所有区域内满足几何关系，在边界上满足给定位移条件的可能位移中，真实位移使系统的总势能取驻值(可证明此驻值为最小值)。此即总势能泛函的极值条件。

2015年3月31日作业

1、给出利用最小势能原理建立单元有限元静力平衡方程的一般推导过程。 解：由单元位移函数：

e u =Nu

式中：N 为插值函数（或称形函数），得到单元内的应变和应力分别为

e ε=Bu e e σ=D ε=DBu =Su

其中D 为弹性矩阵，它完全取决于弹性常数E 和μ。将位移、应力和应变代入势能泛函有

1(d )(d )(d )2

eT

T e eT T eT T u ΩΩΓ∏=

Ω-Ω-Γ???u B DB u N F u N T

根据最小势能原理，势能泛函取驻值的必要条件：

(d )(d )(d )0T e T T eT ΩΩΓ?∏=Ω-Ω-Γ=????B DB u N F N T u

式中d T e Ω

Ω=?B DB K 称为单元刚度矩阵，d T e V Ω

Ω=?N F f 称为体积力等效

结点力，d T e S Γ

Γ=?N T f 称为面力等效结点力，则上式可写成单元平衡方程：

e e e =K u f

此式即单元有限元静力平衡方程。

2、利用最小势能原理推导结构（假设包含N 个单元）有限元静力平衡方程。

解：利用弹性体所有单元的总势能最小，将所有单元进行组装后即可形成整个结构的有限元静力平衡方程。为此先需要建立单元结点位移向量和总体结点位移向量之间的对应关系，即

i i u =T u

式中：i u 表示第i 个单元的结点位移向量；T u 是总体结点位移向量；矩阵i T 是第

i 个单元的结点位移提取矩阵。对于杆单元可以表示为

12[,]T i i

i u u =u 121[,,,]T n u u u +=L u 0100[

]0010

i =L

L L L T 利用提取矩阵对单元矩阵进行合同变换与累加后，可得到总体的刚度矩阵

K 、质量矩阵M 和载荷向量F ，即

1n T

i i i i ==∑K T k T 1n T

i i i i ==∑M T m T 1

n

T i i i ==∑F T f

根据最小势能变分原理即可得到Ku =F ，此式即为总体结构的有限元静力平衡方程。

2015年4月7日作业

1、推导拉(压)杆单元形状函数，利用最小势能原理推导拉(压)杆单元在局部坐标系中的单元刚度方程。

解：（1）如图1表示某一杆单元，现约定附属于该单元的局部坐标系（单元坐标系）为''''o x y z ，i 点为原点，x 轴沿着杆轴线，其正方向为由i 指向j ，其余各轴按右手螺旋规则确定。设i u ，i v ，i w ，j u ，j v ，j w 为杆元结点位移分量xi F ，yi F ，zi F ，xj F ，yj F ，zj F 杆单元结点力分量，一律规定和坐标轴正向一致时为正。设杆的长度为l ，弹性模量为E ，横截面积为A 。

图1 杆单元结点位移结点力分量

对于铰接杆单元，在小变形假设的前提下，与杆垂直方向的位移并不使杆产生应变和应力。于是，对每一个结点只需考虑一个结点位移和结点力，即如图2所示

图2 二力杆单元

二力杆单元的位移函数为12()u x a a x =+。式中1a ，2a 是两个待定常数，可由i ，j 两结点的位移唯一确定。当()0,0i x u u ==；,()j x l u l u ==时，将其代

入位移函数有：1i u a =，12j u a a l =+，从而可以得到

12i j i a u u u a l

=-=

将1a ，2a 的值代入位移函数得

()1j i

i i j u u x x u x u x u u l

l l -??

=+

=-+ ???

或写成

()1i i i j j j u u x

x u x N N u u

l

l ??

??????=-

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

e

Nu 则有

i i j j u N u N u ==+e Nu

1i j x x

N N l l

=-=，

式中i N ，j N 称为i 点、j 点的形状函数，N 成为形函数矩阵。

（2）根据位移函数与应变的定义x du

dx

ε=

得 ()[]1111x d d x d x dx dx l dx l l εε??????==

=-=- ? ?????????

e e e

Nu u u 或写成=e εBu ，其中[]1/11l =-B 称为应变矩阵。

根据应力与应变的关系=σD ε，将应变代入有=e e σD ε=DBu =Su ，其中

S =DB 称为应力矩阵，对于拉（压）杆有[]11E

l

=

-S 。 根据杆单元在局部坐标系单元刚度矩阵的表达式dV Ω

=?e T K B DB 代入矩阵

[]1/11l =-B 与[]11E

l

=-S =DB 得到拉(压)杆单元在局部坐标系中的单元刚度方程为

1111EA l -??

=

??-??

e K

2、给出将局部坐标系中的拉(压)杆单元刚度矩阵通过坐标变换得到总体坐标系中的单元刚度矩阵推导过程。

解：设oxyz 为总体坐标系，''''o x y z 为局部坐标系，如图1所示。

图1 平面杆单元总体坐标系位移在局部轴方向分量

规定由总体坐标系平面x 轴到局部坐标系'x 轴的夹角α逆时针为正。杆单元总体坐标系下的结点位移分量用,u v 表示，局部坐标下的位移分量用','u v 表示。则平面杆单元结点i 在总体坐标系和局部坐标系下的位移分量关系有

''

cos sin sin cos i i i i

i i u u v v u v αααα

=+=-+

同理对于j 结点有

''cos sin sin cos j j j j j j u u v v u v αααα

=+=-+

用矩阵表示为

()''cos sin ,sin cos l l l l u u l i j v v αααα????

??'====?????

?-????

??l l

u λu

对两结点杆单元，当用总体坐标系位移e u 表示局部坐标系中位移()e

u'时有转换关系

()

''''i i i i j j j j u u v v u u v v ??????????????===??????????????????

??

e

e λu'Tu λ 式中T 矩阵称为坐标变换矩阵。对于图1中的杆单元就可以很容易将局部坐标系的刚度矩阵转换为总体坐标系的刚度矩阵为

cos sin 001

010sin cos 00000000cos sin 1

01000sin cos 0000cos sin 00sin cos 0000cos sin 0

sin cos a EA l a α

ααααααααααααα--?????????

??

?=????

--???

?

????

????-????

??

-??

e

K

22222222cos cos sin cos cos sin cos sin sin cos sin sin cos cos sin cos cos sin cos sin sin cos sin sin a a a a EA l a a a

a αααααααααααααα

αα

??

--?

?

--??

=

??--?

?

--??

3、计算下图所示三杆平面桁架结构的总体坐标系下各单元的单元刚度矩阵。杆单元材料弹性模量E=2.0×10^11(N/m^2)，截面积A=1.0×10^(-4)(m^2)。注：物理量单位均采用国际单位。

解：总体坐标系下的单元刚度矩阵为

22222222cos cos sin cos cos sin cos sin sin cos sin sin cos cos sin cos cos sin cos sin sin cos sin sin a a a a EA l a a a

a αααααααααααααα

αα

??

--??

--??

=

??

--?

?

--??

e K

对杆①有：10α= 1sin 0α= 1cos 1α=

7111410101

01000000000510101010100.4

00002.010 1.0000100e

---????

?

??????????=

=?????--???

?

????

K 对杆②有：290α= 2sin 1α= 2cos 0α=

1172400001010010100006.6710

000010100.3

02.010 1.1010000001e --????

???

?-????==?????

-????

-????

???K 对杆③有：33arctan 4α=- 33sin -5α= 34

cos 5

α=

4311716121612161216122525252525252525129129129129252525252525252541016121612161216120.5252525252525252512912912912925252522.01525

2525250 1.010e -????----?????????----???==???????----?????????----

??????K ?

?????????

4、推导平面梁弯曲的挠曲线控制方程。

解：如下图为梁单元，结点为i 和j

因弯曲而引起的轴向位移为d (,)d v u x y y x ??

=- ???

，其中y 是所讨论点离中面的

距离，应变22d d d d x u v y x x ε==-。又根据平衡关系有x yM I

σ=，其中I 是切面惯性矩，2d A

I y A =?。根据胡克定律，x x E σε=有

22d d yM

v E y I x ??=- ???

22d d v M EI x

=-

又根据剪力与弯矩的关系d d M

Q x

=

有 2222

22d d d d d d M v q EI x x x ??

=-=- ???

现在讨论的梁E ，I 为常数，

故有44d d v EI q x =，若梁上无分布剪力，则有44d 0d v

EI x

=，则可以判断v 是x 的三次函数。设

230123()v x b b x b x b x =+++

则梁的转角

2123d ()23d v

x b b x b x x

θ=

=++ 将i ，j 结点位移代入上式有

01

2

3

2

0123123, , 23i i j j v b b v b b l b l b l b b l b l

θθ===+++=++

用矩阵表示为

01232231000010000

123i i j j v b b v l l l b l l b θθ????

??????

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

??

????????

?? 从而得

0122233

2321000010032312

121i i j j v b b v b l l l l b l l l l θθ??

??????????????????=????---??????

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

-????

232

23

2

32100001003231()1,,,2121i i i i j j j j v v v x x x x v v l l l l l l l

l θθθθ??

????????????????????==????---????????

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

v N

单元的形状函数矩阵v N 为

232323232322320101322321,,, i i j j x x x x x x x x x l l l l l l l l H H H H ??

=-+-+--+??

?

???=??

v N 其中

2323

01232

2323

01232

3221, 32, i i j j x x x x H H x l l l l

x x x x H H l l l l

=-+=-+=-=-+ 设x

l

ξ=

，则上式可写成 ()()

2323012

3

3

2

01132, 232, i i j j H H l a

H H l ξξξξξξξξξ

=-+=-+=-=-

以上四个函数即是二结点梁单元的形状函数．由于每结点有两个位移参数，则每结点有两个形状函数。0i H 是指i 结点零阶导数（即i 点位移）对应的形函数，

1i H 指i 结点一阶导数（即i 结点转角）对应的形函数，位移函数用内插多项式表

示为

0101()i i i i i i j j j j j j v v x H v H H v H v θθθθ??

????

=+++=????????

v N

2015年4月14日作业

1、推导平面直梁单元形状函数，利用最小势能原理推导平面直梁单元在局部坐标系中的单元刚度方程。

解：先将轴力与剪力、弯矩分开考虑。对于平面直梁单元的两个节点1,2各有一个挠度和转角，故可设其挠度为：

23120123112324()()()

()()T w dw dw

w x b b x b x b x N x w N x N x w N x dx dx

=+++=+++=N w 其中，(0,1,2,3)i b i =为待定常数，()(1,2,3,4)i N x i =为单元形函数，

1234[()()()

()]T w N x N x N x N x =N ；

12

3[]T b b b b =w

利用Hermit 插值或书中解方程的方法，可求得：

231()132x x N x L L ????=-+ ? ????? ，23

2()2x x N x x L L L L ????

=-+ ? ?????

2

3

3()32x x N x L L ????=- ? ?????，2

3

4()x x N x L L L L ????

=-+ ? ?

???? 对于轴力引起的位移()u x ，设：

01()i j u x a a x u u =+=????

?

?????

u N 1,x x l l ??

=-????

u N

进而，位移函数改写为

11121

2

3

422()10000()0

e u w x

x

u x l

l u v x N N N N w θθ??????????-????

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

??????

??????

Nu 位移函数求得后，可得到应变和应力的表达式：若忽略剪切影响，N ε是拉力引起的应变，b ε是弯曲引起的应变

1212212222d d d d d d d d N b w u u x u w x w y y x w x εεθθ??

?????????

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

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

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

?????????

u N N ε