计算固体力学课程作业.doc

《计算固体力学课程作业.doc》由会员分享，可在线阅读，更多相关《计算固体力学课程作业.doc（58页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、Four short words sum up what has lifted most successful individuals above the crowd: a little bit more.-author-date计算固体力学课程作业计算固体力学课程作业计算固体力学课程作业专 业 固 体 力 学 学 号 1131301009 姓 名 尹亚川 作业1：（一）、，其中，（1） 试用直接迭代法，Newton-Raphson方法，修正Newton-Raphson方法，拟Newton-Raphson方法进行求解并进行比较。（2） 用Euler-Newton法计算，f分2级求解：（1）直接

2、迭代法：（1）（2）于是得近似解（3）重复这一过程，以第i次近似解求出第i+1次近似解的迭代公式为（4）（5）直到（6）变得充分小，即近似解收敛时，终止迭代。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为5次，。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为4次，。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为2次，。（1） Newton-Raphson方法（7）（8）（9）（10）（11）（12）当变得充分小，即近似解收敛时，终止迭代。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为11次，。取，令，运用matla

3、b进行编程求解（代码见附录）。可得迭代次数为4次，。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为1次，。（2） 修正的Newton-Raphson方法将Newton-Raphsom法迭代公式中的改用初始矩阵，就是修正的Newton-Raphsom法。仅第一步迭代需要完全求解一个线性方程组，并将存贮起来，以后的每一步迭代都采用公式（13）当变得充分小，即近似解收敛时，终止迭代。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为122244次，。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为21次，。取，令，运用matlab进行编程求解（代

4、码见附录）。可得迭代次数为1次，。（3） 拟Newton-Raphson方法K的修正要满足一下的拟牛顿方程（14）对于单变量情况，上式中的是导数的近似表达式，实际上就是割线劲度矩阵。（15）（16）（17）（18）（19）当变得充分小，即近似解收敛时，终止迭代。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为11次，。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为4次，。取，令，运用matlab进行编程求解（代码见附录）。可得迭代次数为1次，。根据结果可知，在精度取时，Newton法和拟Newton法迭代次数基本一致，收敛速度较快，而修正的Newton法迭

5、代次数较多，收敛速度较慢。不过，Newton法和拟Newton法计算量较大，而修正Newton法计算量较小。并且，直接迭代法在解决这种简单问题时迭代次数也较少，收敛速度较快。若本题不考虑迭代次数，而对精度要求较高，建议采用Newton法和拟Newton法；若本题对精度要求不高，主要考虑迭代次数，建议采用Newton法和拟Newton；法；若本题对精度和迭代次数要求不高，主要考虑计算量，建议采用修正Newton法。Euler-Newton法在增量步内采用Newton迭代法。现以和分别表示第级载荷增量时的初值和终值，以表示第级增量时的的终值，则由式（11）得第增量步的迭代公式（20）（21）（22

6、）（23）如果每一增量步内只迭代一次，此时（24）（25）则对第m增量步有（26）（27）设，。设，根据Euler-Newton法基本原理运用matlab编程（代码见附录）得，。于是可得。附录：%直接迭代法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*

7、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);y0 n=n+1;nendend%拟newton-raphsom法clear;y0=1;n=0;for

8、 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)*x

9、01+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)

10、0.0000005);format long,x12dn=n+1;nendend（二）、针对软化问题的求解方法参考文献：a local arc-length procedure for strain softening（1）弧长法弧长法的约束方程：；其中为弧长；为现在荷载增量步第i次迭代的总的增量位移。的计算式：以外部荷载系数增量作为未知量，增量位移向量采用Ramm和Crisfeld写成：其中：R为内部力向量；P为外部力向量；U为第i次迭代总的变形向量；为第i次迭代总的荷载系数。和通过第i次迭代后用下式计算：由上述方程，增量荷载系数表示成：约束方程也改写成：（2）局部弧长法May和Duan进一

11、步提出局部弧长法，认为在非线性处用相对位移代替则局部弧长中约束方程为荷载增量表达式为作业2：开挖荷载的求解方法地基开挖时，需要计算开挖荷载，如下图所示：P建立x-y-z坐标系，z方向垂直向外，基坑的长为a，宽为b。在x-y平面内为平面应变问题，在y-z平面内也是平面应变问题。因此有： （1）用位移分量表示的平面应变的形变势能表达式为：（2）根据题意，设位移分量的表达式为： （3）将（3）带入（2）式得到： （4）采用利兹变分方法求解位移分量的系数，有： （5）（4）式带入（5）式，考虑到： （6）有： （7）进而得到： （8）进而有： （9）根据（1）式有： （10）由上式可得区域内任一点的应

12、力值。进而边界处的面力（开外荷载）即可由上式得到。作业3 基于M-C准则下非线性材料弹塑性矩阵的推导重力坝平面应变问题分析报告1. 计算参数的选取与说明某重力坝段平面图如图1所示。h图1 重力坝段几何平面示意图坝体（砼）坝基（岩石）主要的工程计算参数如下：几何参数坝高：坝顶长：坝底长：坝体材料：砼弹性模量：泊松比：摩擦系数：容重：硬化参数：坝基材料：岩石弹性模量：泊松比：摩擦系数：容重：硬化参数：2. 前处理为了能方便地对图1所示结构进行有限元计算分析，本文针对图一所示坝形，开发了一个专门前处理软件。该软件由Visual Basic 语言编写，软件的功能是通过对AutoCAD的二次开发，实现坝

13、体参数的设置、网格的自动剖分与显示、边界条件与荷载的手动施加以及数据文件的生成。本文处理问题的思路可用图2所示的流程图描述。VB程序的界面如图3 所示。网格剖分如图4所示，细部节点编号如图5所示。启动前处理软件，输入控制参数网格自动剖分施加约束与荷载输出单元信息、节点信息、荷载信息、约束信息等文件整合为计算所需要的数据文件读入数据文件计算，输出结果画出位移变形图，结果分析形成报告图2 本文分析流程图图3 VB程序界面图图 4 网格剖分图图5 细部节点编号示意图3. 计算计算条件设置：（1） 荷载：坝体以及地基受有自重荷载。坝体面收静水压力总用，底面受扬压力作用，如图6所示。等效静水高度分别为1

14、50，150，50。（2） 约束：坝基底部以及左右两侧均受法向唯一约束（连杆约束）。（3） 屈服准则：坝体与坝基采用D-P准则，坝体与坝基接触面所在的坝基单元（建基面）采用M-C准则。（4） 计算参数：荷载增量大小为0.7，最大迭代次数为30。计算步骤：（1） 平衡地应力：仅考虑地基部分，施加相应的荷载以及约束，得到各个高斯点的应力值与节点的位移值。输出结果。（2） 施加坝体自重：此时程序首先读入上一步的坝基的应力值，坝基的位移值抛弃不要，施加自重荷载，计算应力与位移。输出结果。（3） 施加水压力：程序先要读取上一步的应力与位移值，然后施加水压力荷载进行计算。输出结果。程序更改说明：计算所采用

15、的程序为教学程序。由于每个计算步骤所读取的信息以及输出的信息都有差别，本文将程序修该成三个版本，每个版本针对于每个计算步。这些版本的主要差别是程序起始出读入文件名称的不同、输出结果文件名称的不同、输出结果的形式的不同（高斯点应力或者单元平均应力）。另外本程序中将屈服准则参数设定为单元信息的的量，以考虑不同单元不同屈服准则的要求。本程序依然只考虑四边形等参元。图 6 静水压力示意图H=150h1=150h2=504 计算结果与后处理。以上三个步骤的计算结果以及程序见附录文件。1. 地应力平衡从结果文件可以看出，按高度方向，每一层单元的应力是相同的，随着地基深度的加大，y方向应力值也在不断增大，这

16、说明地应力平衡的结果是具有参考价值的。2. 坝体自重施加变形图如图7示，变形效果放大10000倍。345号节点的位移值最大，横向为-0.12E-2m,纵向为-0.139E-2m。345图7 施加坝体自重坝体变形图3. 水压力施加变形图如图8所示，变形效果放大500倍。345号节点的位移值最大，横向为-0.484 E -1m,纵向为-0.560E-1m。图8施加水压力坝体变形图水压力施加后建基面的应力分布曲线 从以上结果可以看出：（1） 建基面的三个应力值在沿从坝踵向坝趾方向都在减小。（2） Y方向的应力水平要高于X方向的应力水平以及切应力水平。（3） 应力值在坝踵处的值相对于其他点都非常大，说明此处发生应力集中现象。5 总结。通过本次计算实例的实践，有主要的收获：1、进一步提高了有限元前后处理的能力2、进一步提升了对非线性问题的认，但也也存在不足：对计算结果的分析能力有待加强，前后处理应该进一步学习（等值线图，云图）。-