本發(fā)明屬于結(jié)構(gòu)分析領(lǐng)域，具體涉及基于彈塑性分解的非線性有限元剛度矩陣更新方法。

背景技術(shù)：

在結(jié)構(gòu)分析領(lǐng)域，有限元方法作為一種應(yīng)用廣泛的數(shù)值計(jì)算方法，由于其通用性和有效性，目前已成為計(jì)算機(jī)輔助設(shè)計(jì)及數(shù)值仿真的重要實(shí)現(xiàn)手段。對于極端環(huán)境荷載作用下的結(jié)構(gòu)響應(yīng)計(jì)算，需考慮材料的非線性屬性，由于此時求解方程中剛度矩陣需根據(jù)當(dāng)前的結(jié)構(gòu)非線性狀態(tài)實(shí)時變化，因此通常采用步進(jìn)式的增量格式及離散的線性化方法進(jìn)行求解，這一過程中迭代計(jì)算不可避免，且需不斷進(jìn)行整體剛度矩陣的更新和分解運(yùn)算。當(dāng)計(jì)算模型具有較大規(guī)?；蜉^為精細(xì)的單元劃分時，大規(guī)模剛度矩陣的實(shí)時更新和分解過程將極為耗時，成為制約非線性有限元技術(shù)應(yīng)用的主要因素。

技術(shù)實(shí)現(xiàn)要素：

為了克服上述現(xiàn)有方法的不足，本發(fā)明通過對材料應(yīng)變進(jìn)行分解，提出了一種新的用于材料非線性問題的剛度矩陣更新方法。本發(fā)明可將進(jìn)入非線性變形狀態(tài)的局部結(jié)構(gòu)區(qū)域與整體彈性結(jié)構(gòu)相互分離，求解時僅需對代表局部非線性變形的小規(guī)模矩陣進(jìn)行更新，從而將傳統(tǒng)方法中由局部非線性變形引起的大規(guī)模剛度矩陣更新過程轉(zhuǎn)變?yōu)閷Τ跏紡椥詣偠染仃嚨牡椭葦z動過程。本發(fā)明具有較為廣泛的適用性，可用于眾多領(lǐng)域的非線性結(jié)構(gòu)分析中。

本發(fā)明的技術(shù)方案：

一種基于彈塑性分解的非線性有限元剛度矩陣更新方法，步驟如下：

(1)應(yīng)變分解

對于具有非線性本構(gòu)關(guān)系的材料，在任意狀態(tài)下其應(yīng)力和應(yīng)變分別用σ和ε表示，將應(yīng)變ε分解為線彈性應(yīng)變ε′與非線性應(yīng)變ε"之和的形式，即：

ε＝ε′+ε″(1)

式中，線彈性應(yīng)變ε′定義為按照材料的初始彈性屬性加載至應(yīng)力σ時對應(yīng)的應(yīng)變，非線性應(yīng)變ε"定義為總應(yīng)變ε與線彈性應(yīng)變ε′的差，應(yīng)力、線彈性應(yīng)變和非線性應(yīng)變之間符合如下關(guān)系：

σ＝deε′＝de(ε-ε″)(2)

式中，de代表材料的初始彈性本構(gòu)矩陣，對于單軸材料，其退化為材料初始彈性模量；

(2)使用插值方法建立單元非線性變形機(jī)制

對于待分析結(jié)構(gòu)，建立對應(yīng)的有限元數(shù)值模型，并在單元內(nèi)設(shè)置若干非線性變形插值結(jié)點(diǎn)，從而單元域內(nèi)任意一點(diǎn)的非線性應(yīng)變可通過非線性應(yīng)變場的插值表達(dá)式獲得，此時，單元非線性應(yīng)變場表示為如下插值形式：

ε″＝cε″p(3)

式中，c為非線性插值函數(shù)矩陣；向量ε″p集中了單元內(nèi)所有非線性變形插值點(diǎn)的非線性應(yīng)變；

(3)在任意增量求解步中，僅考慮模型結(jié)構(gòu)中進(jìn)入非線性變形狀態(tài)的局部區(qū)域，并使用增量格式，建立分塊形式的結(jié)構(gòu)控制方程：

式中，ke為模型結(jié)構(gòu)的初始彈性剛度矩陣；向量δf、δx分別代表結(jié)構(gòu)的荷載增量和位移增量；向量δε″pr為模型中進(jìn)入非線性變形狀態(tài)的局部子結(jié)構(gòu)非線性應(yīng)變增量；k″pr為進(jìn)入非線性變形狀態(tài)的局部子結(jié)構(gòu)剛度矩陣；k′r為整體彈性結(jié)構(gòu)與局部非線性結(jié)構(gòu)之間的交互矩陣；

(4)建立模型結(jié)構(gòu)的非線性求解方程，即：

步驟1中，應(yīng)變分解過程可適用于任意形式的非線性材料本構(gòu)關(guān)系，當(dāng)卸載剛度與初始彈性剛度相同時，步驟1中定義的線彈性應(yīng)變與非線性應(yīng)變將分別等于彈性應(yīng)變與塑性應(yīng)變。

步驟2中，由插值方法建立的單元非線性變形機(jī)制不局限于特定的單元類型，但非線性變形插值點(diǎn)的個數(shù)和位置影響單元計(jì)算精度。

步驟3和步驟4中，需首先判斷結(jié)構(gòu)中所有單元非線性變形插值點(diǎn)的非線性狀態(tài)，并在構(gòu)造矩陣k′r和k″pr的過程中使其僅包含處于非線性變形狀態(tài)的單元和插值點(diǎn)信息。

步驟3和步驟4中，矩陣ke代表整體彈性結(jié)構(gòu)，在分析過程中始終不變；向量δε″pr由所有位于局部非線性變形區(qū)域內(nèi)的非線性插值點(diǎn)的非線性應(yīng)變集成得到，其中不包括處于線彈性狀態(tài)的插值點(diǎn)信息；k″pr代表進(jìn)入非線性變形狀態(tài)的局部子結(jié)構(gòu)，包含位于其中的非線性插值點(diǎn)材料非線性本構(gòu)信息，需在分析時根據(jù)材料非線性狀態(tài)實(shí)時更新；矩陣k′r與數(shù)值模型中產(chǎn)生非線性變形局部區(qū)域的位置有關(guān)，而與材料的非線性程度無關(guān)，其中的每個元素僅需計(jì)算一次，無需在非線性求解過程中重復(fù)計(jì)算。

步驟4中，求解方程式(5)通過對式(4)進(jìn)行凝聚獲得，等號左邊括號中的矩陣表達(dá)式為初始彈性剛度矩陣的低秩攝動，與整體結(jié)構(gòu)的切向剛度矩陣等效。

本發(fā)明的有益效果：

1.非線性變形狀態(tài)的分離。對于非線性僅產(chǎn)生于結(jié)構(gòu)局部的有限元計(jì)算問題，本發(fā)明在控制方程中隔離出了代表非線性變形狀態(tài)的小規(guī)模矩陣，進(jìn)而在非線性求解時僅需對這一小規(guī)模矩陣進(jìn)行更新，避免了對大規(guī)模整體剛度矩陣的更新。

2.廣泛的適用范圍。本發(fā)明所提方法在結(jié)構(gòu)非線性分析時不局限于特定的材料本構(gòu)及單元類型，可應(yīng)用到眾多學(xué)科及研究領(lǐng)域，例如土木、水利、機(jī)械、力學(xué)等，因而具有廣泛的適用性。

附圖說明

圖1為實(shí)施例一的應(yīng)變分解示意圖。

圖2為實(shí)施例二的有限元模型。

圖3為實(shí)施例二的單元模型。

圖4為實(shí)施例二的單元內(nèi)任一點(diǎn)的應(yīng)力與應(yīng)變。

具體實(shí)施方式

本發(fā)明提供了一種基于應(yīng)變分解思想的剛度矩陣更新方法，以下結(jié)合附圖實(shí)例和技術(shù)方案，進(jìn)一步說明本發(fā)明的具體實(shí)施方式。附圖和實(shí)施例僅為說明本發(fā)明的實(shí)施方式，不構(gòu)成對本發(fā)明的任何限制。

實(shí)施例一：材料應(yīng)變分解實(shí)施方式

以圖1所示單軸非線性材料的應(yīng)力-應(yīng)變關(guān)系為例說明本發(fā)明中材料應(yīng)變分解的實(shí)施方式。圖中c點(diǎn)代表當(dāng)前的應(yīng)力應(yīng)變狀態(tài)，其應(yīng)力和應(yīng)變分別為σ和ε，ee代表材料初始彈性模量，延長直線oa交應(yīng)力σ于b點(diǎn)，定義b點(diǎn)橫坐標(biāo)為線彈性應(yīng)變ε′，定義ε與線彈性應(yīng)變ε′的差為非線性應(yīng)變ε"，即：

ε″＝ε-ε′

實(shí)施例二：非線性問題求解

以圖2所示某個由三結(jié)點(diǎn)平面應(yīng)力單元組成的有限元模型為例，說明本發(fā)明所提出方法在非線性有限元計(jì)算中的實(shí)施方式。對于平面應(yīng)力單元，單元內(nèi)任意一點(diǎn)的應(yīng)變可由3個分量表示，即ε＝(εxxεyyγxy)^t，其中εxx、εyy為正應(yīng)變，γxy為剪應(yīng)變，相應(yīng)的應(yīng)力分量為σ＝(σxxσyyτxy)^t，其中σxx、σyy為正應(yīng)力，τxy為剪應(yīng)力。此時材料的應(yīng)變分解表達(dá)式為：

式中，ε′＝(ε′xxε′yyγ′xy)^t代表材料線彈性應(yīng)變；ε″＝(ε″xxε″yyγ″xy)^t代表材料非線性應(yīng)變。

在每個單元中設(shè)置一個非線性變形插值點(diǎn)，如圖3所示，則單元非線性應(yīng)變場可表示為如下插值形式：

ε″＝cε″p(2)

式中，c為非線性插值函數(shù)矩陣；ε″p為非線性變形插值點(diǎn)處的非線性應(yīng)變。本例中，由于僅有一個插值點(diǎn)，非線性應(yīng)變在單元內(nèi)按常數(shù)分布，c為單位矩陣。

非線性求解采用newton-raphson迭代方法，在非線性求解前，首先計(jì)算各個單元彈性剛度矩陣ke及交互矩陣k′，對于任意單元，相應(yīng)的ke和k′的計(jì)算公式為：

式中，b為單元應(yīng)變矩陣，代表的單元應(yīng)變分布和結(jié)點(diǎn)位移的關(guān)系；de為材料初始彈性本構(gòu)矩陣。其次，對單元矩陣ke進(jìn)行集成，獲得結(jié)構(gòu)整體彈性剛度矩陣ke。

在非線性求解過程中，任一增量迭代步中的求解方程更新過程可按如下步驟進(jìn)行：

步驟一：進(jìn)行單元狀態(tài)的確定，用于計(jì)算在當(dāng)前位移狀態(tài)下各單元的實(shí)際非線性狀態(tài)，包括計(jì)算非線性插值點(diǎn)處的切線材料本構(gòu)矩陣dt、判斷單元是否進(jìn)入非線性變形狀態(tài)。若進(jìn)入非線性變形狀態(tài)，即dt≠de，則計(jì)算單元的非線性剛度矩陣k″p，本例中，采用高斯積分算法并令非線性變形插值點(diǎn)與高斯積分點(diǎn)重合，則k″p可表示為

k″p＝atde(de-dt)^-1de(5)

其中，a為單元面積；t為單元厚度。上式中de與dt均為非線性變形插值點(diǎn)處的材料屬性。若塑性插值點(diǎn)未進(jìn)入塑性變形狀態(tài)，即dt＝de，則不計(jì)算矩陣k″p。

步驟二：根據(jù)由步驟一所得的各單元的非線性狀態(tài)，更新矩陣k″pr和k′r。其中矩陣k″pr由所有進(jìn)入非線性變形狀態(tài)的單元矩陣k″p集成得到，矩陣k′r由所有進(jìn)入非線性變形狀態(tài)的單元矩陣k′集成得到。

步驟三：更新非線性求解方程并完成求解：