岩土工程数值分析试卷及参考答案

岩土工程数值分析试题
一、简答题（40分）
1．简述梁单元、杆单元、连续梁单元、平面三角形常量单元和四边形等参单元的特点（10分）。

答：1）梁单元是由两个节点组成，每一个节点都具有三个方向的线性移动位移和三个方向的旋转位移，因而每个节点具有6个自由度，梁单元具有拉，压，剪，弯，扭的变形刚度。

计算理论成熟，建模方便,计算量小，在工程结构有限元分析中得到广泛的应用，适用于各种截面形式的杆件分析。

2）由有限个构件以一定方式连接起来所形成的结构，在同一平面内的杆系结构，其所受的外力作用线位于该平面内，在杆系中，每一个杆件可视为一个单元，每个单元的端点成为结点。

3）对于每跨各自等截面的连续梁，以每跨为一个单元。

结点编号和单元编号一般是从连续梁的左端顺序编到右端。

由于连续梁各单元的轴线方向一致，各单元坐标系与结构坐标系的方向相同，因此在矩阵位移法的计算过程中无须进行坐标变换，在单元坐标系和结构坐标系中单元刚度矩阵的表达式是相同的。

4) 平面三角形单元具有适应性强的优点，较容易进行网络划分和逼近边界形状，应用
比较灵活。

其缺点是它的位移模式是线性函数，单元应力和应变都是常数，精度不够理想。

5) 四边形等参单元能更好地反映物体内的应力变化，适应曲线边界，常使用于弹性力
学平平面问题的分析。

八结点单元一共有16个已知的结点位移分量。

2．除有限单元法外，岩土工程常用到哪些数值方法，并对比其优缺点(10分)。

答：岩土工程常用的数值方法包括：有限差分法、边界元法、离散元法、颗粒元法、不连续变形分析法、流形元法、模糊数学方法、概率论与可靠度分析方法、灰色系统理论、人工智能与专家系统、神经网络方法、时间序列分析法。

有限单元法的优缺点：有限单元法的理论基础是虚功原理和基于最小势能的变分原理，它将研究域离散化，对位移场和应力场的连续性进行物理近似。

有限单元法适用性广泛，从理论上讲对任何问题都适用，但计算速度相对较慢。

即，物理概念清晰、灵活、通用、计算速度叫慢。

有限差分法：该方法适合求解非线性大变形问题，在岩土力学计算中有广泛的应用。

有限差分法和有限单元法都产生一组待解方程组。

尽管这些方程是通过不同方式推导出来的，但两者产生的方程是一致。

另外，有限单元程序通常要将单元矩阵组合成大型整体刚度矩阵，而有限差分则无需如此，因为它相对高效地在每个计算步重新生成有限差分方程。

在有限单元法中，常采用隐式、矩阵解算方法，而有限差分法则通常采用“显式”、时间递步法解算代数方程。

边界元法：该方法的理论基础是功互等定理和基本解，它只要离散求解域的边界，因而得到离散代数方程组中的未知量也只是边界上的量。

边界元法化微分方程为边界积分方程，离散划分少，可以考虑远场应力，有降低维数的优点，可以用较少的内存解决较大的问题，便于提高计算速度。

离散元法：离散元法的理论基础是牛顿第二定律并结合不同的本构关系，适用对非连续体如岩体问题求解。

该方法利用岩体的断裂面进行网格划分，每个单元就是被断裂面切割的岩块，视岩块的运动主要受控于岩体节理系统。

它采用显式求解的方法，按照块体运动、弱面产生变形，变形是接触区的滑动和转动，由牛顿定律、运动学方程求解，无需形
成大型矩阵而直接按时步迭代求解，在求解过程中允许块体间开裂、错动，并可以脱离母体而下落。

离散元法对破碎岩石工程，动态和准动态问题能给出较好解答。

颗粒元法：颗粒元方法是通过离散单元方法来模拟圆形颗粒介质的运动及其相互作用，它采用数值方法将物体分为有代表性的多个颗粒单元，通过颗粒间的相互作用来表达整个宏观物体的应力响应，从而利用局部的模拟结果来计算颗粒群群体的运动与应力场特征。

不连续变形分析方法：该方法是并行于有限单元法的一种方法，其不同之处是可以计算不连续面的错位、滑移、开裂和旋转等大位移的静力和动力问题。

此方法在岩石力学中的应用备受关注。

流形元法；该方法是运用现代数学“流形”的有限覆盖技术所建立起来的一种新的数值方法。

有限覆盖是由物理覆盖和数学覆盖所组成的，它可以处理连续和非连续的问题，在统一解决有限单元法、不连续变形分析法和其他数值方法的耦合计算方面，有重要的应用前景。

无单元法：该方法是一种不划分单元的数值计算方法，它采用滑动最小二乘法所产生的光滑函数去近似场函数，而且又保留了有限单元法的一些特点。

它只要求结点处的信息，而不需要也没有单元的信息。

无单元法可以求解具有复杂边界条件的边值问题，如开裂问题，只要加密离散点就可以跟踪裂缝的传播。

它在解决岩石力学非线性、非连续问题等方面具有重要价值和发展前景。

混合法：对于复杂工程问题，可采用混合法，即有限单元法、边界元法、离散元法等两两耦合来求解。

模糊数学方法：模糊理论用隶属函数代替确定论中的特征函数描述边界不清的过渡性问题，模糊模式识别和综合评判理论对多因素问题分析适用。

概率论与可靠度分析方法：运用概率论方法分析事件发生的概率，进行安全和可靠度评价。

对岩土力学而言，包括岩石（土）的稳定性判断、强度预测预报、工程可靠度分析、顶板稳定性分析、地震研究、基础工程稳定性研究等。

灰色系统理论：以“灰色、灰关系、灰数”为特征，研究介于“黑色”和“白色”之间事件的特征，在社会科学及自然科学领域应用广泛。

岩土力学中，用灰色系统理论进行岩体分类、滑坡发生时间预测、岩爆分析与预测、基础工程稳定性、工程结构分析，用灰色关联度分析岩土体稳定性因素主次关系等。

人工智能与专家系统：应用专家的知识进行知识处理、知识运用、搜索、不确定性推理分析复杂问题并给出合理的建议和决策。

岩石力学中，可进行如岩土（石）分类、稳定性分析、支护设计、加固方案优化等研究。

神经网络方法：试图模拟人脑神经系统的组织方式来构成新型的信息处理系统，通过神经网络的学习、记忆和推理过程进行信息处理。

岩石力学中，用于各种岩土力学参数分析、地应力处理、地压预测、岩土分类、稳定性评价与预测等。

时间序列分析法：通过对系统行为的涨落规律统计，用时间序列函数研究系统的动态力学行为。

岩石力学中，用于矿压显现规律研究、岩石蠕变、岩石工程的位移、边坡和硐室稳定性等、基础工程中降水、开挖、沉降变形等与时间相关的问题。

3．简述商用分析软件的基本框架（10分）。

答：运行程序后，整个构架分为两层，一个是起始层，另一个是处理层。

这两层的关系主要是使用命令输入时，要通过起始层进入不同的处理器。

处理器可视为解决问题步骤中的组合命令。

1.起始层
当进入程序后，所在的位置便是起始层。

可以在输入窗口中看到所在的位置。

起始层的功能是：
（1）进入不同的处理器层。

（2）清除当前的所有资料，开始新任务。

（3）改变工作文件名称。

（4）退出。

2.处理层
处理器包含进行分析的操作命令。

的常用处理有： 1）前处理
（1）建立有限元模型所需输入的资料，如单元节点坐标资料、单元节点排列次序。

（2）定义材料属性。

（3）单元切割。

2）求解处理
（1）定义负载条件。

（2）定义所施加的外力。

（3）定义边界条件及进行求解。

3）后处理
（1）从求解计算结果中读取数据。

（2）对计算结果进行各种图形化显示。

（3）对计算结果进行列表显示。

（4）进行各种后续分析。

4．简述平面杆系结构有限元法的解题步骤（10分）。

答：平面杆系结构有限元法的解题步骤：
（1）结构离散化。

对于杆系结构，通常以一根等截面直杆作为一个单元，把整个结构视为由有限个单元组成的集合体。

这相当于建立位移法的基本结构。

（2）单元分析。

对于杆件单元，只要确定了杆件两端截面的内力，其余截面内力便容易确定，因而杆件的内力可用杆端力作代表。

单元分析的任务是确定杆端力和杆端位移之间的关系，建立所谓单元刚度矩阵。

这与位移法中杆件的形常数及转角位移方程相对应。

（3）整体分析。

把各单元集合成整体结构，这就要求在各结点处满足变形协调条件和平衡条件。

整体分析的任务是由单元刚度矩阵按照刚度集成规则直接形成结构刚度矩阵，并建立整体结构的刚度方程，它反映了结构的结点位移与结点荷载之间的关系，可由此方程求解结点位移。

求出了结点位移，再回到单元分析，便可确定各单元的内力。

二、计算题（60分）
1．试验证：在三节点常量单元内任意一点，有
1=++m j i N N N ，x x N x N x N m m j j i i =++，y y N y N y N m m j j i i =++（10分）
解：由于该三节点为常量单元，所以满足线性位移模式，即：
⎩
⎨
⎧++=++=y x y x y y
x y x x 654321),(),(αααααα （1） 式（1）可写为
α0N y x x =⎭⎬⎫
⎩⎨⎧=
其中
]0
0[
0p
p N N N =
]
1
[y x
N p =
{}T 654321ααααααα=
将三角形单元的结点坐标代入式（1），得
αN x e = （2）
其中
{}T
m m
j
j
i i
e y x y x y x x =
⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢
⎢⎢
⎢⎣⎡=m m
m m j j j j
i i i i
y x y x y x y x y x y x N 10
000
11000000110000001 {}T 654321ααααααα=
从式（2）可以求出α：
e x N 1-=α
故得
⎪⎪⎪⎪⎭⎪
⎪⎪⎪⎬⎫
⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∆=
⎪⎪
⎪
⎪
⎭⎪⎪⎪
⎪
⎬⎫
⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧m m j j i i m j
i
m j i m j i m j i
m j i m j i y x y x y x c c c b b b a a a c c c b b b a a a 0
00000000000000021654321αααααα （3）
其中
⎪⎩⎪
⎨⎧-=-=-=-=-=-=-=-=-=i
j m j
i m i j j i m
m i j i m j m i i m j j
m i m
j i j m m j i x x c y y b y
x y x a x x c y y b y x y x a x x c y y b y x y x a （4）
-=∆m
m
j j
i
i y x y x y x 111
2三角形面积 （5） 将式（3）~式（5）代入式（1），得
⎩⎨
⎧=++=++y y N y N y N x
x N x N x N m m j j i
i m m j j i i （6） 其中
∆++=
2),(y
c x b a y x N i i i i
∆
++=
2),(y
c x b a y x N j j j j
∆
++=
2),(y
c x b a y x N m m m m
式（6）可写成
{}[]{}
e
x N y x y y x x x =⎭
⎬⎫⎩⎨⎧=),(),(
其中
⎥⎦
⎤
⎢
⎣⎡=m j
i
m j i
N N N N N N N 0
000][ {}T
m m
j j i i
e y x y x y x x =][
式中：[N]为形函数矩阵；[]为单元结点位移列阵。

形函数有以下两个特点：
1）),,)(,(m j i i y x N i i i =在结点i 处为1，其他结点处为0，即
m
j i i y x N i i i ,,1),(==
i j m j i j y x N j j i ≠==;,,0
),(
2）全部形函数之和等于1，即
m j i i y x N i
,,1
),(==∑
根据这个性质，不管i x 、j x 、m x 为何值，由式（6），均有
m
j i i x y x x i
i i ,,),(=≡
设想单元沿x 轴平移x 0，则将 x 0代入式（6），得到
()00000x N N N x N x N x N x x m j i m j i ++=++==
所以证得
1
=++m j i N N N
2．试用最小势能原理，推导有限元法基本公式（15分）。

解：设单元的节点位移列向量为e
δ，那么单元内任一点（）的位移列向量d 可以利用节点位移值的插值函数来表示，即e
N d δ=。

根据几何方程，单元内的应变也可以用节点位移来表示，即
e B d L δε==
其中矩阵B 称为应变矩阵或几何矩阵，矩阵B 中的元素可能是常数，也可能是的函数。

根据线弹性的广义胡克定律，单元内的应力为
e e S DB D δδεσ===
其中矩阵D 称为弹性常数矩阵，矩阵DB S =称为应力矩阵，S 中的元素可能是常数，也可能是的函数。

单元之间通过节点相互作用的力即为单元节点力列向量e
F 设单元内体积力荷载列向量为
T
bz by bx b F F F F ),,(=
设除去单元节点力e
F 外的表面力荷载列向量为
T
sz sy sx s F F F F ),,(=
在单元上利用最小总势能原理得
]])([[)(0=---=-=⎰⎰⎰S
e T
e T s V
T b F dS d F dV d F d U W U V δδδδ
即
)(][]
21[]
21[=---=---=---=⎰⎰⎰⎰⎰⎰⎰⎰
⎰⎰⎰⎰
eT
S
e s T V
b T e
V
T S
e eT s T eT V b T eT e V
T
eT S
e eT s T V b T V
T
F dS F N dV F N dV DB B F dS F N dV F N dV DB B F dS F d dV F d dV V δ
δδδδδδδδδσεδδ
引入单元刚度矩阵为
dV
DB B k V
T e
⎰=
引入单元等效节点荷载矩阵（列向量）为
⎰⎰⎰+=S
s T V
b T E
e dS
F N dV F N F
在引入单元刚度矩阵和单元等效节点荷载矩载矩阵之后，由于)(eT
δδ的任意性，则可得到单元的平衡方程即单元刚度方程的一般形式为
e
e E
e e k F
F δ=+
3．用有限单元法中的先处理法，计算如图1所示结构。

已知：3×1072 ；梁式杆的0.20m 2,
0.004m 4；桁杆的0.06m 2。

注意：只需完成下列步骤要求（15分）
（1）建立坐标系，离散结构（包括单元编号、节点编号、单元信息和节点坐标）； （2）写出每个单元的定位向量；
（3）分别写出一个杆单元和一个梁单元的刚度矩阵（整体坐标系下）。

（4）写出由单刚集成总刚的方法与步骤； （5）形成总体荷载列阵。

第二大题第一小题图
15kN/m
1m
1m
2m
45
q
图1
解：（1）建立如图所示直角坐标系，结构离散后，单元编号及结点编号如图所示。

（2）各单元的定位向量如下：
{}{}4,3,2,1,0,0)1(=λ，{}{}22,21,20,5,3,2)2(=λ，{}{}9,8,7,6,3,2)3(=λ {}{}17,16,15,10,8,7)4(=λ，{}{}13,0,12,11,8,7)5(=λ，{}{}18,16,15,14,0,12)6(=λ {}{}23,21,20,19,16,15)7(=λ，{}{}24,21,20,1,0,0)8(=λ
（3）在整体坐标下杆单元②的刚度矩阵为：
{}⎪⎪⎭⎪
⎪⎬⎫
⎪⎪⎩⎪⎪⎨
⎧--⨯⨯=101
00000101
000
01065m k
在整体坐标下梁单元⑦的刚度矩阵为：
[]⎥⎥
⎥
⎥⎥⎥
⎥
⎥
⎥
⎦⎤⎢⎢⎢⎢⎢⎢
⎢⎢⎢⎣⎡⨯⨯-⨯⨯⨯-⨯⨯-⨯--⨯⨯-⨯⨯⨯⨯-⨯⨯=44
4
4
4
4446644
4
4
4444661081020
1041020
102103
2
01021032
00010001010410201081020
102103
2
0102103
2
00
100010k （4）编码法，其步骤为：1）将总体刚度矩阵中的全部元素置零0，1,2n ，1,2n 。

2）形成每个单元中的刚度矩阵，并寻求单元的结点编号i 、j 。

3）将单元刚度矩阵中的各元素(1~41~4)分别按局部自由度与总体自由度编号的关系
IJ rs
K J I j j i i s r k →→⎪⎪⎭
⎪⎪⎬⎫
⎪⎪⎩⎪⎪⎨⎧--→→2412322121
组集到总刚矩阵的位置上。

4）依次完成个单元的组集，形成总体刚度矩阵[K]。

必要时，将(1)累加到对应的中。

（5）如图所示，总体荷载矩阵为
{}6
543
21090045450060000000045750⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎪
⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬
⎫⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧---=P
4．用弹性平面问题有限单元法（后处理法）计算图2所示平面问题。

已知：F 1 F 1y =10; F 220；1x107 2，泊松比μ=0.3，γ=243，薄板厚1单位，图中水平与竖向单元边长均为1m 。

完成下列步骤要求（20分）。

（1）建立整体坐标系；
（2）计算单元①，②单元刚度矩阵； （3）计算结点荷载列向量；
（4）按照赵更新主编教材P223程序的输入模式，写出文件的数据。

图2
解：（1）建立如图所示整体直角坐标系。

（2）单元①的刚度矩阵为：
[])
1(333231232221131211)
1(][][][][][][][][][⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡=K K K K K K K K K K
由单元结点坐标得：1/2，b 1=0，b 21，b 3=1，c 1=1，c 21，c 3=0。

由公式
)
,,;,,(21212
121)1(4][2m j i s m j i r b b c c c b b c b c c b c c b b A Et K s r s r s
r s r s r s r s r s r rs ==⎥
⎥⎥⎦
⎤
⎢⎢⎢⎣
⎡-+-+-+
-+-=μ
μμμμμμ
可得⎥⎦⎤⎢⎣⎡⨯⨯=10035.010915][811K ，⎥⎦⎤⎢⎣⎡----⨯⨯=13.035.035.010915][8
12K ⎥⎦⎤⎢⎣⎡⨯⨯=
03.035.0010915][813K ，⎥⎦
⎤⎢⎣⎡----⨯⨯=135.03.035.010915][8
21K
⎥⎦⎤⎢⎣⎡⨯⨯=35.165.065.035.110915][822K ，⎥⎦⎤⎢⎣⎡----⨯⨯=35.03.035.0110915][8
23K
⎥⎦⎤⎢⎣⎡⨯⨯=
035.03.0010915][831K ，⎥⎦
⎤⎢⎣⎡----⨯⨯=35.035.03.0110915][8
32K
⎥⎦
⎤⎢⎣⎡⨯⨯=
35.000110915
][833K 所以综合得：
[])
1(8
)
1(35.0035.035.0035.0013.013.0035.03.035.165.0135.035.0165.035.13.035.003.013.01035.0035.035.0035
.010915⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------⨯⨯=K
同理可得单元②的刚度矩阵为：
[])
2(8
)
2(35.0035.035.0035.0013.013.0035.03.035.165.0135.035.0165.035.13.035.003.013.01035.0035.035.0035.010915⎥⎥⎥⎥⎥
⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-------------------⨯⨯=K
（3）各单元自重的等效结点荷载计算如下：
kN tA 42
1
1243131=⨯⨯⨯=γ 321404040}{)1(⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧---=E P ，542404040}{)2(⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧---=E P ，532404040}{)3(⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧---=E P ，6
53404040}{)4(⎪⎪⎪⎪⎭⎪
⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧---=E P
65432
1401204012012040654321404440404440444040}{⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎪⎪⎪⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧------=⎪⎪⎪⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧------------=E P 形成结构的综合结点荷载列阵{P}。

654321402805430120122061065432140120401201204065432100400503000020101065432140120401201204065432100000201010}{}{}{6544⎪⎪⎪⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧-----=⎪⎪⎪⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧------+⎪⎪⎪⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧--=⎪⎪⎪⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧------+⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎪⎪⎪⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧=+=y y y x E J R R R R P P P （4）4，6，3，3，0，0，
1.0E7，0.3，1，24
0.0，0.0，1.0，0.0，1.0，2.0
2.0，1.0，1.0，0.0，0.0，0.0
1，2，3
2，4，5
2，5，3
3，5，6
4，1，1
5，0，1
6，0，1
1，1，10
1，2，10
2，1，20。