彭文明，段云嶺

(1.中國電建集團成都勘測設(shè)計研究院有限公司勘測設(shè)計分公司，成都，610072；2.清華大學水利系，北京，100084；3.水電水利規(guī)劃設(shè)計總院，北京，100120)

1 前言

針對RCCD仿真模擬的問題，朱伯芳[1-2]首次提出了仿真分析的并層算法，把混凝土材料屬性接近的澆筑層合并后統(tǒng)一劃分網(wǎng)格，變成均質(zhì)單元，以節(jié)省單元數(shù)和節(jié)點數(shù)，降低有限元分析中線性代數(shù)方程的個數(shù)。此后，王建江等[3]提出了“非均勻單元方法”；凌道盛等[4]構(gòu)造出“虛擬層合單元”；朱岳明等[5]提出“非均質(zhì)層合單元法”。這些方法在一定程度上解決了規(guī)模太大的問題，但是在計算過程中，單元均質(zhì)化會帶來誤差。

非均質(zhì)層合單元(以下統(tǒng)一簡稱“層合單元”)與常規(guī)均質(zhì)單元采用一樣的插值函數(shù)。當單元內(nèi)各層材料屬性相差較大的時候，這些插值函數(shù)描述的位移場與事實不盡相符。圖1(a)為一個具有3層材料的層合單元，在均布荷載的作用下，變形如(b)中實線所示。這種變形后的位移場，用常規(guī)4節(jié)點等參元的雙線性插值函數(shù)已經(jīng)無法正確描述，如不做修正，將影響計算精度。

(a)初始單元 (b)變形后的單元

針對層合單元中不同層的位移場特點，本文嘗試對常規(guī)層合單元(standard laminated element，簡稱為SLE)的位移插值函數(shù)進行適當修正，以便提高計算精度。修正后的層合單元記為“修正層合單元”(modified laminated element，簡稱MLE)，以示區(qū)別。

2 層合單元及其形函數(shù)

對于如圖2所示層合單元，單元內(nèi)含有n層材料，層厚分別為ti，材料均為各向同性。

圖2 層合單元

層合單元的形函數(shù)與常規(guī)四節(jié)點等參元一致，如下表示：

(1)

式中：Ni(i=1，2，3，4)為形函數(shù)；ξ和η為局部坐標系。

(2)

式中：tj為第j層材料的厚度；ηi-1和ηi分別為第i層材料界面的局部坐標值；n為材料層數(shù)。

所以分層積分方法共取2×n個積分點，權(quán)系數(shù)為：

Hij=HiHj(i=1,2；j=1,2,…n)

(3)

3 層合單元的形函數(shù)修正方法

層合單元的位移沿層厚方向不是線性變化，在層間的交界面上存在轉(zhuǎn)折點。當各層材料屬性相差較大時，單元內(nèi)的位移場不適合用式(1)描述。在下面的修正中，我們將強制位移場滿足層間交界面上的位移，體現(xiàn)各層材料對位移場的影響。

3.1 求解各層材料分界點的形函數(shù)值

為方便，假定泊松比都一樣，設(shè)為ν。對于第i種材料，其彈性模量為aiE0(E0為常量)。假定單元內(nèi)的平均應(yīng)力為σ，每一層內(nèi)應(yīng)變?yōu)棣舏。對于平面應(yīng)變問題有：

(4)

其中σz=ν(σx+σy)。

(5)

νAi可表示為：

νAi=αiν1+βiν4

(6)

其中αi+βi=1，而ν4=ν41+ν1，所以式(6)可寫成

νAi=βiν41+ν1

(7)

(8)

由式(8)可求出N1在Ai的取值分別為αi(i=1，2，…，n-1)，N4在Ai的取值分別為βi(i=1，2，…，n-1)；如果各層材料在單元內(nèi)厚度處處相等，可得N2在βi的取值為αi(i=1，2，…，n-1)，N3在βi的取值則為βi(i=1，2，…，n-1)。

3.2 采用高階曲線修正插值函數(shù)

上節(jié)求解得到各層材料界面點的位移，可知層合單元位移曲線應(yīng)為折線形。顯然，用常規(guī)四節(jié)點等參元的雙線性函數(shù)是無法準確多點折線的位移模式，本文采用高階曲線修正插值函數(shù)。

圖3 層合單元變形示意

設(shè)插值函數(shù)為

(9)

對于Ni來說，fi(ξ，η)應(yīng)滿足：①在邊12、23、34上取值均為0；②在點Ai的取值使得Ni等于ai(i=1，2，…n-1)。

由條件①，可令

f1(ξ，η)=(1-ξ)(1-η2)g(ξ，η)

(10)

由于條件②中ξ為常量，則g(ξ，η)=g(η)，設(shè)為g(η)為η的多項式，如下：

(11)

其中sj為多項式系數(shù)。所以式(10)為：

(12)

由于N1在Ai中的取值已知，共有n-1個等式，最多能求解n-1個未知數(shù)，所以取m=n-2。設(shè)點Ai局部坐標為(-1，ηAi)，由式(8)-式(12)可得

(13)

[β]{S}={C}

(14)

其中：

{S}=[s0s1s2…sn-2]T

{C}=[c1c2c3…cn-1]T

由于ηAi均不相等，B的行向量線性無關(guān)，于是|B|≠0，則有：

{S}=[B]-1{C}

(15)

由式(15)求出式(12)中的未知數(shù){S}，f1(ξ，η)可以唯一確定。同理可求fi(ξ，η)(i=2，3，4)。fi(ξ，η)(i=1，2，3，4)有如下關(guān)系：

3.3 層合單元的積分格式

對于常規(guī)4節(jié)點等參元而言，其單元剛度矩陣[K]按式(16)求解：

(16)

其中：B為應(yīng)變矩陣；D為彈性矩陣；J為雅克比矩陣行列式；h為單元厚度。

(17)

4 懸臂梁例題

圖4為一個懸臂梁模型，尺寸為15m寬、4.5m高、1m厚，為平面應(yīng)變問題。懸臂梁左側(cè)固支，在模型右側(cè)上方加一豎向集中荷載100kN，不計自重。用C語言編制層合單元有限元計算程序，分別采用常規(guī)層合單元和修正層合單元計算。模型網(wǎng)格劃分為30個單元，44個節(jié)點。

圖4 懸臂梁模型

為檢驗修正層合單元函數(shù)的有效性，下面分兩種情況計算：(1)每個單元內(nèi)沿Y向包含5層材料(厚度均等)，由下而上材料的彈性模量分別為50GPa、45GPa、40GPa、35GPa和30GPa；(2)每個單元內(nèi)沿Y向包含10層材料(厚度均等)，由下至上各層材料彈性模量等差遞減：50GPa、47GPa、……、23GPa。材料的泊松比都取0.16。

為比較，還使用MSC.PATRAN &NASTRAN大型有限元商業(yè)軟件對上述例子建模計算，商業(yè)軟件沒有層合單元計算模塊，只能根據(jù)材料特性分別劃分網(wǎng)格，5層材料時150個單元，176個節(jié)點；10層材料時600個單元，651個節(jié)點。單元網(wǎng)格分別是層合單元的5倍、10倍。

4.1 位移計算結(jié)果

圖5-圖6是用上述三種方法計算得到C截面和D截面位移結(jié)果。圖中的P結(jié)果(5)、原結(jié)果(5)和修正結(jié)果(5)分別表示PATRAN、常規(guī)層合單元和修正層合單元的計算結(jié)果，括號中的“5”表示5層材料。

(a)水平位移

(a)水平位移

從上述各圖可以看出，由于層合單元內(nèi)材料屬性差異較大，如不修正插值函數(shù)，計算結(jié)果偏差比較大，嚴重影響計算精度；使用修正層合單元計算分析，精度得到明顯提高。

4.2 應(yīng)力計算結(jié)果

圖7為D截面和E截面的X向應(yīng)力計算結(jié)果。由于施加的外荷載不大，5層材料和10層材料在相同斷面的σx應(yīng)力值不大，其中在梁的上下面差值相對較大。

圖7 D和E截面X向應(yīng)力結(jié)果

4.3 計算成果對比

典型斷面節(jié)點的位移和應(yīng)力對比分析如表1所示。與商業(yè)軟件對各材料劃分網(wǎng)格的計算成果相比，常規(guī)層合單元的計算誤差平均值在5%～6%，修正層合單元的計算誤差平均值控制在1.5%～2.0%，計算精度有較大提高。

表1 計算結(jié)果對比

5 總結(jié)

層合單元對每層材料進行高斯積分，使得單元內(nèi)材料屬性可以不相同，為成層結(jié)構(gòu)的仿真分析提供了一條切實可行的途徑。由于各層材料屬性不同，層合單元沿層厚方向的變形呈折線狀，而它的形函數(shù)為雙線性插值函數(shù)，該函數(shù)只能描述單向線性變化的位移場，所以在有限元分析時，層合單元只是用線性變化的位移場去逼近實際位移場。當材料屬性相差較大時，上述逼近將產(chǎn)生很大誤差。

本文在常規(guī)層合單元的基礎(chǔ)上，根據(jù)其變形的特點，對層合單元的插值函數(shù)進行適當修正，修正后的插值函數(shù)可以表示為連續(xù)函數(shù)，充分考慮了單元內(nèi)各層材料屬性的差異，所以描述的位移場更能反映實際。算例表明，修正層合單元在不增加單元和節(jié)點的情況下，有效地提高了計算結(jié)果的精度，具有一定的工程應(yīng)用價值。