计算固体力学课程作业
《固体物理》第一章作业题
解 以 H2 为基团，组成 fcc 结构的晶体，如略去动能，分子间按 Lennard—Jones 势相互作
用，则晶体的总相互作用能为：
U = 2N i
Pij −12
R
12
−
j
Pij
−6
R
6
.
Pij−6 = 14.45392; Pij−12 = 12.13188,
→ →→→
c = a1+ a2 − a3
晶列
→
a+
→
b−
2
→
c
可化为
→
a+
→
b−
2
→
c
=
−2
→
a1
+
→
a2
−
2
→
a3
由上式可知，AC晶列在原胞坐标系中的指数为 112
题4.对于晶格常数为a的简单立方晶格，考虑晶格中的一
个晶面(hkl)，证明该晶面所属的晶面族的面间距：
a2 dhkl = h2 + k 2 + l 2
b−
→
c)
2
2
→
BC
=
→
OC −
→
OB
=
→c +
1 2
→
(a+
→b )
−
1 2
→
(b+
→
c)
=
1 2
→
(a+
→
c)
→ → 1 → → → 1→ → a→ → →
BA BC = (2 a+ b− c) (a+ c) = (a− 3 b− c)
固体物理作业2
（a）在t=0时刻，选一电子，令P(t)为该电子在之后的t s内不受
碰撞的几率，电子是时间间隔t~t+dt 内收到碰撞的几率为：
P(t)-P(t+dt)=P(t)dt/τ
()
=−
()
积分得：P(t)=exp(-t/τ)
设P1(t)为该电子在t=0之前的t s内不受碰撞的几率，同理可得：
（c）由（b）知，电子连续两次碰撞时间间隔处于t和t+dt之间
的几率为(dt/τ)*exp(-t/τ)
所以电子连续两次碰撞的平均时间为：
∞
= න
令 = /, =
∞
Байду номын сангаас 0
0
−
−
=
dt/τ
11.20
2）说明并解释金属、半导体电导率随温度变化的差异。
对k空间态密度有：
又因为
代入上式
=
ℏ
=
2
3
(2)
1
=
2
3
(2)
，所以 =
3
ℏ3
1 23
=
()
3
3
2 ℏ
考虑到热平衡状态下体系的均匀性， 0 (, )与r无关。
所以：
23
0 (, ) = 3 ()
亚铁磁：
= +
1
0
−
+ ′
−Θ
； ′为渐近居里点
; Θ为亚铁居里温度
12.4
Fe具有体心立方晶格，原子量=55.85，密度=7.86 g/cm3，
材料固体力学答案1-6章
第一章习题1 证明δ-e 恒等式jtks kt js ist ijke e δδδδ-=[证明]()()()jtks kt js ktjs jtks jtks ktjs jtks kt js itjs jtis ki it ks ktis ji jtks kt js ii ktks ki jtjsjiitis ii ist ijk e e δδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδδ-=-++--=-+---==33习题2 证明若jiij ji ijb b a a -==;，则0=ij ijb a[证明]jiij jiijbb aa-==; jiji ij ij b a b a -=∴，0=+=+∴pq pq ij ij ji ji ij ijb a b a b a b a又因为所有的指标都是哑指标，ijij pq pqb a b a =，所以02=aijbij，即0=ij ijb a习题3 已知某一点的应力分量xxσ，yyσ，zzσ，xyσ不为零，而0==yzxzσσ，试求过该点和z 轴，与x 轴夹角为α的面上的正应力和剪应力。
[解] 如图1.1，过该点和z 轴，与x 轴夹角为α的面的法线，其与x 轴，y 轴和z 轴的方向余弦分别为cos α，sin α，0，则由斜面应力公式的分量表达式，iji jσνσν=)(，可求得该面上的应力为ασασσνσνsin cos 11)(xyxxj j +== ασασσνσνsin cos 22)(yyyxjj +== 033==j j v σνσ)(由斜面正应力表达式ji ij nννσσ=，可求得正应力为ασαασασσ22sinsin cos 2cosyyxyxxn++=剪应力为ασασσστ2cos 2sin )(2122)()(xyxx yynn n +-=-=-=σσσn习题4 如已知物体的表面由0),,(=z y x f 确定，沿物体表面作用着与其外法线方向一致分布载荷()z y x p ,,。
固体力学复习题
固 体 力 学 复 习 题1、 概念1）弹性与塑性弹性：物体在引起形变的外力被除去以后能恢复原形的这一性质。
塑性：物体在引起形变的外力被除去以后有部分变形不能恢复残留下来的这一性质。
2）应力和应力状态应力：受力物体某一截面上一点处的内力集度。
应力状态：某点处的9个应力分量组成的新的二阶张量∑。
3）球张量和偏量球张量：球形应力张量，即σ=000000m m m σσσ⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦，其中()13m x y z σσσσ=++偏量：偏斜应力张量，即x m xy xz ij yxy m yz zx zy z m S σστττσστττσσ⎡⎤-⎢⎥=-⎢⎥⎢⎥-⎣⎦，其中()13m x y z σσσσ=++4）转动张量：表示刚体位移部分，即110221102211022u v u w y x z x v u v w ij x y z y w u w v x z y z W ⎡⎤⎛⎫⎛⎫∂∂∂∂--⎢⎥ ⎪ ⎪ ⎪∂∂∂∂⎢⎥⎝⎭⎝⎭⎢⎥⎛⎫⎛⎫∂∂∂∂⎢⎥=-- ⎪⎪⎢⎥ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎢⎥⎢⎥⎛⎫⎛⎫∂∂∂∂⎢⎥-- ⎪ ⎪ ⎪⎢⎥∂∂∂∂⎝⎭⎝⎭⎣⎦5）应变张量：表示纯变形部分，即112211221122u u v u w x y x z x v u v v w ij x y yz y w u w v wx z y z zε⎡⎤⎛⎫⎛⎫∂∂∂∂∂++⎢⎥ ⎪ ⎪ ⎪∂∂∂∂∂⎢⎥⎝⎭⎝⎭⎢⎥⎛⎫⎛⎫∂∂∂∂∂⎢⎥=++ ⎪⎪⎢⎥ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭⎢⎥⎢⎥⎛⎫⎛⎫∂∂∂∂∂⎢⎥++ ⎪ ⎪ ⎪⎢⎥∂∂∂∂∂⎝⎭⎝⎭⎣⎦6）应变协调条件：物体变形后必须仍保持其整体性和连续性，因此各应变分量之间，必须要有一定得关系，即应变协调条件。
22222y xyx y x x yεγε∂∂∂+=∂∂∂∂。
7）圣维南原理：如作用在弹性体表面上某一不大的局部面积上的力系，为作用在同一局部面积上的另一静力等效力所代替，则荷载的这种重新分布，只造离荷载作用处很近的地方，才使应力的分布发生显著变化，在离荷载较远处只有极小的影响。
计算固体力学第三章_1
8. 可处理大变形和非线形材料带来的非线形问题.
TSINGHUA UNIVERSITY
TSINGHUA UNIVERSITY
3 协调模型分析
1. 建立协调模型的一般方法
大部分有限单元,都是根据虚功原理, 或由它导出的能量 原理建立的, 这类单元统称为“协调模型”或“相容模 型”（Conforming model)。
每个节点有三个转动 分量和三个位移分量.
TSINGHUA UNIVERSITY
TSINGHUA UNIVERSITY
如图1.4, 用120个节点和297个平面应变三角形单 元模拟. 将对称性应用于整个杆端的一半. 此分析 的目的是找出杆端应力集中最高的位置.
TSINGHUA UNIVERSITY
TSINGHUA UNIVERSITY
有限元法无论对什么样的结构（杆系，平面，三维， 板壳）分析过程是一样的，一般为：
有限元法基本步骤:
TSINGHUA UNIVERSITY
有限元法基本步骤
将物体划分为具体有相关节点的等价系统,选择最适当 的单元类型来最接近的模拟实际的物理性能. 所用的单元总 数和给顶物体内单元大小和类型的变化是需要工程判断的 主要问题. 单元必须小到可以给出有用的结果,又必须足够大以节省 计算费用.
一点的位移列阵： 一点的应变列阵：
一点的应力列阵：
一点的体积力列阵： 一点的表面力列阵：
边界外法线方向余弦矩阵：
其中：
平衡方程：（内力与体积力的关系方程）
写成矩阵形式：
其中
A - 微分算子矩阵
几何方程：(应变与位移的关系方程)
写成矩阵形式：
物理方程（应力与应变的关系方程)
北航计算固体力学大作业,网格收敛性分析
Homework of Computational Solid MechanicsElement Type : PLANE 183 with 8 nodesThick of the plane : 0.01mMaterial Properties :EX=2×105Mpa PRXY=0.3 Density=7800kg/m3Case1:Concentrated loads F respectively applied on A and B ,F=10kN Case2:Uniformly distributed load q applied on line AB ,q=20kN/m1、1×1 meshingA andB have the same displacements .Displacements of A、B is shown in table 1 .Table 1Figure 1: X-Component of displacement in Case 1Figure 2: X-Component of displacement in Case 2 2、5×5 meshingDisplacements of A、B is shown in table 2 .Table 2Figure 3: X-Component of displacement in Case 1Figure 4: X-Component of displacement in Case 2 3、10×10 meshingDisplacements of A、B is shown in table 3.Table 34、The first 10 frequenciesWe also get the first 10 frequencies shown in table 4 .Table 4Figure 5: X-Component of displacement in Case 1Figure 6: X-Component of displacement in Case 2Reorganize the data in table 1 ,table 2 and table 3 we get a new table shown as table 5 .Data in table 5 is the displacements in X direction of A,B in different cases and meshings .Table 5Analysis:1、Analyze the data in table 5 ,we can know ,the displacement of A ,B is larger in case 1 than that in case 2 .That means it has a larger displacement when a concentrate load is applied .2、Analyze the data in table 5 , we can also know that the numbers with more meshes are a little larger .In fact , the data will be more accurate if we use more meshes .3、Look at figure 3, figure 4,figure 5and figure 6,we can know that the nodes on line AB have the same displacements when a uniformly distributed load is applied .It’s not suitable for case 1 when a concentrate load is applied .4、Analyze the data in table 4,the frequencies get smaller when the number of meshes is bigger .。
固体物理作业及答案
固体物理作业2.1 光子的波长为20 nm ，求其相应的动量与能量。
答：由λhP =，υh E =得：动量12693410313.3102010626.6----⋅⋅⨯=⨯⋅⨯==m s J ms J hP λ 能量J ms m s J chh E 189183410932.9102010998.210626.6----⨯=⨯⋅⨯⨯⋅⨯===λυ2.2 作一维运动的某粒子的波函数可表达为：, 求归一化常数A? 粒子在何处的几率最大？答：再由2)()(x x ψω=得：222)()(x a x A x -=ω 其中 a x ≤≤0；322222462)(x A x aA x A a dx x d +-=ω 令0)(=dx x d ω得：2,21a x a x ==而a x =1时，0)(=x ω，显然不是最大； 故当22ax =时，粒子的几率最大。
3.1 晶体中原子间的排斥作用和吸引作用有何关系？在什么情况下排斥力和吸引力分别起主导作用？ 答：在原子由分散无规的中性原子结合成规则排列的晶体过程中, 吸引力起到了主要作用. 在吸引力的作用下, 原子间的距离缩小到一定程度, 原子间才出现排斥力. 当排斥力与吸引力相等时, 晶体达到稳定结合状态. 可见, 晶体要达到稳定结合状态, 吸引力与排斥力缺一不可. 设此时相邻原子间的距离为0r , 当相邻原子间的距离0r r 时, 吸引力起主导作用;当相邻原子间的距离0r r 时, 排斥力起主导作用。
3.2已知某晶体中相邻两原子间的相互作用势能可表达为：(1) 求出平衡时两原子间的距离；(2) 平衡时的结合能；(3) 若取m=2, n=10，两原子间的平衡距离为3 Å，晶体的结合能为4 eV/atom 。
请计算出A 和B 的值。
答：设平衡时原子间的距离为0r 。
达到平衡时，相互作用势能应具有最小值，即)(r u 满足：0)(0=∂∂r rr u ，求得mn AmBn r -=10)(……（1） 将0r 代入，得平衡时的结合能mn mn m AmBn AmBn A r u --+-=n 0)(B )()( (2)当m=2，n=10时，由（1）式得5B=A 0r 8，再由0r =3Å，)(0r u -=4eV 代人（2）式可得： 109610001090.54)(m eV r r u B ⋅⨯=-=- 2192000100201050.4)(45)(m eV r r u r u r r A ⋅⨯=-=⎥⎥⎦⎤⎢⎢⎣⎡-=-B4.1 一定温度下，一个光学波的声子数目多，还是声学波的声子数目多？ 答：频率为的格波的(平均) 声子数为：.因为光学波的频率比声学波的频率高, ()大于(), 所以在温度一定情况下, 一个光学波的声子数目少于一个声学波的声子数目.4.2 爱因斯坦模型和德拜模型的主要近似分别是什么？简述德拜温度及其物理意义。
清华大学计算固体力学全套课件
TSINGHUA UNIVERSITY
全套课件
计算固体力学
TSINGHUA UNIVERSITY
第1章 绪论
计算固体力学课程体系
TSINGHUA UNIVERSITY
全面介绍非线性有限元的前沿性内容，使学习 者能进入这一领域的前沿，应用非线性有限元方法 求解弹塑性材料、几何大变形和接触碰撞这些非线 性力学的主要问题，增强工程结构中非线性计算和 虚拟仿真的能力，提高非线性有限元的教学和科研 水平。
TSINGHUA UNIVERSITY
计算固体力学课程体系
教学内容：
1. 绪论：非线性有限元的基本概念，发展历史，工程应用， 标记方法，网格表述和偏微分方程的分类。（2） 2. 一维L有限元：TL和UL格式的控制方程。E有限元：E公式 的控制方程，弱形式与强形式。（4） 3. 连续介质力学：变形和运动，应力－应变的度量，守恒 方程，框架不变性。（4） 4. L网格：UL有限元离散，编制程序，旋转公式。（4） 5. 材料本构模型：一维弹性，非线性弹性，如次弹性和超 弹性。一维塑性，多轴塑性，超弹－塑性（橡胶和泡沫 模型），粘弹性（蠕变和松弛等），经验本构模型，如 J-C方程等。应变硬化和软化。（4） 6. 求解方法：应力更新算法，平衡解答和隐式时间积分 （N-R求解等），显示时间积分（中心差分等） ，波的 传播问题。（4） TSINGHUA UNIVERSITY
Engineering Science- is the systematic acquisition of knowledge for the purpose of applying it to the solution of problems effecting the needs and well-being of human kind. SBES- engineering science and science that employs the principles and methods of modeling and computer simulation to acquire and apply knowledge for the benefit of human kind.
计算固体计算力学 - 第四章 几何非线性问题
。对于某一固定时刻t这种变换可以表示为
* 拉各朗日（Lagrange）描述
t
xi t xi ( 0 x1, 0 x2 , 0 x3 )
基于变形前的构型表述变形后的构型。以变形前的各点坐标 为基本未知数，描述各个量。 根据变形的连续性要求，这种变化必须一一对应，即变换是 单值连续的。同时变换应有唯一的逆变换，也是单值连续的 * 欧拉（Eular）描述
t
t 0 ji
T
t 0
S ji
15
计算固体计算力学
各种应力张量之间的关系： （1）由质量守恒：
t
0
0
V
dV
t
t
V
dV
t
0
V
det( 0 t xi , j ) dV
0 det( t xi , j ) 0 t 0 t 0 t t t t （2） 0 Tji t t x j , m mi ， ji 0 0 x j , m 0Tmi t 0 t t t t 0 0 t （3 ） t S x x Smn x x ji 0 j , m 0 i , n 0 0 0 ji t j ,m t i,n mn ， t
其中：
不能求解
uk
--现实位移分量的变分； --应变的变分； --在现实位形内度量的面积载荷 --在现实位形内度量的体积载荷
t
17
计算固体计算力学
第三节 大变形情况下的本构关系
等温、绝热条件下的小变形线弹性情况，可以用三 种等效的方法描述应力和应变之间的关系
ij Dijkl kl
W ij ij
W
1 Dijkl ij kl 2
应 用 固 体 力 学 大作业
该图为标定长度内的伸长量-载荷曲线，颈缩点坐标为（1838,4.26）。 c） Find the failure time and the total engineering strain (et). Q点是标定长度里最危险的一点，故考虑Q点的应变即可！由于失效已经是在 颈缩点之后， 故不能按照塑性体积不变来求有效应变，这里暂且采用工程应变来 判断。由abaqus输出的Q点的位移除以初始长度D0得到Q点的名义应变,当增量步 由74到75时，该名义应变由0.7136达到0.849，为简化，取增量步75进行考虑。此 时时间time=1.1173，而et=0.228。 d）Calculate the local true stress and true strain in the center area. (Here, let’s define a very small central area CDQE as shown in the above mesh figure) 工程应变：
Problem #1 Build the CAE model and run the simulation. a）Output the Von Mises Stress Contours at 0.9, 1.0 and 1.1 second. 经过计算，该材料的屈服应力为236MPa，并且在总应变达到1时，弹性应变 仍然占不到1%，故可以把上述关系式看成真实塑性应变和真实应力的关系，输 入ABAQUS的材料属性里面。 并且， 为了计算更加准确， 应变每隔0.02取一个值， 以减少abaqus内部的插值导致的误差，结果表明这样得到的Mises应力值要略大 于每隔0.1取的值,这可以从下面的三幅图对比中看出,并且后面的应力应变曲线 也将变得更加光滑。 报告中所有应力的单位为MPa，长度单位为mm，时间为s。 下面三幅图分别是time=0.9, 1.0 and 1.1秒时的Mises Stress云图对比， 坐标是应 变间隔为0.1，右边为0.02：
中科大研究生高等计算固体力学作业参考答案
高等计算固体力学作业参考答案 *解答: 设332210)(x a x a x a a x +++=φ, 余量)()(22x Q dxd x R +=φ由边界条件0)0(=φ, 可得00=a ;由10==Lx dxd φ可得010322321=-++L a L a a(1)(a) 配点法: 取x=L/3和2L/3为配点, 要求:0)3/(=L R (2) 0)3/2(=L R(3)解方程组(1)-(3),可得La a L a 21 ,1 ,2/10321=-=+= (b) 子域法: 取2/0L x ≤≤和L x L ≤≤2/为子域, 则0)(2/0=⎰dx x R L (4) 0)(2/=⎰dx x R LL(5)*任何问题请email to: yqhuang@解方程组(1),(4),(5),可得La a L a 31 ,4/3 ,2/10321=-=+= (b) 伽辽金法. 取权函数33221,,x W x W x W ===,则0)10()(101=--=⎰L x Ldx d W dx x R W φ(6) 0)10()(202=--=⎰L x Ldx d W dx x R W φ(7) 0)10()(303=--=⎰Lx Ldx d W dx x R W φ(8)解方程组(6)-(8),可得La a L a 165,32/23 ,321710321=-=+=解答: 微分算子为) ()() () (2222c y x L +∂∂+∂∂=,取任意函数u, v ,dsnvu ds n u v dxdy u vL ds n v u ds n y u v n x u v dxdy cu y u x u v ds n y v u n x v u dxdy cuv y v y u x v x u dxdy cv y v x v u dxdy v uL y x y x ⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰∂∂+∂∂-=∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-⎪⎪⎭⎫ ⎝⎛+∂∂+∂∂=⎪⎪⎭⎫⎝⎛∂∂+∂∂+⎪⎪⎭⎫ ⎝⎛+∂∂∂∂-∂∂∂∂-=⎪⎪⎭⎫⎝⎛+∂∂+∂∂=)()(22222222故算子是自伴随的.原问题等价于: (假设φ已满足1Γ上的边界条件)0212121212121)()(222222222222222=⎥⎥⎦⎤⎢⎢⎣⎡+⎪⎪⎭⎫ ⎝⎛++⎪⎪⎭⎫ ⎝⎛∂∂-⎪⎭⎫ ⎝⎛∂∂-=+⎪⎪⎭⎫ ⎝⎛++⎪⎪⎭⎫ ⎝⎛∂∂-⎪⎭⎫ ⎝⎛∂∂-=+⎪⎪⎭⎫⎝⎛++⎪⎪⎭⎫ ⎝⎛∂∂∂∂-⎪⎭⎫ ⎝⎛∂∂∂∂-=⎪⎭⎫⎝⎛-∂∂-∂∂+⎪⎪⎭⎫ ⎝⎛++∂∂∂∂-∂∂∂∂-=⎪⎭⎫⎝⎛-∂∂-⎪⎪⎭⎫ ⎝⎛++∂∂+∂∂⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰ΓΓΓΓΓds q dxdy Q c y x ds q dxdy Q c y x ds q dxdy Q c y y x x ds q n ds n dxdy Q c y y x x ds q n dxdy Q c y x φφφφφδδφφφφφδδφδφφδφφδφφδφφδφφδφδφφδφδφφδφφφδφφφφδφ等价的自然变分原理为:()⎰⎰⎰Γ+⎪⎪⎭⎫⎝⎛++⎪⎪⎭⎫ ⎝⎛∂∂-⎪⎭⎫ ⎝⎛∂∂-=∏2222212121ds q dxdy Q c y x φφφφφφ 或()⎰⎰⎰Γ-⎪⎪⎫ ⎛--⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂=∏2221222ds q dxdy Q c y x φφφφφφ解答: 此时问题的变分原理简化为()⎰⎰⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂=∏dxdy y x φφφφ42122 将近似函数代入可以得到：截面的扭矩⎰⎰2Tφ=dxdy解答:()()()()()()0=Γ⎪⎭⎫⎝⎛∂∂+Γ⎪⎭⎫ ⎝⎛+-∂∂+Ω⎥⎦⎤⎢⎣⎡-⎪⎪⎭⎫ ⎝⎛∂∂∂∂-⎪⎭⎫ ⎝⎛∂∂∂∂-=Γ--Γ⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+Ω⎥⎦⎤⎢⎣⎡-⎪⎪⎭⎫⎝⎛∂∂∂∂-⎪⎭⎫ ⎝⎛∂∂∂∂-=Γ--Ω⎥⎦⎤⎢⎣⎡-∂∂∂∂+∂∂∂∂=Γ--Ω⎥⎦⎤⎢⎣⎡-⎪⎪⎭⎫⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂=∏⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰Γ-ΓΓΩΓΓΩΓΩΓΩqq qqq d n k d q n k d Q y k y x k x d q d n y k n x k d Q y k y x k x d q d Q y y k x x k d q d Q y y k x x ky x δφφδφαφφδφφφδφαφδφδφφδφφδφδφφδφφδφαφδφδφδφφδφφδφαφδφδφφδφφδφφδ由变分δφ的任意性,可得相应的欧拉方程和边界条件:0=+⎪⎪⎭⎫⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂Q y k y x k x φφ, Ω内 q Γ上的自然边界条件: 0=+-∂∂q nkαφφq Γ-Γ上的强迫边界条件: 0=δφ,或φφ=解答: 方法1：设A,B 两点的坐标为(x 1,y 1,z 1), (x 2,y 2,z 2), 并设x=x(s),y=y(s),z=z(s), 则dsdxx x ds dy 21--=问题的泛函可以表示为：),()()(11222222z x L ds dsdz ds dx x dz dy dx ds L BAB ABA =+-=++==⎰⎰⎰ 问题转化为求泛函L(x,z)在满足端点条件下的最小值问题。
固体物理作业题
Chapter 1Problem 1.1:Compute the packing fraction f for the bcc lattice.Problem 1.2:(a) Show that the packing fraction f for the diamond lattice is π 3 /16 .(b) What is the packing fraction and coorination number of the honeycomb lattice?Problem 1.3:Consider the hexagonal close packed lattice. (a) Show that c = a 8 3 = 1.633a . Frequently a crystalstructure is called hcp even c is not exactly equal to the ideal value. (b) Show that the packing fraction forthe ideal hcp lattice is π 2 / 6 = 0.7405Problem 1.4:The ionic compound A+B- crystallizes in the NaCl structure. Plot the packing fraction as a function of theratio + −ζ= r / r . Assume that ζ< 1.Problem 1.5:Repeat the calculation of problem 1.4 for the CsCl structure.Problem 1.6:Use the information in the textbook to calculate the densities (in kgm-3) of the following solids: (a)Aluminum, (b) Iron, (c) Silicon and (d) Zinc. Atomic weights of some common elements are listed in thetextbook.Problem 1.7:SrTiO3 crystallizes in the perovskite structure. The strontium atoms are at the corners of the cube with sidea, the titanium atoms are at the body center, while the oxygen atoms occupy the cube faces. (a) What is the Bravais lattice type?2(b) Verify that the primitive unit cell contains one Sr, one Ti and three O atoms.(c) Write down a set of primitive lattice vectors and basis vectors for the perovskite structure.Problem 1.8:The primitive lattice vectors of a certain Bravais lattice can be writtenR n n ax n by n z v r r r1 2 2 1 3( 2 ) 12= 1 + + +What is the lattice type?Problem 1.9:In each of the following cases indicate whether the structure is a Bravais lattice. If it is, give three primitivelattice vectors. If it is not describe it as a Bravais lattice with as small as possible basis. In all cases thelength of the side of the unit cube is a.(a) Base centered cubic (simple cubic with additional points in the centers of the horizontal faces of thecubic cell).(b) Side centered cubic (simple cubic with additional points in the centers of the vertical faces of the cubiccell).(c) Edge centered cubic (simple cubic with additional points at the midpoints of the lines joining nearestneighbors).Problem 1.10:指出体心立方晶格(111) 面与(100) 面，(111) 面与(110) 的交线的晶向。
固体物理作业1
������2������2������+1 ������������2
= ������
������2������+2 + 10������2������ − 11������2������+1
方程组有下列形式的格波解 ������2������ = A������������[������������− 2������ ������������]
11������ − ������������2 ������ − ������−������������������ + 10������������������������ ������������ = 0
ቐ −
10������−������������������
+ ������������������������
������2 = 2 ���������Ԧ��� − 2 ���Ԧ���
������3 = ���������Ԧ���
S1(0,0,0)
Mo1(233 ������, 0, ������)
������2
������11
S2(0,0,2z)
x ������2
S3(233 ������, 0,2������ + ������) Mo2(0,0,3z+w) S4(233 ������, 0,4������ + ������)
23
所以
������100: ������110: ������111 = 6: 3: 2 2 面心立方原子堆积最密集的面是(111)面。
4、在晶体中，由于平移对称性的限制，证明
旋转对称轴只能是1, 2, 3, 4和6重轴，对称元
固体物理习题课_01
02/16
固体物理学——习题课一
面心立方格子原胞基矢: 面心立方格子原胞基矢
倒格子基矢
2π b1 = (i + j + k ) a
2π b3 = (i + j k ) a
2π 同理 b2 = (i j + k ) a
可见由
为基矢构成的格子为体心立方格子
03/16
固体物理学——习题课一
11/16
维格纳 —— 塞茨原胞 —— 14面体 面体 —— 八个面正 六边形 —— 六个面正 四边形
固体物理学——习题课一
12/16
补充例题二:与晶列[l1l2l3]垂直的倒格面的面指数是什么? [解答 解答] 解答 正格子与倒格子互为倒格子. 正格子晶面(h1h2h3)与倒格矢
固体物理学——习题课一
08/16
简单立方晶格 塞茨原胞为原点和6个近邻格点连线的垂直 维格纳 —— 塞茨原胞为原点和 个近邻格点连线的垂直 平分面围成的立方体
固体物理学——习题课一
09/16
面心立方格子 塞茨原胞为原点和12个近邻格点连线的垂直 维格纳 —— 塞茨原胞为原点和 个近邻格点连线的垂直 平分面围成的正十二面体
固体物理学——习题课一
14/16
补充例题四:在晶体衍射中,为什么不能用可见光? [解答 解答] 解答 晶体中原子间距的数量级为 的衍射光栅,光波的波长应小于 米,要使原子晶格成为光波 米. 但可见光的波长为
米, 是晶体中原子间距的1000倍. 因此, 在晶体 衍射中,不能用可见光.
固体物理学——习题课一
固体物理学——习题课一
1.6 如果基矢 的面间距为: 的面间距为:
固体力学作业薄板的振动的固有频率与振型
固体⼒学作业薄板的振动的固有频率与振型固体⼒学作业薄板的振动的固有频率与振型1、问题矩形薄板的参数如下33150,100,5,210,0.3,7.9310/a mm b mm h mm E GPa v kg m ρ======?求矩形薄板在（1）四边简⽀（2）四边固⽀条件下的固有频率和振型2、薄板振动微分⽅程薄板是满⾜⼀定假设的理想⼒学模型，⼀般根据实际的尺⼨和受⼒特点来将某个实际问题简化为薄板模型，如厚度要⽐长、宽的尺⼨⼩得的结构就可以采⽤薄板模型。
薄板在上下表⾯之间存在⼀个对称平⾯，此平⾯称为中⾯，且假定：（1）板的材料由各向同性弹性材料组成；（2）振动时薄板的挠度要⽐它的厚度要⼩；（3）⾃由⾯上的应⼒为零；（4）原来与中⾯正交的横截⾯在变形后始终保持正交，即薄板在变形前中⾯的法线在变形后仍为中⾯的法线。
为了建⽴应⼒、应变和位移之间的关系，取空间直⾓坐标Oxyz ，且坐标原点及xOy 坐标⾯皆放在板变形前的中⾯位置上，如图 1所⽰。
设板上任意⼀点a 的位置，将由变形前的坐标x 、y 、z 来确定。
图 1 薄板模型根据假定（2），板的横向变形和⾯内变形u 、v 是相互独⽴的。
为此，其弯曲变形可由中⾯上各点的横向位移(,,)w x y t 所决定。
根据假定（4），剪切应变分量为零。
由薄板经典理论，可以求得板上任意⼀点(,,)a x y z 沿,,x y z 三个⽅向的位移分量,,a a a u v w 的表达式分别为()a a a w u zx wv zy w w ?=-??=-?=+⾼阶⼩量 (1.1)根据应变与位移的⼏何关系可以求出各点的三个主要是应变分量为222222a x a y a a xyu w z x x v w z y yu v w z y x x yεεγ??==-==-?=+=- (1.2)胡克定律，从⽽获得相对应的三个主要应⼒分量为：2222222222222()()11()()111x x y y y x xy xyE Ez w w y xEz wG x yσεµεµµµσεµεµµµτγµ??=+=-+--=+=-+--???==-+?? (1.3)现画薄板微元的受⼒图如图 2所⽰。
计算固体力学算例
SIGY= -0.100000E+01 TXY = 0.344589E-07 SGMI= -0.100000E+01 CETA= 0.599667E+02 SIGY= -0.100000E+01 TXY = 0.000000E+00 SGMI= -0.100000E+01 CETA= 0.000000E+00 SIGY= -0.100000E+01 TXY = -0.549480E-07 SGMI= -0.100000E+01 CETA= 0.107082E+03 SIGY= -0.100000E+01 TXY = -0.344589E-07 SGMI= -0.100000E+01 CETA= 0.120033E+03
Psw.dat (数据输出)
FINITE ELEMENT ANALYSIS IN PLANE STRESS/STRAIN PROBLEM NJ 6 NE 4 NS 5 NPJ IPS 3 0
PLANE STRESS PROBLEM E 0.200000E+04 PR 0.000000E+00 T 0.100000E+01 V 0.000000E+00
3 2
（a）图
psr.dat (数据输入) 6,4,5,3,0, 2000,0,1,0, 1,2,3, 2,4,5, 2,5,3, 3,5,6, 0.0,2.0, 0.0,1.0, 1.0,1.0, 0.0,0.0, 1.0,0.0, 2.0,0.0, 1,1,0, 2,1,0, 4,1,1, 5,0,1, 6,0,1, 1,-0.5,-0.5, 3,-1.0,-1.0, 6,-0.5,-0.5,
计算固体力学课程作业
计算固体力学课程作业专 业 固 体 力 学 学 号 1131301009 姓 名 尹亚川作业1：（一）、0=+f H ϕ，其中10=f ，)1(108ϕe H +=（1） 试用直接迭代法，Newton-Raphson 方法，修正Newton-Raphson 方法，拟Newton-Raphson 方法进行求解并进行比较。
（2） 用Euler-Newton 法计算，f 分2级求解：（1）直接迭代法： 0=+f H ϕ（1）)(00ϕH H =（2）于是得近似解)()(101f H -=-ϕ（3）重复这一过程，以第i 次近似解求出第i +1次近似解的迭代公式为)(i i H H ϕ=（4） )()(11f H i i -=-+ϕ（5）直到i i ϕϕϕ-=∆+1（6）变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为5次，-0.9996645=ϕ。
取00=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为4次，-0.9996644=ϕ。
取10-=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为2次，-0.9996642=ϕ。
（1） Newton-Raphson 方法 0=+f H ϕ（7） 0)()(≠+=-≡=f H R F i i i i i ϕϕϕψψ（8）ii iT i T K K ⎪⎪⎭⎫ ⎝⎛∂∂≡=ϕψϕ)( （9）)()()(11f H K K i i i T i i T i +=-=∆--ϕψϕ（10） ii iT H K ⎪⎪⎭⎫ ⎝⎛∂∂=⎪⎪⎭⎫ ⎝⎛∂∂=ϕϕϕψ （11）i i i ϕϕϕ∆+=+1（12）当i ϕ∆变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
计算流体力学课程作业
计算流体力学大作业——有限差分法解Poisson 方程五点格式解区域内Poisson 方程摘要：本文结合计算流体力学课上所学知识，采用数值解法中的有限差分法求解Poisson 方程（偏微分方程中椭圆型方程的一种），并用其五点格式采用高斯—塞德尔（Gauss-Seidel ）迭代求解。
并比较了数值近似解与真实解，以及不同步长情况下误差的大小，得到了一定的结论。
关键词：Poisson 方程 有限差分法 五点格式一、计算流体流体力学的特点计算流体力学中许多问题求解最终都会变成偏微分方程的求解，而在数学上，除了几种极少数情况外，要求出它们精确解是很难的。
计算机技术的发展使得这一难题的一很好地解决。
二、偏微分方程的种类2.1、 椭圆型偏微分方程椭圆型偏微分方程的一般形式为()(,)div c u au f x t -∇+= 其中：若12(,,,,)(,)n u u x x x t u x t ==，u ∇为u 的梯度，则其定义为 12,,,n u u x x x ⎡⎤∂∂∂∇=⎢⎥∂∂∂⎣⎦ 散度()div v 的定义为12()n div v v x x x ⎛⎫∂∂∂=+++ ⎪∂∂∂⎝⎭这样，()div c u ∇可以更明确地表示为1122()n n u u u div c u c c c x x x x x x ⎡⎤⎛⎫⎛⎫⎛⎫∂∂∂∂∂∂∇=+++⎢⎥ ⎪ ⎪ ⎪∂∂∂∂∂∂⎝⎭⎝⎭⎝⎭⎣⎦若c 为常数，则进一步化简为 22222212()n div c u c u c u x x x ⎛⎫∂∂∂∇=+++=∆ ⎪∂∂∂⎝⎭其中，∆又称为Laplace 算子。
这样椭圆型偏微分方程可以简单地写为22222212(,)n c u au f x t x x x ⎛⎫∂∂∂-++++= ⎪∂∂∂⎝⎭2.2、抛物型偏微分方程抛物型偏微分方程的一般形式为 ()(,)u d div c u au f x t t∂-∇+=∂ 根据上面叙述，若c 为常数，则该方程可以更简单地写为22222212(,)n u d c u au f x t t x x x ⎛⎫∂∂∂∂-++++= ⎪∂∂∂∂⎝⎭ 2.3、双曲型偏微分方程双曲型偏微分方程的一般形式为22()(,)u d div c u au f x t t∂-∇+=∂ 若c 为常数，则可以将该方程简化为2222222212(,)n u d c u au f x t t x x x ⎛⎫∂∂∂∂-++++= ⎪∂∂∂∂⎝⎭三类方程的直接的区别在于u 对t 的导数的阶次。
ME7001《应用固体力学》课程教学大纲-上海交通大学机械与
ME7001《应用固体力学》课程教学大纲课程名称：应用固体力学课程代码：ME7001学分/学时：3学分/48学时开课学期：春季学期适用专业：机械工程及自动化先修课程：理论力学，弹性力学，有限元开课单位：机械与动力工程学院一、课程性质和教学目标课程介绍：固体力学是开展机械工程相关科学基础研究和工程技术应用需要掌握的重要理论基础，对于提高机械工程专业博士研究生的力学理论基础及其工程应用能力具有重要作用。
本课程面向机械工程博士研究生在科学研究中的固体力学分析需求，讲授连续介质力学基本理论，包括张量分析基础、弹塑性理论、非线性有限元方法，及其在结构和工艺分析中的应用。
教学目标：学生通过学习本课程，可以掌握固体力学的一些基本概念，了解机械工程问题中数学和力学建模、求解的一般原理，初步具备对机械工程中结构和工艺问题进行建模和计算的应用能力，从而为从事机械工程科研工作奠定基础。
具体目标包括：（1）掌握材张量分析理论的基本概念、技术术语。
（2）掌握连续介质力学的基本概念和基本原理。
（3）培养应用固体力学原理解决工程问题和设计满足要求的构件或系统的能力。
二、课程教学内容及学时分配1.固体力学及应用概论（1学时）主要讲述固体力学涉及的理论内容概述，固体力学在机械工程领域科研和工程实践中的应用基本情况。
2.张量分析基础（6学时）主要讲述欧式空间中的矢量和张量、张量和矩阵的几种记法、矢量和张量分析、张量函数的导数、坐标变换、二阶张量及其不变量、Cayley-Hamilton定理、各向同性张量等内容。
3.线弹性问题（6学时）主要讲述各向同性线弹性材料的应力-应变关系、各向异性弹性固体材料的应力-应变关系、弹性刚度张量的对称性、线弹性理论中的变分方法、不变原理和最小势能原理、有限元方法理论、单元插值函数、单元应变、应力、刚度矩阵、边界加载、位移边界条件的引入等。
4.大变形问题基本方程（6学时）主要讲述无限小应变的适用性、物体的变形分析、物体的运动分析、物体的应变度量、物体的应力度量、静力平衡与能量原理、大变形弹性本构方程等。
2017《固体力学原理》作业参考答案(学生版本)
第二章习题 2、证明下列等式（1）123123123123123123ijk i j k j k i k i j j i k k j i i k j εδδδδδδδδδδδδδδδδδδ=++--- 证明：采用穷举法i=1时 右端=3232j k k j δδδδ- j 或k=1 右端=0=左端 j=2，k=2 右端=0=左端 j=2，k=3 右端=1=左端 j=3，k=2 右端= -1=左端 j=3，k=3 右端=0=左端 同理i=2，i=3时也有 右端=0=左端 故 右端=左端 成立 证法2：111222333det[]l m nl m n lmn ij lmna a a a a a a a a a ε=，令ij ij a δ=，得 111222333l m n lmnl m n l m nδδδεδδδδδδ= 推论il im in ijk lmnjl jm jn jl jm jnδδδεεδδδδδδ=（2）ijk ilm jl km jm kl εεδδδδ=-（εδ-关系式）证明：采用穷举法j=1时 k=1 左边=0=右边 k=2 l 或m=3 左边=0=右边 l =1 m=1 左边=0=右边 m=2 左边=1=右边 l =2 m=2 左边=0=右边 m=1 左边= -1=右边 同理k=3 左边=右边同理 j=2，j=3 左边=右边 故 右端=左端 成立 证法2： 运用in il imijk nlm jn jl jm jn jl jmδδδεεδδδδδδ=，令n=i, ，展开行列式证法3：运用恒等式()()()a b c a c b b c a ⨯⨯=⋅-⋅，证明上式（3）il ljk ijk klj ijk δεεεε2==证明：klj kij klj ijk εεεε=il il il lj ij jj il δδδδδδδ23=-=-=3、化简（1）ijk kji εε=-6 （2）mi ijk mkj δεε=-6 （3）ij jk ki δδδ=3 提示：可以直接用上题结论。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算固体力学课程作业专 业 固 体 力 学 学 号 1131301009 姓 名 尹亚川作业1：（一）、0=+f H ϕ，其中10=f ，)1(108ϕe H +=（1） 试用直接迭代法，Newton-Raphson 方法，修正Newton-Raphson 方法，拟Newton-Raphson 方法进行求解并进行比较。
（2） 用Euler-Newton 法计算，f 分2级求解：（1）直接迭代法： 0=+f H ϕ（1）)(00ϕH H =（2）于是得近似解)()(101f H -=-ϕ（3）重复这一过程，以第i 次近似解求出第i +1次近似解的迭代公式为)(i i H H ϕ=（4） )()(11f H i i -=-+ϕ（5）直到i i ϕϕϕ-=∆+1（6）变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为5次，-0.9996645=ϕ。
取00=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为4次，-0.9996644=ϕ。
取10-=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为2次，-0.9996642=ϕ。
（1） Newton-Raphson 方法 0=+f H ϕ（7） 0)()(≠+=-≡=f H R F i i i i i ϕϕϕψψ（8）ii iT i T K K ⎪⎪⎭⎫ ⎝⎛∂∂≡=ϕψϕ)( （9）)()()(11f H K K i i i T i i T i +=-=∆--ϕψϕ（10） ii iT H K ⎪⎪⎭⎫ ⎝⎛∂∂=⎪⎪⎭⎫ ⎝⎛∂∂=ϕϕϕψ （11）i i i ϕϕϕ∆+=+1（12）当i ϕ∆变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为11次，-0.99966411=ϕ。
取00=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为4次，-0.9996644=ϕ。
取10-=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为1次，-0.9996641=ϕ。
（2） 修正的Newton-Raphson 方法将Newton-Raphsom 法迭代公式中的i T K 改用初始矩阵)(00ϕT T K K =，就是修正的Newton-Raphsom 法。
仅第一步迭代需要完全求解一个线性方程组，并将TK 存贮起来，以后的每一步迭代都采用公式)()()(1010f H K K i i T i T i +=-=∆--ϕψϕ（13）当i ϕ∆变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为122244次，-0.986216122244=ϕ。
取00=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为21次，-0.99966321=ϕ。
取10-=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为1次，-0.9996641=ϕ。
（3） 拟Newton-Raphson 方法K 的修正要满足一下的拟牛顿方程)()()(111i i i i i K ϕψϕψϕϕ-=-+++（14）对于单变量情况，上式中的1+i K 是导数()i ϕϕϕψ=∂∂的近似表达式，实际上就是割线劲度矩阵。
)()()(100100f H K K i i +=-=∆--ϕψϕ（15） 001ϕϕϕ∆+=（16） 010101011)(ψψϕϕψψϕ--=-∆=-K（17） )()(11111f H K +=∆-ϕϕ （18）ii ii i i i Kψψϕϕψϕ--=∆∆=++-+1111)( （19）当i ϕ∆变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为11次，-0.99966411=ϕ。
取00=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为4次，-0.9996644=ϕ。
取10-=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为1次，-0.9996641=ϕ。
根据结果可知，在精度取0000005.0<∆ϕ时，Newton 法和拟Newton 法迭代次数基本一致，收敛速度较快，而修正的Newton 法迭代次数较多，收敛速度较慢。
不过，Newton 法和拟Newton 法计算量较大，而修正Newton 法计算量较小。
并且，直接迭代法在解决这种简单问题时迭代次数也较少，收敛速度较快。
若本题不考虑迭代次数，而对精度要求较高，建议采用Newton 法和拟Newton 法；若本题对精度要求不高，主要考虑迭代次数，建议采用Newton 法和拟Newton ；法；若本题对精度和迭代次数要求不高，主要考虑计算量，建议采用修正Newton 法。
Euler-Newton 法在增量步内采用Newton 迭代法。
现以0m δ和m δ分别表示第m 级载荷增量时δ的初值和终值，以m R 表示第m 级增量时的R 的终值，则由式（11）得第m 增量步的迭代公式 10-=m m δδ（20） R R m m λ=（21）()())()(11,11,i m m im T i m m im T i m R K F R K ψλλδ-∆=-=∆---- （22）im i m i m δδδ∆+=+1（23）如果每一增量步内只迭代一次，此时 m m δδ=1（24）m m δδ∆=∆0（25）则对第m 增量步有 ())(011,m m m T m R K ψλδ-∆=∆-- （26）m m m δδδ∆+=-1（27）设00=λ，00=R ，00=δ，5.01=∆λ，5.01=λ，5.02=∆λ，12=λ。
设0000005.0<∆m δ，根据Euler-Newton 法基本原理运用matlab 编程（代码见附录）得001=δ，75.011-=δ，775099.212-=δ，-1.00000022=δ，-0.99966432=δ。
于是可得999664.0-=ϕ。
附录：%直接迭代法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,y0n=n+1;nendend%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,y0n=n+1;nendend%修正的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);y0n=n+1;nendend%拟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);y0n=n+1;nendend%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,x12dn=n+1;nendend（二）、针对软化问题的求解方法 参考文献：《a local arc-length procedure for strain softening 》 （1）弧长法弧长法的约束方程：21(1,2,3,...)T i l i δδ∆⋅==；其中l ∆为弧长；i δ为现在荷载增量步第i 次迭代的总的增量位移。
i δ的计算式：1()(1,2,3,...)ii j j U i δ∆===∑以外部荷载系数增量∆λ作为未知量，增量位移向量采用Ramm 和Crisfeld 写成：()i F i P U U U ∆∆λ=+⋅1(())(())F i i i U K U R U P λ-=-⋅-⋅ 1(())P i U K U P -=-⋅其中：R 为内部力向量；P 为外部力向量；U 为第i 次迭代总的变形向量；i λ为第i 次迭代总的荷载系数。
i U 和i λ通过第i 次迭代后用下式计算：11(),(1,2,3...)i i ii i iU U U i ∆λλ∆λ-+=+⎧=⎨=+⎩ 由上述方程，增量荷载系数表示成：1T PPlU U ∆∆λ=⋅11()(2,3,4,...)T P i F i TP PU U i U U δ∆λ∆λ-⋅+=-=⋅ 约束方程也改写成：211(1,2,3,...)T l i δδ∆⋅==21(2,3,4,...)T i il i δδ∆-⋅== 1TPPlU U ∆∆λ=⋅2111()(2,3,4,...)T Ti i F i Ti Pl U i U ∆δδ∆λδ----⋅+==⋅（2）局部弧长法May 和Duan 进一步提出局部弧长法，认为在非线性处用相对位移代替1,21,32,1,[...,]n n n ∆δδδδδδδδδ-=----则局部弧长中约束方程为211()()mT e i e e l ∆δ∆δ∆=⋅=∑ 荷载增量表达式为11()()mTP e P ee lUU ∆∆λ∆∆==⋅∑1111()()(2,3,4,...)()()mT P ei F ee i m TP e P ee UU i U U ∆∆δ∆∆λ∆λ∆∆-==⋅+=-=⋅∑∑作业2：开挖荷载的求解方法地基开挖时，需要计算开挖荷载，如下图所示：建立x-y-z 坐标系，z 方向垂直向外，基坑的长为a ，宽为b 。