基于線性規(guī)劃模型的極限分析上限有限元程序的malab存儲(chǔ)

1上限法具有目前，利用邊界分析法進(jìn)行巖石工程穩(wěn)定性分析已成為研究的熱點(diǎn)之一。利用該法可進(jìn)行邊坡穩(wěn)定性、地基承載力、隧道土壓力等諸多課題的研究假定破壞模式的方式是利用上限法求解問(wèn)題的傳統(tǒng)途徑。該法直觀明了、概念清晰,對(duì)特定問(wèn)題的研究非常有效。然而,當(dāng)模型的邊界條件、力學(xué)參數(shù)等因素復(fù)雜時(shí),精確的破壞模式難以構(gòu)建,限制了該方法的應(yīng)用。與此同時(shí),采用與有限單元相結(jié)合的上限法近年來(lái)發(fā)展迅速。該方法無(wú)需假定破壞模式,可考慮復(fù)雜的邊界條件、不同的力學(xué)參數(shù),使得其應(yīng)用前景較為廣闊。Sloan等采用非線性規(guī)劃模型可減少所需優(yōu)化變量的數(shù)目,但求解過(guò)程相當(dāng)耗時(shí)。而線性規(guī)劃問(wèn)題求解相對(duì)容易,但實(shí)際應(yīng)用中,為了提高計(jì)算精度,其模型包含的優(yōu)化變量目數(shù)相當(dāng)龐大,將導(dǎo)致計(jì)算和存儲(chǔ)兩方面的困難。為此,本文在Sloan等2滑塊內(nèi)部塑性應(yīng)變的求解極限分析上限定理表述為:在所有的運(yùn)動(dòng)許可的塑性應(yīng)變率場(chǎng)ε式(1)左邊為極限荷載采用傳統(tǒng)的途徑求解時(shí),一般需假定由剛性滑塊體系組成的破壞模式。為計(jì)算簡(jiǎn)便,滑塊內(nèi)部不產(chǎn)生塑性應(yīng)變,耗散功率只產(chǎn)生于速度間斷線上,且速度間斷線上相對(duì)速度方向也預(yù)先確定。而采用上限有限元法,單元內(nèi)部也可發(fā)生塑性應(yīng)變,且速度間斷線上的相對(duì)速度方向未知。3上限和有限面積3.1應(yīng)變率單元應(yīng)用在上限有限元采用的單元中,三角形常應(yīng)變率單元應(yīng)用最廣。如圖1所示,三角形單元每個(gè)節(jié)點(diǎn)有2個(gè)速度分量,每個(gè)三角形單元包含p個(gè)塑性乘子式中:u3.2剛塑性體流動(dòng)法則上限法要求破壞區(qū)域內(nèi)的速度場(chǎng)滿足相關(guān)聯(lián)流動(dòng)法則和速度相容條件,對(duì)于平面應(yīng)變條件下的理想彈塑性或剛塑性體,相關(guān)聯(lián)流動(dòng)法則表示為式中:塑性乘子式中:c和?分別為土體黏結(jié)力和內(nèi)摩擦角,以拉應(yīng)力為正。如文獻(xiàn)[8-9],上式所表示的屈服函數(shù)在以X=σ其中:而α對(duì)于每個(gè)三角形單元,由相關(guān)聯(lián)流動(dòng)法則和速度相容條件,需要施加3個(gè)等式約束,單元內(nèi)部的約束條件可表示為其中:式中:A為單元面積;y3.3速度間斷線上的速度相對(duì)速度三角形常應(yīng)變率單元屬于低階單元,僅利用低階單元離散破壞區(qū)域,常難以達(dá)到理想的精度。為彌補(bǔ)此缺陷,需要在單元之間引入速度間斷線,允許速度發(fā)生跳躍。此時(shí),應(yīng)保證速度間斷線上的相對(duì)速度滿足相關(guān)聯(lián)流動(dòng)法則和相容約束條件。如圖2所示,節(jié)點(diǎn)1、2和3、4所在的速度間斷線與x軸夾角為θ,對(duì)于摩爾-庫(kù)侖屈服準(zhǔn)則,速度間斷線上的法向和切線相對(duì)速度?v和?u需滿足:由于絕對(duì)值的存在難以形成標(biāo)準(zhǔn)的線性規(guī)劃模型。為此,Sloan在速度間斷線相鄰節(jié)點(diǎn)引入新的優(yōu)化變量u盡管采用將速度間斷線節(jié)點(diǎn)相對(duì)速度分量?v和?u由絕對(duì)速度分量表示如下:于是,對(duì)于每條速度間斷線,需要施加塑性流動(dòng)約束條件如下:其中:3.4剛性擋土墻體的速度邊界條件速度場(chǎng)應(yīng)滿足已知的速度邊界條件。設(shè)節(jié)點(diǎn)i位于與x軸夾角為θ的速度邊界上,已知該邊界上的切向和法向速度分別為寫(xiě)成矩陣形式如下:其中:上述速度邊界條件可用于剛性結(jié)構(gòu)施加于土體的速度邊界條件,如剛性基礎(chǔ)或剛性擋土墻等結(jié)構(gòu)。對(duì)于柔性結(jié)構(gòu)作用于土體表面的速度邊界條件,可以對(duì)土體表面施加如下約束:式中:3.5單元內(nèi)部的耗散能量每個(gè)單元內(nèi)部耗散功率將式(7)代入式(24)可得如果假定土體的黏聚力c在單元內(nèi)部也按照線性分布,c3.6速度控制在線延遲的耗散功率每條速度間斷線上的耗散功率式中:|?u|=|u由于速度間斷線l上任一點(diǎn)的u其中:3.7元內(nèi)部線性分布體力如自重在每個(gè)單元上產(chǎn)生的外力功率對(duì)于三角形常應(yīng)變率單元,速度在單元內(nèi)部線性分布,有:寫(xiě)成矩陣形式:對(duì)于面力所作的功率,可按照柔性結(jié)構(gòu)作用于土體表面施加相應(yīng)的速度邊界條件或直接進(jìn)行計(jì)算;當(dāng)作用的面力的分布形式較復(fù)雜時(shí),可將面力沿邊界積分求解相應(yīng)的外力功率,這里不再詳述。4稀疏矩陣存儲(chǔ)上節(jié)已闡述了上限有限元的離散、約束條件的形成以及耗散功率與外力功率的求解等問(wèn)題。接下來(lái)根據(jù)功率平衡方程(1)可得到目標(biāo)函數(shù),即耗散功率減去外力功率的差值。其優(yōu)化變量為節(jié)點(diǎn)速度分量u最小值為式中:大寫(xiě)字母表示各量為全局變量,即考慮所有單元、節(jié)點(diǎn)以及邊界條件。其中X約束條件為實(shí)際計(jì)算中為保證精度,需采用大量三角形單元,并增加每個(gè)單元的塑性乘子由于單元內(nèi)部約束條件和相鄰單元并無(wú)關(guān)系,同樣速度間斷線約束條件和相鄰速度間斷線并無(wú)關(guān)系,這就決定了等式約束的系數(shù)矩陣具有高度稀疏性,而且隨單元數(shù)目的增加,其稀疏性相應(yīng)增加。為此,在編程中采用稀疏矩陣存儲(chǔ)的方式存儲(chǔ)等式約束的系數(shù)矩陣,即僅存儲(chǔ)非0元素的數(shù)值及在矩陣中的行號(hào)和列號(hào),不存儲(chǔ)其余的0元素,采用該法能大大降低計(jì)算過(guò)程中的內(nèi)存占用量。采用MATLAB編制上限有限元程序,以稀疏矩陣存儲(chǔ)的方式在程序中形成相應(yīng)系數(shù)矩陣時(shí)予以考慮。通過(guò)調(diào)用優(yōu)化工具箱linprog函數(shù)求解線性規(guī)劃問(wèn)題。該函數(shù)采用大型算法時(shí),能很好地支持稀疏矩陣,并且無(wú)需設(shè)置優(yōu)化變量的初始值。由此可解決較大規(guī)模線性規(guī)劃問(wèn)題,其中優(yōu)化變量個(gè)數(shù)可達(dá)幾十萬(wàn),滿足一般問(wèn)題的精度要求。以下以條形基礎(chǔ)地基承載力課題為例,進(jìn)行說(shuō)明。5模型網(wǎng)格結(jié)構(gòu)條形基礎(chǔ)地基承載力問(wèn)題是巖土工程中的經(jīng)典課題。這里采用上限有限元求解條形基礎(chǔ)地基承載力系數(shù)N采用3種不同的網(wǎng)格對(duì)計(jì)算模型進(jìn)行劃分,如圖4~6所示。圖4和5示意的模型網(wǎng)格1和3均為規(guī)則網(wǎng)格單元,其區(qū)別在于網(wǎng)格密度不同?？紤]到剛性基礎(chǔ)邊緣為奇異點(diǎn),該點(diǎn)附近區(qū)域速度場(chǎng)變化劇烈,于是,圖6示意的模型網(wǎng)格在基礎(chǔ)邊緣進(jìn)行加密,方法為網(wǎng)格大小沿著x和y方向設(shè)置一定的梯度,形成漸變的網(wǎng)格。表1為3種模型網(wǎng)格的參數(shù)。從表1可以看出,3種不同密度的網(wǎng)格的參數(shù)相差很大。當(dāng)p=36時(shí),網(wǎng)格1和網(wǎng)格3的優(yōu)化變量數(shù)目(節(jié)點(diǎn)速度和塑性乘子之和)分別為122880和240000,而網(wǎng)格2更達(dá)到768000(此處的優(yōu)化變量包括了模型邊界上已知的節(jié)點(diǎn)速度,數(shù)目較少,在計(jì)算前可去除)。在這些優(yōu)化變量中,塑性乘子占相當(dāng)大的比例?？梢?jiàn),將屈服函數(shù)線性化引入了大量?jī)?yōu)化變量。計(jì)算過(guò)程等式約束的系數(shù)矩陣相當(dāng)龐大。如上所述,通過(guò)引入稀疏矩陣存儲(chǔ)的方式,能大大降低計(jì)算過(guò)程中的內(nèi)存需求。同時(shí),由于采用線性規(guī)劃模型,其求解速度較快,該計(jì)算模型可在目前主流PC(配備2G內(nèi)存)上進(jìn)行。表2、3和圖7、8為?=30從表2、3可以看出,對(duì)于網(wǎng)格1和2,計(jì)算結(jié)果的精度隨著網(wǎng)格數(shù)目和塑性乘子數(shù)目p的增加而增加。對(duì)于網(wǎng)格3,其網(wǎng)格數(shù)目小于網(wǎng)格2,然而其基礎(chǔ)邊緣網(wǎng)格加密,得到的N從表2、3和圖7、8還可看出,塑性乘子數(shù)目p對(duì)計(jì)算結(jié)果影響較大。對(duì)于N為了詳細(xì)說(shuō)明經(jīng)過(guò)優(yōu)化得到地基的破壞模式,將速度場(chǎng)計(jì)算結(jié)果示于圖9、10。該速度場(chǎng)的繪制通過(guò)MATLAB計(jì)算結(jié)果的數(shù)據(jù)經(jīng)處理后導(dǎo)入CAD進(jìn)行繪圖得到。圖9為p=24時(shí),利用3種網(wǎng)格計(jì)算N圖10為p=24時(shí),利用3種網(wǎng)格計(jì)算N6在一定程度上解決了等式約束的系數(shù)矩陣的高度稀疏性本文在Sloan等工作的基礎(chǔ)上,基于線性規(guī)劃模