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

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

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

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

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

5、求不高，主要考慮迭代次數(shù)，建議采用Newton法和擬Newton；法；若本題對(duì)精度和迭代次數(shù)要求不高，主要考慮計(jì)算量，建議采用修正Newton法。Euler-Newton法在增量步內(nèi)采用Newton迭代法?，F(xiàn)以和分別表示第級(jí)載荷增量時(shí)的初值和終值，以表示第級(jí)增量時(shí)的的終值，則由式（11）得第增量步的迭代公式（20）（21）（22）（23）如果每一增量步內(nèi)只迭代一次，此時(shí)（24）（25）則對(duì)第m增量步有（26）（27）設(shè)，。設(shè)，根據(jù)Euler-Newton法基本原理運(yùn)用matlab編程（代碼見(jiàn)附錄）得，。于是可得。附錄：%直接迭代法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（二）、針對(duì)軟化問(wèn)題的求解方法參考文獻(xiàn)：a local arc-length procedure for strain softening（1）弧長(zhǎng)法

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

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

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

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

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

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