計算固體力學(xué)課程作業(yè)

1、計算固體力學(xué)課程作業(yè)專 業(yè) 固 體 力 學(xué) 學(xué) 號 1131301009 姓 名 尹亞川 作業(yè)1：（一）、，其中，（1） 試用直接迭代法，Newton-Raphson方法，修正Newton-Raphson方法，擬Newton-Raphson方法進行求解并進行比較。（2） 用Euler-Newton法計算，f分2級求解：（1）直接迭代法：（1）（2）于是得近似解（3）重復(fù)這一過程，以第i次近似解求出第i+1次近似解的迭代公式為（4）（5）直到（6）變得充分小，即近似解收斂時，終止迭代。取，令，運用matlab進行編程求解（代碼見附錄）。可得迭代次數(shù)為5次，。取，令，運用matlab進行編程求解（

2、代碼見附錄）?？傻玫螖?shù)為4次，。取，令，運用matlab進行編程求解（代碼見附錄）?？傻玫螖?shù)為2次，。（1） Newton-Raphson方法（7）（8）（9）（10）（11）（12）當變得充分小，即近似解收斂時，終止迭代。取，令，運用matlab進行編程求解（代碼見附錄）?？傻玫螖?shù)為11次，。取，令，運用matlab進行編程求解（代碼見附錄）。可得迭代次數(shù)為4次，。取，令，運用matlab進行編程求解（代碼見附錄）?？傻玫螖?shù)為1次，。（2） 修正的Newton-Raphson方法將Newton-Raphsom法迭代公式中的改用初始矩陣，就是修正的Newton-Raphsom法

3、。僅第一步迭代需要完全求解一個線性方程組，并將存貯起來，以后的每一步迭代都采用公式（13）當變得充分小，即近似解收斂時，終止迭代。取，令，運用matlab進行編程求解（代碼見附錄）。可得迭代次數(shù)為122244次，。取，令，運用matlab進行編程求解（代碼見附錄）?？傻玫螖?shù)為21次，。取，令，運用matlab進行編程求解（代碼見附錄）?？傻玫螖?shù)為1次，。（3） 擬Newton-Raphson方法K的修正要滿足一下的擬牛頓方程（14）對于單變量情況，上式中的是導(dǎo)數(shù)的近似表達式，實際上就是割線勁度矩陣。（15）（16）（17）（18）（19）當變得充分小，即近似解收斂時，終止迭代。取，令，

4、運用matlab進行編程求解（代碼見附錄）?？傻玫螖?shù)為11次，。取，令，運用matlab進行編程求解（代碼見附錄）?？傻玫螖?shù)為4次，。取，令，運用matlab進行編程求解（代碼見附錄）?？傻玫螖?shù)為1次，。根據(jù)結(jié)果可知，在精度取時，Newton法和擬Newton法迭代次數(shù)基本一致，收斂速度較快，而修正的Newton法迭代次數(shù)較多，收斂速度較慢。不過，Newton法和擬Newton法計算量較大，而修正Newton法計算量較小。并且，直接迭代法在解決這種簡單問題時迭代次數(shù)也較少，收斂速度較快。若本題不考慮迭代次數(shù)，而對精度要求較高，建議采用Newton法和擬Newton法；若本題對精度要

5、求不高，主要考慮迭代次數(shù)，建議采用Newton法和擬Newton；法；若本題對精度和迭代次數(shù)要求不高，主要考慮計算量，建議采用修正Newton法。Euler-Newton法在增量步內(nèi)采用Newton迭代法?，F(xiàn)以和分別表示第級載荷增量時的初值和終值，以表示第級增量時的的終值，則由式（11）得第增量步的迭代公式（20）（21）（22）（23）如果每一增量步內(nèi)只迭代一次，此時（24）（25）則對第m增量步有（26）（27）設(shè)，。設(shè)，根據(jù)Euler-Newton法基本原理運用matlab編程（代碼見附錄）得，。于是可得。附錄：%直接迭代法clear;y0=1;n=0;for i=1:100;y1=-1

6、0/(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=

7、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 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*

8、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;k0

9、2=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）弧長法

10、弧長法的約束方程：；其中為弧長；為現(xiàn)在荷載增量步第i次迭代的總的增量位移。的計算式：以外部荷載系數(shù)增量作為未知量，增量位移向量采用Ramm和Crisfeld寫成：其中：R為內(nèi)部力向量；P為外部力向量；U為第i次迭代總的變形向量；為第i次迭代總的荷載系數(shù)。和通過第i次迭代后用下式計算：由上述方程，增量荷載系數(shù)表示成：約束方程也改寫成：（2）局部弧長法May和Duan進一步提出局部弧長法，認為在非線性處用相對位移代替則局部弧長中約束方程為荷載增量表達式為作業(yè)2：開挖荷載的求解方法地基開挖時，需要計算開挖荷載，如下圖所示：P建立x-y-z坐標系，z方向垂直向外，基坑的長為a，寬為b。在x-y平面內(nèi)為

11、平面應(yīng)變問題，在y-z平面內(nèi)也是平面應(yīng)變問題。因此有： （1）用位移分量表示的平面應(yīng)變的形變勢能表達式為：（2）根據(jù)題意，設(shè)位移分量的表達式為： （3）將（3）帶入（2）式得到： （4）采用利茲變分方法求解位移分量的系數(shù)，有： （5）（4）式帶入（5）式，考慮到： （6）有： （7）進而得到： （8）進而有： （9）根據(jù)（1）式有： （10）由上式可得區(qū)域內(nèi)任一點的應(yīng)力值。進而邊界處的面力（開外荷載）即可由上式得到。作業(yè)3 基于M-C準則下非線性材料彈塑性矩陣的推導(dǎo)重力壩平面應(yīng)變問題分析報告1. 計算參數(shù)的選取與說明某重力壩段平面圖如圖1所示。h圖1 重力壩段幾何平面示意圖壩體（砼）壩基（巖石

12、）主要的工程計算參數(shù)如下：幾何參數(shù)壩高：壩頂長：壩底長：壩體材料：砼彈性模量：泊松比：摩擦系數(shù)：容重：硬化參數(shù)：壩基材料：巖石彈性模量：泊松比：摩擦系數(shù)：容重：硬化參數(shù)：2. 前處理為了能方便地對圖1所示結(jié)構(gòu)進行有限元計算分析，本文針對圖一所示壩形，開發(fā)了一個專門前處理軟件。該軟件由Visual Basic 語言編寫，軟件的功能是通過對AutoCAD的二次開發(fā)，實現(xiàn)壩體參數(shù)的設(shè)置、網(wǎng)格的自動剖分與顯示、邊界條件與荷載的手動施加以及數(shù)據(jù)文件的生成。本文處理問題的思路可用圖2所示的流程圖描述。VB程序的界面如圖3 所示。網(wǎng)格剖分如圖4所示，細部節(jié)點編號如圖5所示。啟動前處理軟件，輸入控制參數(shù)網(wǎng)格自

13、動剖分施加約束與荷載輸出單元信息、節(jié)點信息、荷載信息、約束信息等文件整合為計算所需要的數(shù)據(jù)文件讀入數(shù)據(jù)文件計算，輸出結(jié)果畫出位移變形圖，結(jié)果分析形成報告圖2 本文分析流程圖圖3 VB程序界面圖圖 4 網(wǎng)格剖分圖圖5 細部節(jié)點編號示意圖3. 計算計算條件設(shè)置：（1） 荷載：壩體以及地基受有自重荷載。壩體面收靜水壓力總用，底面受揚壓力作用，如圖6所示。等效靜水高度分別為150，150，50。（2） 約束：壩基底部以及左右兩側(cè)均受法向唯一約束（連桿約束）。（3） 屈服準則：壩體與壩基采用D-P準則，壩體與壩基接觸面所在的壩基單元（建基面）采用M-C準則。（4） 計算參數(shù)：荷載增量大小為0.7，最大迭

14、代次數(shù)為30。計算步驟：（1） 平衡地應(yīng)力：僅考慮地基部分，施加相應(yīng)的荷載以及約束，得到各個高斯點的應(yīng)力值與節(jié)點的位移值。輸出結(jié)果。（2） 施加壩體自重：此時程序首先讀入上一步的壩基的應(yīng)力值，壩基的位移值拋棄不要，施加自重荷載，計算應(yīng)力與位移。輸出結(jié)果。（3） 施加水壓力：程序先要讀取上一步的應(yīng)力與位移值，然后施加水壓力荷載進行計算。輸出結(jié)果。程序更改說明：計算所采用的程序為教學(xué)程序。由于每個計算步驟所讀取的信息以及輸出的信息都有差別，本文將程序修該成三個版本，每個版本針對于每個計算步。這些版本的主要差別是程序起始出讀入文件名稱的不同、輸出結(jié)果文件名稱的不同、輸出結(jié)果的形式的不同（高斯點應(yīng)力或

15、者單元平均應(yīng)力）。另外本程序中將屈服準則參數(shù)設(shè)定為單元信息的的量，以考慮不同單元不同屈服準則的要求。本程序依然只考慮四邊形等參元。圖 6 靜水壓力示意圖H=150h1=150h2=504 計算結(jié)果與后處理。以上三個步驟的計算結(jié)果以及程序見附錄文件。1. 地應(yīng)力平衡從結(jié)果文件可以看出，按高度方向，每一層單元的應(yīng)力是相同的，隨著地基深度的加大，y方向應(yīng)力值也在不斷增大，這說明地應(yīng)力平衡的結(jié)果是具有參考價值的。2. 壩體自重施加變形圖如圖7示，變形效果放大10000倍。345號節(jié)點的位移值最大，橫向為-0.12E-2m,縱向為-0.139E-2m。345圖7 施加壩體自重壩體變形圖3. 水壓力施加變形圖如圖8所示，變形效果放大500倍。345號節(jié)點的位移值最大，橫向為-0.484 E -1m,縱向為-0.560E-1