计算固体力学课程作业
专 业 固 体 力 学 学 号 1131301009 姓 名 尹亚川
作业1：
（一）、H ϕ+f =0，其中f =10，H =10(1+e 8ϕ)
（1） 试用直接迭代法，Newton-Raphson 方法，修正Newton-Raphson 方法，拟
Newton-Raphson 方法进行求解并进行比较。 （2） 用Euler-Newton 法计算，f 分2级
求解：
（1）直接迭代法：
于是得近似解
H ϕ+f =0
（1） （2）
H 0=H (ϕ0)
ϕ1=(H 0) -1(-f )
重复这一过程，以第i 次近似解求出第i +1次近似解的迭代公式为
（3）
直到
H i =H (ϕi )
（4） （5）
ϕi +1=(H i ) -1(-f )
∆ϕ=ϕi +1-ϕi
（6）
变得充分小，即近似解收敛时，终止迭代。
取ϕ0=1，令∆ϕ
取ϕ0=0，令∆ϕ
取ϕ0=-1，令∆ϕ
H ϕ+f =0
（7） （8） （9）
ψi =ψ(ϕi ) ≡F (ϕi ) -R =H i ϕi +f ≠0
⎛∂ψi i
K T =K T (ϕi ) ≡ ∂ϕ
⎝
⎫⎪⎪ ⎭
i
i -1i i -1
∆ϕi =-(K T ) ψ=(K T ) (H i ϕi +f )
（10） （11） （12）
⎛∂ψi
K T = ∂ϕ
⎝⎫⎛∂H ϕ⎫⎪⎪= ∂ϕ⎪⎪ ⎭⎝⎭
i i
ϕi +1=ϕi +∆ϕi
当∆ϕi 变得充分小，即近似解收敛时，终止迭代。
取ϕ0=1，令∆ϕ
取ϕ0=0，令∆ϕ
取ϕ0=-1，令∆ϕ
i 0
将Newton-Raphsom 法迭代公式中的K T 改用初始矩阵K T =K T (ϕ0) ，就是
修正的Newton-Raphsom 法。仅第一步迭代需要完全求解一个线性方程组，并将
存贮起来，以后的每一步迭代都采用公式 K T
0-1i 0-1
∆ϕi =-(K T ) ψ=(K T ) (H i ϕi +f )
（13）
当∆ϕi 变得充分小，即近似解收敛时，终止迭代。
取ϕ0=1，令∆ϕ
取ϕ0=0，令∆ϕ
取ϕ0=-1，令∆ϕ
可得迭代次数为1次，ϕ1=-0.999664。 （3） 拟Newton-Raphson 方法
K 的修正要满足一下的拟牛顿方程
K i +1(ϕi +1-ϕi ) =ψ(ϕi +1) -ψ(ϕi )
（14）
对于单变量情况，上式中的K i +1是导数(∂ψ∂ϕ)ϕ=ϕi 的近似表达式，实际上就是割线劲度矩阵。
1-1
∆ϕ0=-(K 0) -1ψ0=(K 0) -1(H i ϕi +f )
（15） （16） （17） （18） （19）
ϕ1=ϕ0+∆ϕ0
∆ϕ0ϕ1-ϕ0
(K ) =1=
ψ-ψ0ψ1-ψ0
∆ϕ1=(K 1) -1(H 1ϕ1+f ) (K
i +1-1
∆ϕi ϕi +1-ϕi
) ==i +1 i i
∆ψψ-ψ
当∆ϕi 变得充分小，即近似解收敛时，终止迭代。
取ϕ0=1，令∆ϕ
取ϕ0=0，令∆ϕ
取ϕ0=-1，令∆ϕ
根据结果可知，在精度取∆ϕ
若本题对精度要求不高，主要考虑迭代次数，建议采用Newton 法和拟Newton ；法；若本题对精度和迭代次数要求不高，主要考虑计算量，建议采用修正Newton 法。
Euler-Newton 法
在增量步内采用Newton 迭代法。现以δm 和δm 分别表示第m 级载荷增量时δ
的初值和终值，以R m 表示第m 级增量时的R 的终值，则由式（11）得第m 增量步的迭代公式
i i
∆δm =K T , m -1
δm =δm -1
（20） （21）
R m =λm
()
-1
i i
(λm -F m ) =K T , m -1
()
-1
i （22） (∆λm -ψm )
i +1i i
δm =δm +∆δm
（23）
如果每一增量步内只迭代一次，此时
则对第m 增量步有
∆δm =K T , m -1
1
δm =δm 0∆δm =∆δm
（24） （25）
()
-1
0 (∆λm -ψm )
（26） （27）
δm =δm -1+∆δm
设λ0=0，R 0=0，δ0=0，∆λ1=0. 5，λ1=0. 5，∆λ2=0. 5，λ2=1。设
∆δm
13
=-2. 775099=-0.999664得δ10=0，δ11=-0. 75，δ2，δ22=-1.000000，δ2。于是
可得ϕ=-0. 999664。
附录：
%直接迭代法 clear; y0=1; n=0; for i=1:100;
y1=-10/(10*(1+exp(8*y0))); d=y1-y0; y0=y1;
if (abs(d)>0.0000005); format long,y0 n=n+1; n end end
%newton-raphsom法 clear; y0=1; n=0; for i=1:100;
k=10+10*exp(8*y0)+80*y0*exp(8*y0); f=10*y0+10*y0*exp(8*y0); d=1/k*(-10-f); y0=y0+d; if (abs(d)>0.0000005); format long,y0 n=n+1; n end end
%修正的newton-raphsom 法 clear; y0=1;
a=1; n=0;
for i=1:100000000;
k=10+10*exp(8*a)+80*a*exp(8*a); f=10*y0+10*y0*exp(8*y0); d=1/k*(-10-f); y0=y0+d; if (abs(d)>0.0000005); y0 n=n+1; n end end
%拟newton-raphsom 法 clear; y0=1; n=0; for i=1:100;
f0=10*y0+10*y0*exp(8*y0); k0=10+10*exp(8*y0)+80*y0*exp(8*y0); b=1/k0*(-10-f0); a0=f0+10; y1=y0+b;
f1=10*y1+10*y1*exp(8*y1); a1=f1+10; k1=(a1-a0)/b; d=y1-y0; k0=k1; y0=y1; f0=f1;
if (abs(d)>0.0000005); y0 n=n+1;
n end end
%Euler_Newton法 clear; x0=0; n=0; r=-10; a1=0.5; a2=1; x01=x0;
f01=10*(1+exp(8*x01))*x01+10; k01=10+10*exp(8*x0)+80*x0*exp(8*x0); d11=1/k01*(a1*r-f01); x11=x01+d11; x1=x11; x02=x1;
f02=10*(1+exp(8*x02))*x02+10; k02=10+10*exp(8*x1)+80*x1*exp(8*x1); d12=1/k02*(a2*r-f01); x12=x02+d12; for i=1:100;
k=10+10*exp(8*x12)+80*x12*exp(8*x12); f=10*x12+10*x12*exp(8*x12); d=1/k*(-10-f); x12=x12+d; if (abs(d)>0.0000005); format long,x12 d n=n+1; n end end
（二）、针对软化问题的求解方法 参考文献：《a local arc-length procedure for strain softening》 （1）弧长法
弧长法的约束方程：δ1T ⋅δi =∆l 2(i =1,2,3,...) ；其中∆l 为弧长；δi 为现在荷载增量步第i 次迭代的总的增量位移。
i
δi 的计算式：δi =∑(∆U ) j (i =1,2,3,...)
j =1
以外部荷载系数增量∆λ作为未知量，增量位移向量采用Ramm 和Crisfeld 写成：
(∆U ) i =U F +∆λi ⋅U P
U F =-(K (U i )) -1⋅(R (U i ) -λi ⋅P ) U P =-(K (U i )) -1⋅P
其中：R 为内部力向量；P 为外部力向量；U 为第i 次迭代总的变形向量；λi 为第i 次迭代总的荷载系数。U i 和λi 通过第i 次迭代后用下式计算：
⎧U i =U i -1+(∆U ) i
,(i =1,2,3...) ⎨
⎩λi +1=λi +∆λi
由上述方程，增量荷载系数表示成：
∆λ1=
T U P ⋅(δi -1+U F )
∆λi =∆λ1-(i =2,3,4,...) T
U P ⋅U P
约束方程也改写成：
δ1T ⋅δ1=∆l 2(i =1,2,3,...)
2
δi T ⋅δ=∆l (i =
2,3,4,...) -1i
∆λ1=
T
∆l 2-δi T ⋅(δ-1i -1+U F ) ∆λi =(i =2,3,4,...) T
δi -1⋅U P
（2）局部弧长法
May 和Duan 进一步提出局部弧长法，认为在非线性处用相对位移代替
∆δ=[δ1-δn , δ2-δ1, δ3-δ2, ..., δn -δn -1, ]
则局部弧长中约束方程为
T 2
(∆δ) ⋅(∆δ) =∆l ∑1e i e e =1m
荷载增量表达式为
∆λ1=
m
T P e
∆λi =∆λ1-
∑(∆U
) ⋅(∆δi -1+∆U F ) e
i =2,3,4,...)
作业2：开挖荷载的求解方法
地基开挖时，需要计算开挖荷载，如下图所示：
建立x-y-z 坐标系，z 方向垂直向外，基坑的长为a ，宽为b 。在x-y 平面内为平面应变问题，在y-z 平面内也是平面应变问题。因此有：
⎧E ⎡υ∂u ⎤⎪σx =1+υ⎢1-2υθ+∂x ⎥
⎣⎦⎪
⎪E ⎡υ∂v ⎤σ=θ+⎪y ⎢1-2υ⎥1+υ∂y ⎣⎦⎪⎪E ⎡υ∂w ⎤⎪σ=θ+⎨z ⎢⎥1+υ1-2υ∂z ⎣⎦⎪ （1） ⎪∂u ∂v ∂w ⎪θ=++⎪∂x ∂y ∂z ⎪∂u ∂w
=0⎪=
∂x ∂z ⎪⎩
用位移分量表示的平面应变的形变势能表达式为：
⎡υ⎛∂u ∂v ⎫2⎛∂u ⎫2⎛∂v ⎫21⎛∂v ∂u ⎫2⎤E
V ε=⎢ +⎪+ ⎪+ ⎪+ +⎪⎥dxdy （2）⎰⎰21+υ⎢1-2υ∂x ∂y ∂x ∂y 2∂x ∂y ⎝⎭⎝⎭⎝⎭⎝⎭⎥⎣⎦
根据题意，设位移分量的表达式为：
⎧u =0
⎨2 （3）
v =B y +B y ⎩12
将（3）带入（2）式得到：
V ε=
E 22⎤⎡υ
B +2B y +B +B y dxdy （4） ()()1212⎰⎰⎢⎥21+υ⎣1-2υ⎦
采用利兹变分方法求解位移分量的系数，有：
∂V ε
=⎰⎰f y v m dxdy +⎰f y v m ds （5） ∂B m
（4）式带入（5）式，考虑到：
⎧v 1=y ⎨2 （6） v =y ⎩2
有：
E 1-υ⎧∂V ε2
=(B +2B y ) dxdy =-ρgydxdy =-ρgab /22⎰⎰⎪∂B (1+υ) 1-2υ⎰⎰1
⎪1⎨ （7） ⎪∂V ε=2E 1-υ(B +2B y ) ydxdy =-ρgy 2dxdy =-ρgab 3/3
12⎰⎰⎰⎰⎪∂B (1+υ) 1-2υ⎩2
进而得到：
1+υ)(1-2υ)⎧(⎪B 1=-ρgb
E 1-υ⎪
⎨
⎪B =ρg (1+υ)(1-2υ)⎪22E 1-υ⎩
进而有：
（8）
1+υ)(1-2υ)(∂v
εy ==-ρg (b -y ) （9）
∂y E 1-υ
根据（1）式有：
E υ∂v υ⎧
⎪σx =1+υ1-2υ∂y =-ρg 1-υ(b -y ) ⎪⎨ （10） E 1-υ∂v ⎪σ=-ρg (b -y ) y =⎪1+υ1-2υ∂y ⎩
由上式可得区域内任一点的应力值。进而边界处的面力（开外荷载）即可由上
式得到。
作业3 基于M-C 准则下非线性材料弹塑性矩阵的推导
重力坝平面应变问题分析报告
1. 计算参数的选取与说明 某重力坝段平面图如图1所示。
主要的工程计算参数如下： 几何参数
坝高(m )：H =150 坝顶长(m )：h =10 坝底长(m )：L =100 坝体材料：砼
弹性模量：E c =2⨯105MPa 泊松比：u =
摩擦系数：f =1.2 容重：γ=2.4T /m 3 硬化参数：H '=0
坝基材料：岩石
弹性模量：E c =4⨯105MPa
泊松比：u =0.25 摩擦系数：f =1.0 容重：γ=2.7T /m 3 硬化参数：H '=0
2. 前处理
为了能方便地对图1所示结构进行有限元计算分析，本文针对图一所示坝形，开发了一个专门前处理软件。该软件由Visual Basic 语言编写，软件的功能是通过对AutoCAD 的二次开发，实现坝体参数的设置、网格的自动剖分与显示、边界条件与荷载的手动施加以及数据文件的生成。本文处理问题的思路可用图2所示的流程图描述。
VB 程序的界面如图3 所示。
网格剖分如图4所示，细部节点编号如图5所示。
图3 VB程序界面图
图 4 网格剖分图
图5 细部节点编号示意图
3. 计算
计算条件设置：
（1） 荷载：坝体以及地基受有自重荷载。坝体面收静水压力总用，底面
受扬压力作用，如图6所示。等效静水高度分别为150，150，50。
（2） 约束：坝基底部以及左右两侧均受法向唯一约束（连杆约束）。
（3） 屈服准则：坝体与坝基采用D-P 准则，坝体与坝基接触面所在的坝
基单元（建基面）采用M-C 准则。
（4） 计算参数：荷载增量大小为0.7，最大迭代次数为30。
计算步骤：
（1） 平衡地应力：仅考虑地基部分，施加相应的荷载以及约束，得到各个
高斯点的应力值与节点的位移值。输出结果。
（2） 施加坝体自重：此时程序首先读入上一步的坝基的应力值，坝基的位
移值抛弃不要，施加自重荷载，计算应力与位移。输出结果。
（3） 施加水压力：程序先要读取上一步的应力与位移值，然后施加水压力
荷载进行计算。输出结果。
程序更改说明：计算所采用的程序为教学程序。由于每个计算步骤所读取的信息以及输出的信息都有差别，本文将程序修该成三个版本，每个版本针对于每个计算步。这些版本的主要差别是程序起始出读入文件名称的不同、输出结果文件名称的不同、输出结果的形式的不同（高斯点应力或者单元平均应力）。另外本程序中将屈服准则参数设定为单元信息的的量，以考虑不同单元不同屈服准则的要求。本程序依然只考虑四边形等参元。
4 计算结果与后处理。
以上三个步骤的计算结果以及程序见附录文件。
1. 地应力平衡
从结果文件可以看出，按高度方向，每一层单元的应力是相同的，随着地基深度的加大，y 方向应力值也在不断增大，这说明地应力平衡的结果是具有参考价值的。
2. 坝体自重施加
变形图如图7示，变形效果放大10000倍。345号节点的位移值最大，横向为-0.12E-2m, 纵向为-0.139E-2m 。
图7 施加坝体自重坝体变形图
3. 水压力施加
变形图如图8所示，变形效果放大500倍。345号节点的位移值最大，横向为-0.484 E -1m,纵向为-0.560E-1m 。
图8施加水压力坝体变形图
水压力施加后建基面的应力分布曲线
从以上结果可以看出：
（1） 建基面的三个应力值在沿从坝踵向坝趾方向都在减小。
（2） Y 方向的应力水平要高于X 方向的应力水平以及切应力水平。
（3） 应力值在坝踵处的值相对于其他点都非常大，说明此处发生应力集中现
象。
5 总结。
通过本次计算实例的实践，有主要的收获：
1、进一步提高了有限元前后处理的能力
2、进一步提升了对非线性问题的认，但也也存在不足：对计算结果的分析能力有待加强，前后处理应该进一步学习（等值线图，云图）。
计算固体力学课程作业
专 业 固 体 力 学 学 号 1131301009 姓 名 尹亚川
作业1：
（一）、H ϕ+f =0，其中f =10，H =10(1+e 8ϕ)
（1） 试用直接迭代法，Newton-Raphson 方法，修正Newton-Raphson 方法，拟
Newton-Raphson 方法进行求解并进行比较。 （2） 用Euler-Newton 法计算，f 分2级
求解：
（1）直接迭代法：
于是得近似解
H ϕ+f =0
（1） （2）
H 0=H (ϕ0)
ϕ1=(H 0) -1(-f )
重复这一过程，以第i 次近似解求出第i +1次近似解的迭代公式为
（3）
直到
H i =H (ϕi )
（4） （5）
ϕi +1=(H i ) -1(-f )
∆ϕ=ϕi +1-ϕi
（6）
变得充分小，即近似解收敛时，终止迭代。
取ϕ0=1，令∆ϕ
取ϕ0=0，令∆ϕ
取ϕ0=-1，令∆ϕ
H ϕ+f =0
（7） （8） （9）
ψi =ψ(ϕi ) ≡F (ϕi ) -R =H i ϕi +f ≠0
⎛∂ψi i
K T =K T (ϕi ) ≡ ∂ϕ
⎝
⎫⎪⎪ ⎭
i
i -1i i -1
∆ϕi =-(K T ) ψ=(K T ) (H i ϕi +f )
（10） （11） （12）
⎛∂ψi
K T = ∂ϕ
⎝⎫⎛∂H ϕ⎫⎪⎪= ∂ϕ⎪⎪ ⎭⎝⎭
i i
ϕi +1=ϕi +∆ϕi
当∆ϕi 变得充分小，即近似解收敛时，终止迭代。
取ϕ0=1，令∆ϕ
取ϕ0=0，令∆ϕ
取ϕ0=-1，令∆ϕ
i 0
将Newton-Raphsom 法迭代公式中的K T 改用初始矩阵K T =K T (ϕ0) ，就是
修正的Newton-Raphsom 法。仅第一步迭代需要完全求解一个线性方程组，并将
存贮起来，以后的每一步迭代都采用公式 K T
0-1i 0-1
∆ϕi =-(K T ) ψ=(K T ) (H i ϕi +f )
（13）
当∆ϕi 变得充分小，即近似解收敛时，终止迭代。
取ϕ0=1，令∆ϕ
取ϕ0=0，令∆ϕ
取ϕ0=-1，令∆ϕ
可得迭代次数为1次，ϕ1=-0.999664。 （3） 拟Newton-Raphson 方法
K 的修正要满足一下的拟牛顿方程
K i +1(ϕi +1-ϕi ) =ψ(ϕi +1) -ψ(ϕi )
（14）
对于单变量情况，上式中的K i +1是导数(∂ψ∂ϕ)ϕ=ϕi 的近似表达式，实际上就是割线劲度矩阵。
1-1
∆ϕ0=-(K 0) -1ψ0=(K 0) -1(H i ϕi +f )
（15） （16） （17） （18） （19）
ϕ1=ϕ0+∆ϕ0
∆ϕ0ϕ1-ϕ0
(K ) =1=
ψ-ψ0ψ1-ψ0
∆ϕ1=(K 1) -1(H 1ϕ1+f ) (K
i +1-1
∆ϕi ϕi +1-ϕi
) ==i +1 i i
∆ψψ-ψ
当∆ϕi 变得充分小，即近似解收敛时，终止迭代。
取ϕ0=1，令∆ϕ
取ϕ0=0，令∆ϕ
取ϕ0=-1，令∆ϕ
根据结果可知，在精度取∆ϕ
若本题对精度要求不高，主要考虑迭代次数，建议采用Newton 法和拟Newton ；法；若本题对精度和迭代次数要求不高，主要考虑计算量，建议采用修正Newton 法。
Euler-Newton 法
在增量步内采用Newton 迭代法。现以δm 和δm 分别表示第m 级载荷增量时δ
的初值和终值，以R m 表示第m 级增量时的R 的终值，则由式（11）得第m 增量步的迭代公式
i i
∆δm =K T , m -1
δm =δm -1
（20） （21）
R m =λm
()
-1
i i
(λm -F m ) =K T , m -1
()
-1
i （22） (∆λm -ψm )
i +1i i
δm =δm +∆δm
（23）
如果每一增量步内只迭代一次，此时
则对第m 增量步有
∆δm =K T , m -1
1
δm =δm 0∆δm =∆δm
（24） （25）
()
-1
0 (∆λm -ψm )
（26） （27）
δm =δm -1+∆δm
设λ0=0，R 0=0，δ0=0，∆λ1=0. 5，λ1=0. 5，∆λ2=0. 5，λ2=1。设
∆δm
13
=-2. 775099=-0.999664得δ10=0，δ11=-0. 75，δ2，δ22=-1.000000，δ2。于是
可得ϕ=-0. 999664。
附录：
%直接迭代法 clear; y0=1; n=0; for i=1:100;
y1=-10/(10*(1+exp(8*y0))); d=y1-y0; y0=y1;
if (abs(d)>0.0000005); format long,y0 n=n+1; n end end
%newton-raphsom法 clear; y0=1; n=0; for i=1:100;
k=10+10*exp(8*y0)+80*y0*exp(8*y0); f=10*y0+10*y0*exp(8*y0); d=1/k*(-10-f); y0=y0+d; if (abs(d)>0.0000005); format long,y0 n=n+1; n end end
%修正的newton-raphsom 法 clear; y0=1;
a=1; n=0;
for i=1:100000000;
k=10+10*exp(8*a)+80*a*exp(8*a); f=10*y0+10*y0*exp(8*y0); d=1/k*(-10-f); y0=y0+d; if (abs(d)>0.0000005); y0 n=n+1; n end end
%拟newton-raphsom 法 clear; y0=1; n=0; for i=1:100;
f0=10*y0+10*y0*exp(8*y0); k0=10+10*exp(8*y0)+80*y0*exp(8*y0); b=1/k0*(-10-f0); a0=f0+10; y1=y0+b;
f1=10*y1+10*y1*exp(8*y1); a1=f1+10; k1=(a1-a0)/b; d=y1-y0; k0=k1; y0=y1; f0=f1;
if (abs(d)>0.0000005); y0 n=n+1;
n end end
%Euler_Newton法 clear; x0=0; n=0; r=-10; a1=0.5; a2=1; x01=x0;
f01=10*(1+exp(8*x01))*x01+10; k01=10+10*exp(8*x0)+80*x0*exp(8*x0); d11=1/k01*(a1*r-f01); x11=x01+d11; x1=x11; x02=x1;
f02=10*(1+exp(8*x02))*x02+10; k02=10+10*exp(8*x1)+80*x1*exp(8*x1); d12=1/k02*(a2*r-f01); x12=x02+d12; for i=1:100;
k=10+10*exp(8*x12)+80*x12*exp(8*x12); f=10*x12+10*x12*exp(8*x12); d=1/k*(-10-f); x12=x12+d; if (abs(d)>0.0000005); format long,x12 d n=n+1; n end end
（二）、针对软化问题的求解方法 参考文献：《a local arc-length procedure for strain softening》 （1）弧长法
弧长法的约束方程：δ1T ⋅δi =∆l 2(i =1,2,3,...) ；其中∆l 为弧长；δi 为现在荷载增量步第i 次迭代的总的增量位移。
i
δi 的计算式：δi =∑(∆U ) j (i =1,2,3,...)
j =1
以外部荷载系数增量∆λ作为未知量，增量位移向量采用Ramm 和Crisfeld 写成：
(∆U ) i =U F +∆λi ⋅U P
U F =-(K (U i )) -1⋅(R (U i ) -λi ⋅P ) U P =-(K (U i )) -1⋅P
其中：R 为内部力向量；P 为外部力向量；U 为第i 次迭代总的变形向量；λi 为第i 次迭代总的荷载系数。U i 和λi 通过第i 次迭代后用下式计算：
⎧U i =U i -1+(∆U ) i
,(i =1,2,3...) ⎨
⎩λi +1=λi +∆λi
由上述方程，增量荷载系数表示成：
∆λ1=
T U P ⋅(δi -1+U F )
∆λi =∆λ1-(i =2,3,4,...) T
U P ⋅U P
约束方程也改写成：
δ1T ⋅δ1=∆l 2(i =1,2,3,...)
2
δi T ⋅δ=∆l (i =
2,3,4,...) -1i
∆λ1=
T
∆l 2-δi T ⋅(δ-1i -1+U F ) ∆λi =(i =2,3,4,...) T
δi -1⋅U P
（2）局部弧长法
May 和Duan 进一步提出局部弧长法，认为在非线性处用相对位移代替
∆δ=[δ1-δn , δ2-δ1, δ3-δ2, ..., δn -δn -1, ]
则局部弧长中约束方程为
T 2
(∆δ) ⋅(∆δ) =∆l ∑1e i e e =1m
荷载增量表达式为
∆λ1=
m
T P e
∆λi =∆λ1-
∑(∆U
) ⋅(∆δi -1+∆U F ) e
i =2,3,4,...)
作业2：开挖荷载的求解方法
地基开挖时，需要计算开挖荷载，如下图所示：
建立x-y-z 坐标系，z 方向垂直向外，基坑的长为a ，宽为b 。在x-y 平面内为平面应变问题，在y-z 平面内也是平面应变问题。因此有：
⎧E ⎡υ∂u ⎤⎪σx =1+υ⎢1-2υθ+∂x ⎥
⎣⎦⎪
⎪E ⎡υ∂v ⎤σ=θ+⎪y ⎢1-2υ⎥1+υ∂y ⎣⎦⎪⎪E ⎡υ∂w ⎤⎪σ=θ+⎨z ⎢⎥1+υ1-2υ∂z ⎣⎦⎪ （1） ⎪∂u ∂v ∂w ⎪θ=++⎪∂x ∂y ∂z ⎪∂u ∂w
=0⎪=
∂x ∂z ⎪⎩
用位移分量表示的平面应变的形变势能表达式为：
⎡υ⎛∂u ∂v ⎫2⎛∂u ⎫2⎛∂v ⎫21⎛∂v ∂u ⎫2⎤E
V ε=⎢ +⎪+ ⎪+ ⎪+ +⎪⎥dxdy （2）⎰⎰21+υ⎢1-2υ∂x ∂y ∂x ∂y 2∂x ∂y ⎝⎭⎝⎭⎝⎭⎝⎭⎥⎣⎦
根据题意，设位移分量的表达式为：
⎧u =0
⎨2 （3）
v =B y +B y ⎩12
将（3）带入（2）式得到：
V ε=
E 22⎤⎡υ
B +2B y +B +B y dxdy （4） ()()1212⎰⎰⎢⎥21+υ⎣1-2υ⎦
采用利兹变分方法求解位移分量的系数，有：
∂V ε
=⎰⎰f y v m dxdy +⎰f y v m ds （5） ∂B m
（4）式带入（5）式，考虑到：
⎧v 1=y ⎨2 （6） v =y ⎩2
有：
E 1-υ⎧∂V ε2
=(B +2B y ) dxdy =-ρgydxdy =-ρgab /22⎰⎰⎪∂B (1+υ) 1-2υ⎰⎰1
⎪1⎨ （7） ⎪∂V ε=2E 1-υ(B +2B y ) ydxdy =-ρgy 2dxdy =-ρgab 3/3
12⎰⎰⎰⎰⎪∂B (1+υ) 1-2υ⎩2
进而得到：
1+υ)(1-2υ)⎧(⎪B 1=-ρgb
E 1-υ⎪
⎨
⎪B =ρg (1+υ)(1-2υ)⎪22E 1-υ⎩
进而有：
（8）
1+υ)(1-2υ)(∂v
εy ==-ρg (b -y ) （9）
∂y E 1-υ
根据（1）式有：
E υ∂v υ⎧
⎪σx =1+υ1-2υ∂y =-ρg 1-υ(b -y ) ⎪⎨ （10） E 1-υ∂v ⎪σ=-ρg (b -y ) y =⎪1+υ1-2υ∂y ⎩
由上式可得区域内任一点的应力值。进而边界处的面力（开外荷载）即可由上
式得到。
作业3 基于M-C 准则下非线性材料弹塑性矩阵的推导
重力坝平面应变问题分析报告
1. 计算参数的选取与说明 某重力坝段平面图如图1所示。
主要的工程计算参数如下： 几何参数
坝高(m )：H =150 坝顶长(m )：h =10 坝底长(m )：L =100 坝体材料：砼
弹性模量：E c =2⨯105MPa 泊松比：u =
摩擦系数：f =1.2 容重：γ=2.4T /m 3 硬化参数：H '=0
坝基材料：岩石
弹性模量：E c =4⨯105MPa
泊松比：u =0.25 摩擦系数：f =1.0 容重：γ=2.7T /m 3 硬化参数：H '=0
2. 前处理
为了能方便地对图1所示结构进行有限元计算分析，本文针对图一所示坝形，开发了一个专门前处理软件。该软件由Visual Basic 语言编写，软件的功能是通过对AutoCAD 的二次开发，实现坝体参数的设置、网格的自动剖分与显示、边界条件与荷载的手动施加以及数据文件的生成。本文处理问题的思路可用图2所示的流程图描述。
VB 程序的界面如图3 所示。
网格剖分如图4所示，细部节点编号如图5所示。
图3 VB程序界面图
图 4 网格剖分图
图5 细部节点编号示意图
3. 计算
计算条件设置：
（1） 荷载：坝体以及地基受有自重荷载。坝体面收静水压力总用，底面
受扬压力作用，如图6所示。等效静水高度分别为150，150，50。
（2） 约束：坝基底部以及左右两侧均受法向唯一约束（连杆约束）。
（3） 屈服准则：坝体与坝基采用D-P 准则，坝体与坝基接触面所在的坝
基单元（建基面）采用M-C 准则。
（4） 计算参数：荷载增量大小为0.7，最大迭代次数为30。
计算步骤：
（1） 平衡地应力：仅考虑地基部分，施加相应的荷载以及约束，得到各个
高斯点的应力值与节点的位移值。输出结果。
（2） 施加坝体自重：此时程序首先读入上一步的坝基的应力值，坝基的位
移值抛弃不要，施加自重荷载，计算应力与位移。输出结果。
（3） 施加水压力：程序先要读取上一步的应力与位移值，然后施加水压力
荷载进行计算。输出结果。
程序更改说明：计算所采用的程序为教学程序。由于每个计算步骤所读取的信息以及输出的信息都有差别，本文将程序修该成三个版本，每个版本针对于每个计算步。这些版本的主要差别是程序起始出读入文件名称的不同、输出结果文件名称的不同、输出结果的形式的不同（高斯点应力或者单元平均应力）。另外本程序中将屈服准则参数设定为单元信息的的量，以考虑不同单元不同屈服准则的要求。本程序依然只考虑四边形等参元。
4 计算结果与后处理。
以上三个步骤的计算结果以及程序见附录文件。
1. 地应力平衡
从结果文件可以看出，按高度方向，每一层单元的应力是相同的，随着地基深度的加大，y 方向应力值也在不断增大，这说明地应力平衡的结果是具有参考价值的。
2. 坝体自重施加
变形图如图7示，变形效果放大10000倍。345号节点的位移值最大，横向为-0.12E-2m, 纵向为-0.139E-2m 。
图7 施加坝体自重坝体变形图
3. 水压力施加
变形图如图8所示，变形效果放大500倍。345号节点的位移值最大，横向为-0.484 E -1m,纵向为-0.560E-1m 。
图8施加水压力坝体变形图
水压力施加后建基面的应力分布曲线
从以上结果可以看出：
（1） 建基面的三个应力值在沿从坝踵向坝趾方向都在减小。
（2） Y 方向的应力水平要高于X 方向的应力水平以及切应力水平。
（3） 应力值在坝踵处的值相对于其他点都非常大，说明此处发生应力集中现
象。
5 总结。
通过本次计算实例的实践，有主要的收获：
1、进一步提高了有限元前后处理的能力
2、进一步提升了对非线性问题的认，但也也存在不足：对计算结果的分析能力有待加强，前后处理应该进一步学习（等值线图，云图）。