基于APDL的ANSYS嵌入式配筋實(shí)現(xiàn)

杜永峰+張玉星+朱前坤

摘要： 以ANSYS為平臺(tái)，應(yīng)用APDL實(shí)現(xiàn)鋼筋混凝土有限元模型中鋼筋單元嵌入混凝土單元的功能.應(yīng)用該方法可以在ANSYS中方便地建立復(fù)雜的鋼筋混凝土分離式模型.利用NewtonRaphson方法求解等參數(shù)單元的坐標(biāo)插值函數(shù)中被約束節(jié)點(diǎn)在自然坐標(biāo)系下的坐標(biāo)值，得到自然坐標(biāo)系下節(jié)點(diǎn)對(duì)應(yīng)形函數(shù)值，并將對(duì)應(yīng)的形函數(shù)值作為約束方程中混凝土單元節(jié)點(diǎn)位移的分項(xiàng)因數(shù)，從而實(shí)現(xiàn)嵌入節(jié)點(diǎn)與被嵌入單元位移協(xié)調(diào).計(jì)算結(jié)果與Marc的算例結(jié)果吻合良好，說(shuō)明該方法合理.

關(guān)鍵詞： 鋼筋混凝土； 嵌入式配筋； 分離式模型； NewtonRaphson方法； 約束方程； ANSYS； APDL

中圖分類(lèi)號(hào)： TU311； P315.9文獻(xiàn)標(biāo)志碼： B

0引言

通用有限元程序ANSYS和Marc等可以比較方便地實(shí)現(xiàn)嵌入式配筋.ANSYS在分析鋼筋混凝土構(gòu)件時(shí)有2種方法：一種為彌散式配筋，另一種為分離式配筋.彌散式配筋的優(yōu)點(diǎn)為建立模型簡(jiǎn)單：將鋼筋本構(gòu)矩陣D乘以體積配筋率，疊加在混凝土應(yīng)力、應(yīng)變矩陣上，形成彌散單元的本構(gòu)矩陣.這種方法需要計(jì)算單元的體積配筋率，如果單元?jiǎng)澐植灰?guī)則或配筋復(fù)雜，不容易對(duì)單元設(shè)定相應(yīng)的配筋率，并且在后處理可視化中不容易觀察鋼筋的受力情況.分離式配筋在ANSYS中實(shí)現(xiàn)有3種方法，分別為實(shí)體切分法、節(jié)點(diǎn)耦合法和約束方程法.[1]這3種方法的實(shí)現(xiàn)需要建立的混凝土單元節(jié)點(diǎn)與鋼筋節(jié)點(diǎn)在空間中的距離較為接近：實(shí)體切分需要完全共用節(jié)點(diǎn)，節(jié)點(diǎn)耦合是將距離最近的鋼筋節(jié)點(diǎn)與混凝土節(jié)點(diǎn)耦合在一起，約束方程將鋼筋節(jié)點(diǎn)位移用距離最近的單元面上的節(jié)點(diǎn)位移通過(guò)建立線性約束方程表示.這些方法在ANSYS中的實(shí)現(xiàn)都要考慮鋼筋的位置，否則建立的模型計(jì)算結(jié)果偏差較大.本文利用有限元等參數(shù)單元的特點(diǎn)，應(yīng)用NewtonRaphson方法求解被約束節(jié)點(diǎn)在自然坐標(biāo)系下的坐標(biāo)值，進(jìn)而求得形函數(shù)在該點(diǎn)的值，并將所求值作為約束方程中的節(jié)點(diǎn)位移分項(xiàng)系數(shù)，以此實(shí)現(xiàn)嵌入節(jié)點(diǎn)與被嵌入單元的位移協(xié)調(diào)，實(shí)現(xiàn)單元嵌入的功能.

1計(jì)算方法

1.1等參數(shù)單元形函數(shù)

在有限元法中，一般用等參數(shù)單元對(duì)所要分析的模型進(jìn)行離散計(jì)算，將物理坐標(biāo)系下單元的坐標(biāo)值和位移值用相同的形函數(shù)插值表示，形函數(shù)由自然坐標(biāo)系下的坐標(biāo)值表示.對(duì)于八節(jié)點(diǎn)六面體單元，其形函數(shù)基本形式[2]為Ni=18（1±s）（1±t）（1±r）（1）形函數(shù)通過(guò)物理坐標(biāo)系下的單元節(jié)點(diǎn)坐標(biāo)和單元節(jié)點(diǎn)位移，建立自然坐標(biāo)系中標(biāo)準(zhǔn)單元與物理坐標(biāo)系中實(shí)際劃分單元之間的映射關(guān)系，這種映射關(guān)系一一對(duì)應(yīng)，插值函數(shù)的雅克比矩陣的行列式值在求解域內(nèi)不能為0，保證求解非線性方程組坐標(biāo)值的唯一性.等參元的優(yōu)點(diǎn)是保證整體坐標(biāo)系中不規(guī)則網(wǎng)格邊界位移的連續(xù)性，并通過(guò)映射關(guān)系將物理坐標(biāo)系中的不定積分域轉(zhuǎn)化為自然坐標(biāo)系中的定域積分.[3]

單元坐標(biāo)插值函數(shù)為f1=ni=1Nixi-x=0

f2=ni=1Niyi-y=0

f3=ni=1Nizi-z=0 （2）單元位移插值函數(shù)為u=ni=1Niui

v=ni=1Nivi

w=ni=1Niwi（3）通過(guò)求解式（2）中被約束節(jié)點(diǎn)在自然坐標(biāo)系中的s，t和r坐標(biāo)值，應(yīng)用式（3）求解形函數(shù)在約束節(jié)點(diǎn)的值，將所求值作為所要建立約束方程中各節(jié)點(diǎn)位移的分項(xiàng)系數(shù)，實(shí)現(xiàn)嵌入節(jié)點(diǎn)與被嵌入單元位移的協(xié)調(diào).

1.2NewtonRaphson方法

NewtonRaphson方法是求解非線性方程的方法，該方法通過(guò)迭代求解線性化的方程解決非線性方程的求解問(wèn)題.[4]其基本步驟為：將非線性方程略去高階項(xiàng)并展開(kāi)為一階線性方程，確定求解的初始值，在初始值基礎(chǔ)上迭代求解不平衡項(xiàng)，當(dāng)計(jì)算增量在設(shè)定許可范圍時(shí)即可停止迭代；對(duì)非線性方程組，需要求解方程組的1階偏導(dǎo)數(shù)矩陣即雅可比矩陣.設(shè)定初值為自然坐標(biāo)下的坐標(biāo)原點(diǎn)，八節(jié)點(diǎn)六面體單元的雅可比矩陣為2次的，同時(shí)NewtonRaphson方法具有2階收斂速度[4]，所以對(duì)選定的形函數(shù)只需迭代1次即可獲得足夠精度的解.通過(guò)求解結(jié)果向量的2范數(shù)判斷收斂.該方法的迭代格式為 f1sf1tf1r

f2sf2tf2r

f3sf3tf3rnΔs

Δt

Δr=-f1（sn）

f2（tn）

f3（rn），

sn+1=sn+Δs

tn+1=tn+Δt

rn+1=rn+Δr （4） 1.3約束方程

通過(guò)約束方程將被約束節(jié)點(diǎn)自由度用約束節(jié)點(diǎn)自由度線性表示，求解過(guò)程將整體剛度矩陣中的被約束節(jié)點(diǎn)自由度凝聚，只保留其剛度和節(jié)點(diǎn)載荷的影響，約束方程的基本形式為Cons t=ni=1（Coefficient（i）·u（i）） （5）Coefficient（i）為節(jié)點(diǎn)位移u（i）的分項(xiàng)系數(shù)，分項(xiàng)系數(shù)為被約束節(jié)點(diǎn)在自然坐標(biāo)系下形函數(shù)的值.[1]

1.4APDL實(shí)現(xiàn)步驟

實(shí)現(xiàn)過(guò)程應(yīng)用到ANSYS程序內(nèi)部的函數(shù)命令，其中ENEARN是獲得節(jié)點(diǎn)臨近單元號(hào)的內(nèi)部函數(shù)命令，*MOPER為求解線性方程組的命令.

基本步驟如下：

（1）首先選擇嵌入單元和被嵌入單元，應(yīng)用CM命令建立集合.

（2）應(yīng)用*GET命令獲取嵌入節(jié)點(diǎn)的數(shù)量，并建立嵌入節(jié)點(diǎn)號(hào)數(shù)組.

（3）應(yīng)用*DO循環(huán)和內(nèi)部函數(shù)命令ENEARN得到每個(gè)嵌入節(jié)點(diǎn)臨近的嵌入單元號(hào)，并存入數(shù)組.

（4）獲得嵌入單元節(jié)點(diǎn)與被嵌入單元節(jié)點(diǎn)在物理坐標(biāo)系下的坐標(biāo)值，應(yīng)用*IF命令判斷單元形函數(shù)所屬類(lèi)型，選擇相應(yīng)的形函數(shù).

（5）將被嵌入單元節(jié)點(diǎn)與嵌入單元節(jié)點(diǎn)坐標(biāo)值代入雅克比矩陣和位移插值函數(shù)形成的向量，應(yīng)用*MOPER，*DO和*IF命令語(yǔ)句建立NewtonRaphson求解過(guò)程.[5]endprint

（6）設(shè)定初值0，求解s，t和r值代入位移插值函數(shù)形成各位移的分項(xiàng)系數(shù)，應(yīng)用CE命令建立節(jié)點(diǎn)約束方程.

實(shí)現(xiàn)流程見(jiàn)圖1.

圖 1實(shí)現(xiàn)流程

Fig.1Implementation process

NR迭代APDL過(guò)程，應(yīng)用*MOPER計(jì)算線性方程組，需要將雅可比矩陣和不平衡項(xiàng)FFZERO的第一項(xiàng)寫(xiě)入命令相應(yīng)位置.[5]可以設(shè)定循環(huán)求解的次數(shù)，超出求解迭代次數(shù)認(rèn)為不收斂，對(duì)本問(wèn)題可以取值為3，觀察最終的ITERTIME（I），即每個(gè)嵌入節(jié)點(diǎn)最終的求解迭代次數(shù)，結(jié)果顯示均為1，由此可以說(shuō)明NR方法的2階收斂速度.

*DO，M，1，3

*MOPER，RES（1），JACOBI（1，1），SOLV，F(xiàn)FZERO（1）

S（I）=S（I）+RES（1）

T（I）=T（I）+RES（2）

R（I）=R（I）+RES（3）

*IF，RES（1）**2+RES（2）**2+RES（3）**2，LE，1E6，THEN

*EXIT

*ELSE

*ENDIF

ITERTIME（I）=ITERTIME（I）+1

*ENDDO

2算例分析

2.1計(jì)算模型和參數(shù)

為驗(yàn)證計(jì)算方法的合理性，將應(yīng)用Marc計(jì)算的結(jié)果與本文方法應(yīng)用ANSYS計(jì)算的結(jié)果進(jìn)行比較，通過(guò)設(shè)定具有相同的幾何模型和物理參數(shù)的懸臂柱，對(duì)柱頂進(jìn)行位移加載，比較柱頂?shù)姆戳ξ灰魄€.

柱截面取為300 mm×400 mm，保護(hù)層厚度取為20 mm，混凝土采用我國(guó)《混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范》[6]規(guī)定的混凝土軸心受壓應(yīng)力應(yīng)變關(guān)系，混凝土材料強(qiáng)度取C30，鋼筋為Q235鋼材，箍筋間距為100 mm，鋼筋混凝土構(gòu)件尺寸見(jiàn)圖2.有限元模型見(jiàn)圖3.

圖 2混凝土懸臂柱尺寸， mm

Fig.2Concrete cantilever column size， mm

（a） Marc模型（b） ANSYS模型圖 3有限元模型

Fig.3Finite element model

在有限元計(jì)算中，Marc需要輸入減去彈性應(yīng)變的等效塑性應(yīng)力應(yīng)變關(guān)系.[7]對(duì)于小應(yīng)變問(wèn)題，工程應(yīng)變與真實(shí)應(yīng)變的計(jì)算結(jié)果基本一致.本文輸入為工程應(yīng)力應(yīng)變關(guān)系，輸入的本構(gòu)關(guān)系見(jiàn)圖4.設(shè)定隨動(dòng)強(qiáng)化和等向強(qiáng)化材料屬性，對(duì)單向加載中單元無(wú)卸載的情況計(jì)算結(jié)果沒(méi)有影響.本文設(shè)定為隨動(dòng)強(qiáng)化模型，屈服準(zhǔn)則均設(shè)定為von Mises屈服準(zhǔn)則，鋼筋采用理想彈塑性模型.

圖 4混凝土單軸受壓本構(gòu)關(guān)系

Fig.4Concrete uniaxial compression constitutive relation

2.2計(jì)算結(jié)果分析

分3組模型進(jìn)行計(jì)算：第一組為不嵌入鋼筋的素混凝土，本構(gòu)關(guān)系設(shè)定為不考慮開(kāi)裂的彈塑性本構(gòu)模型；第二組為考慮嵌入鋼筋不設(shè)定開(kāi)裂影響的計(jì)算模型；第三組為考慮嵌入鋼筋和考慮開(kāi)裂的計(jì)算模型.ANSYS混凝土為SOLID65單元，鋼筋為L(zhǎng)INK8單元，Marc實(shí)體選用八節(jié)點(diǎn)7號(hào)單元，桿單元為二節(jié)點(diǎn)9號(hào)單元.

2.2.1第一、二組計(jì)算對(duì)比

第一、二組計(jì)算結(jié)果為不考慮或考慮鋼筋嵌入計(jì)算結(jié)果.通過(guò)將柱頂?shù)墓?jié)點(diǎn)進(jìn)行耦合，對(duì)保留節(jié)點(diǎn)施加水平方向40 mm位移進(jìn)行計(jì)算，結(jié)果見(jiàn)圖5.

圖 5第一組與第二組結(jié)果對(duì)比

Fig.5Result comparison of group 1 and group 2

對(duì)沒(méi)有嵌入鋼筋的素混凝土，2個(gè)程序的計(jì)算結(jié)果相同；對(duì)嵌入鋼筋的計(jì)算結(jié)果，ANSYS的峰值點(diǎn)略高于Marc的計(jì)算結(jié)果，在峰值之前兩者的計(jì)算結(jié)果相同.由此可以說(shuō)明，通過(guò)本文的方法可以實(shí)現(xiàn)在ANSYS程序中嵌入配筋的計(jì)算.

2.2.2第三組計(jì)算對(duì)比

為進(jìn)一步計(jì)算考慮混凝土加載過(guò)程受拉開(kāi)裂的結(jié)果，設(shè)定Marc開(kāi)裂的極限拉應(yīng)力為2 MPa，軟化模量為彈性模量的1/10，極限壓應(yīng)變?yōu)?.003 3，開(kāi)裂單元的剪力傳遞因數(shù)為0.3.ANSYS通過(guò)設(shè)定破壞準(zhǔn)則設(shè)定相同的開(kāi)裂拉應(yīng)力，張開(kāi)裂縫剪力傳遞因數(shù)取0.3，閉合裂縫剪力傳遞因數(shù)取0.95.關(guān)閉壓碎，同時(shí)設(shè)定SOLID65單元的單元選項(xiàng)為1和7，不考慮位移形函數(shù)附加項(xiàng)和考慮開(kāi)裂后的拉應(yīng)力釋放，并將收斂準(zhǔn)則CNVTOL設(shè)定為位移判定，保證計(jì)算收斂.[8]第三組結(jié)果對(duì)比見(jiàn)圖6.在沒(méi)有開(kāi)裂前的彈性段，2個(gè)程序的計(jì)算結(jié)果相同；在開(kāi)裂之后，計(jì)算的反力ANSYS要小于Marc，計(jì)算的結(jié)果趨勢(shì)相近.同時(shí)比較開(kāi)裂與不考慮開(kāi)裂的剛度和峰值，結(jié)果相差較大，開(kāi)裂峰值反力在20 kN左右，不開(kāi)裂峰值反力在80 kN左右.

圖 6第三組結(jié)果對(duì)比

Fig.6Result comparison of group 3

3結(jié)束語(yǔ)

通過(guò)有限元理論，應(yīng)用等參數(shù)單元的坐標(biāo)插值函數(shù)，利用NR方法求解自然坐標(biāo)系中的節(jié)點(diǎn)坐標(biāo)值，并將求解結(jié)果代入位移插值函數(shù)，應(yīng)用線性約束方程將嵌入單元節(jié)點(diǎn)的位移用被嵌入單元節(jié)點(diǎn)的位移表示，應(yīng)用ANSYS實(shí)現(xiàn)實(shí)體單元嵌入桿單元的功能.與Marc的計(jì)算結(jié)果進(jìn)行對(duì)比，結(jié)果表明：在不考慮混凝土開(kāi)裂情況下，兩者計(jì)算結(jié)果吻合良好，驗(yàn)證本文方法的合理性.2種軟件對(duì)開(kāi)裂的處理有一定的差別，考慮開(kāi)裂后計(jì)算結(jié)果雖然趨勢(shì)相近，但計(jì)算結(jié)果仍有一定差距.參考文獻(xiàn)：

[1]王新敏. ANSYS工程結(jié)構(gòu)數(shù)值分析[M]. 北京： 人民交通出版社， 2007： 479498.

[2]王新敏， 李義強(qiáng)， 許宏偉. ANSYS結(jié)構(gòu)分析單元與應(yīng)用[M]. 北京： 人民交通出版社， 2011： 187195.

[3]李人憲. 有限元法基礎(chǔ)[M]. 2版. 北京： 國(guó)防工業(yè)出版社， 2004.

[4]殷有泉. 非線性有限元基礎(chǔ)[M]. 北京： 北京大學(xué)出版社， 2007： 18.

[5]博弈創(chuàng)作室. APDL參數(shù)化有限元分析技術(shù)及其應(yīng)用實(shí)例[M]. 北京： 中國(guó)水利水電出版社， 2004： 9094.

[6]GB 50010—2010混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范[S].

[7]陸新征， 葉列平， 繆志偉， 等. 建筑抗震彈塑性分析[M]. 北京： 中國(guó)建筑工業(yè)出版社， 2009： 169202.

[8]陸新征， 江見(jiàn)鯨. 用ANSYS SOLID65單元分析混凝土組合構(gòu)件復(fù)雜應(yīng)力[J]. 建筑結(jié)構(gòu)， 2003， 33（6）： 2224.

LU Xinzheng， JIANG Jianjing. ANSYS SOLID65 element analysis in concrete composite components with complex stress state[J]. Building Struct， 2003， 33（6）： 2224.

（編輯武曉英）endprint

（6）設(shè)定初值0，求解s，t和r值代入位移插值函數(shù)形成各位移的分項(xiàng)系數(shù)，應(yīng)用CE命令建立節(jié)點(diǎn)約束方程.

實(shí)現(xiàn)流程見(jiàn)圖1.

圖 1實(shí)現(xiàn)流程

Fig.1Implementation process

NR迭代APDL過(guò)程，應(yīng)用*MOPER計(jì)算線性方程組，需要將雅可比矩陣和不平衡項(xiàng)FFZERO的第一項(xiàng)寫(xiě)入命令相應(yīng)位置.[5]可以設(shè)定循環(huán)求解的次數(shù)，超出求解迭代次數(shù)認(rèn)為不收斂，對(duì)本問(wèn)題可以取值為3，觀察最終的ITERTIME（I），即每個(gè)嵌入節(jié)點(diǎn)最終的求解迭代次數(shù)，結(jié)果顯示均為1，由此可以說(shuō)明NR方法的2階收斂速度.

*DO，M，1，3

*MOPER，RES（1），JACOBI（1，1），SOLV，F(xiàn)FZERO（1）

S（I）=S（I）+RES（1）

T（I）=T（I）+RES（2）

R（I）=R（I）+RES（3）

*IF，RES（1）**2+RES（2）**2+RES（3）**2，LE，1E6，THEN

*EXIT

*ELSE

*ENDIF

ITERTIME（I）=ITERTIME（I）+1

*ENDDO

2算例分析

2.1計(jì)算模型和參數(shù)

為驗(yàn)證計(jì)算方法的合理性，將應(yīng)用Marc計(jì)算的結(jié)果與本文方法應(yīng)用ANSYS計(jì)算的結(jié)果進(jìn)行比較，通過(guò)設(shè)定具有相同的幾何模型和物理參數(shù)的懸臂柱，對(duì)柱頂進(jìn)行位移加載，比較柱頂?shù)姆戳ξ灰魄€.

柱截面取為300 mm×400 mm，保護(hù)層厚度取為20 mm，混凝土采用我國(guó)《混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范》[6]規(guī)定的混凝土軸心受壓應(yīng)力應(yīng)變關(guān)系，混凝土材料強(qiáng)度取C30，鋼筋為Q235鋼材，箍筋間距為100 mm，鋼筋混凝土構(gòu)件尺寸見(jiàn)圖2.有限元模型見(jiàn)圖3.

圖 2混凝土懸臂柱尺寸， mm

Fig.2Concrete cantilever column size， mm

（a） Marc模型（b） ANSYS模型圖 3有限元模型

Fig.3Finite element model

在有限元計(jì)算中，Marc需要輸入減去彈性應(yīng)變的等效塑性應(yīng)力應(yīng)變關(guān)系.[7]對(duì)于小應(yīng)變問(wèn)題，工程應(yīng)變與真實(shí)應(yīng)變的計(jì)算結(jié)果基本一致.本文輸入為工程應(yīng)力應(yīng)變關(guān)系，輸入的本構(gòu)關(guān)系見(jiàn)圖4.設(shè)定隨動(dòng)強(qiáng)化和等向強(qiáng)化材料屬性，對(duì)單向加載中單元無(wú)卸載的情況計(jì)算結(jié)果沒(méi)有影響.本文設(shè)定為隨動(dòng)強(qiáng)化模型，屈服準(zhǔn)則均設(shè)定為von Mises屈服準(zhǔn)則，鋼筋采用理想彈塑性模型.

圖 4混凝土單軸受壓本構(gòu)關(guān)系

Fig.4Concrete uniaxial compression constitutive relation

2.2計(jì)算結(jié)果分析

分3組模型進(jìn)行計(jì)算：第一組為不嵌入鋼筋的素混凝土，本構(gòu)關(guān)系設(shè)定為不考慮開(kāi)裂的彈塑性本構(gòu)模型；第二組為考慮嵌入鋼筋不設(shè)定開(kāi)裂影響的計(jì)算模型；第三組為考慮嵌入鋼筋和考慮開(kāi)裂的計(jì)算模型.ANSYS混凝土為SOLID65單元，鋼筋為L(zhǎng)INK8單元，Marc實(shí)體選用八節(jié)點(diǎn)7號(hào)單元，桿單元為二節(jié)點(diǎn)9號(hào)單元.

2.2.1第一、二組計(jì)算對(duì)比

第一、二組計(jì)算結(jié)果為不考慮或考慮鋼筋嵌入計(jì)算結(jié)果.通過(guò)將柱頂?shù)墓?jié)點(diǎn)進(jìn)行耦合，對(duì)保留節(jié)點(diǎn)施加水平方向40 mm位移進(jìn)行計(jì)算，結(jié)果見(jiàn)圖5.

圖 5第一組與第二組結(jié)果對(duì)比

Fig.5Result comparison of group 1 and group 2

對(duì)沒(méi)有嵌入鋼筋的素混凝土，2個(gè)程序的計(jì)算結(jié)果相同；對(duì)嵌入鋼筋的計(jì)算結(jié)果，ANSYS的峰值點(diǎn)略高于Marc的計(jì)算結(jié)果，在峰值之前兩者的計(jì)算結(jié)果相同.由此可以說(shuō)明，通過(guò)本文的方法可以實(shí)現(xiàn)在ANSYS程序中嵌入配筋的計(jì)算.

2.2.2第三組計(jì)算對(duì)比

為進(jìn)一步計(jì)算考慮混凝土加載過(guò)程受拉開(kāi)裂的結(jié)果，設(shè)定Marc開(kāi)裂的極限拉應(yīng)力為2 MPa，軟化模量為彈性模量的1/10，極限壓應(yīng)變?yōu)?.003 3，開(kāi)裂單元的剪力傳遞因數(shù)為0.3.ANSYS通過(guò)設(shè)定破壞準(zhǔn)則設(shè)定相同的開(kāi)裂拉應(yīng)力，張開(kāi)裂縫剪力傳遞因數(shù)取0.3，閉合裂縫剪力傳遞因數(shù)取0.95.關(guān)閉壓碎，同時(shí)設(shè)定SOLID65單元的單元選項(xiàng)為1和7，不考慮位移形函數(shù)附加項(xiàng)和考慮開(kāi)裂后的拉應(yīng)力釋放，并將收斂準(zhǔn)則CNVTOL設(shè)定為位移判定，保證計(jì)算收斂.[8]第三組結(jié)果對(duì)比見(jiàn)圖6.在沒(méi)有開(kāi)裂前的彈性段，2個(gè)程序的計(jì)算結(jié)果相同；在開(kāi)裂之后，計(jì)算的反力ANSYS要小于Marc，計(jì)算的結(jié)果趨勢(shì)相近.同時(shí)比較開(kāi)裂與不考慮開(kāi)裂的剛度和峰值，結(jié)果相差較大，開(kāi)裂峰值反力在20 kN左右，不開(kāi)裂峰值反力在80 kN左右.

圖 6第三組結(jié)果對(duì)比

Fig.6Result comparison of group 3

3結(jié)束語(yǔ)

通過(guò)有限元理論，應(yīng)用等參數(shù)單元的坐標(biāo)插值函數(shù)，利用NR方法求解自然坐標(biāo)系中的節(jié)點(diǎn)坐標(biāo)值，并將求解結(jié)果代入位移插值函數(shù)，應(yīng)用線性約束方程將嵌入單元節(jié)點(diǎn)的位移用被嵌入單元節(jié)點(diǎn)的位移表示，應(yīng)用ANSYS實(shí)現(xiàn)實(shí)體單元嵌入桿單元的功能.與Marc的計(jì)算結(jié)果進(jìn)行對(duì)比，結(jié)果表明：在不考慮混凝土開(kāi)裂情況下，兩者計(jì)算結(jié)果吻合良好，驗(yàn)證本文方法的合理性.2種軟件對(duì)開(kāi)裂的處理有一定的差別，考慮開(kāi)裂后計(jì)算結(jié)果雖然趨勢(shì)相近，但計(jì)算結(jié)果仍有一定差距.參考文獻(xiàn)：

[1]王新敏. ANSYS工程結(jié)構(gòu)數(shù)值分析[M]. 北京： 人民交通出版社， 2007： 479498.

[2]王新敏， 李義強(qiáng)， 許宏偉. ANSYS結(jié)構(gòu)分析單元與應(yīng)用[M]. 北京： 人民交通出版社， 2011： 187195.

[3]李人憲. 有限元法基礎(chǔ)[M]. 2版. 北京： 國(guó)防工業(yè)出版社， 2004.

[4]殷有泉. 非線性有限元基礎(chǔ)[M]. 北京： 北京大學(xué)出版社， 2007： 18.

[5]博弈創(chuàng)作室. APDL參數(shù)化有限元分析技術(shù)及其應(yīng)用實(shí)例[M]. 北京： 中國(guó)水利水電出版社， 2004： 9094.

[6]GB 50010—2010混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范[S].

[7]陸新征， 葉列平， 繆志偉， 等. 建筑抗震彈塑性分析[M]. 北京： 中國(guó)建筑工業(yè)出版社， 2009： 169202.

[8]陸新征， 江見(jiàn)鯨. 用ANSYS SOLID65單元分析混凝土組合構(gòu)件復(fù)雜應(yīng)力[J]. 建筑結(jié)構(gòu)， 2003， 33（6）： 2224.

LU Xinzheng， JIANG Jianjing. ANSYS SOLID65 element analysis in concrete composite components with complex stress state[J]. Building Struct， 2003， 33（6）： 2224.

（編輯武曉英）endprint

（6）設(shè)定初值0，求解s，t和r值代入位移插值函數(shù)形成各位移的分項(xiàng)系數(shù)，應(yīng)用CE命令建立節(jié)點(diǎn)約束方程.

實(shí)現(xiàn)流程見(jiàn)圖1.

圖 1實(shí)現(xiàn)流程

Fig.1Implementation process

NR迭代APDL過(guò)程，應(yīng)用*MOPER計(jì)算線性方程組，需要將雅可比矩陣和不平衡項(xiàng)FFZERO的第一項(xiàng)寫(xiě)入命令相應(yīng)位置.[5]可以設(shè)定循環(huán)求解的次數(shù)，超出求解迭代次數(shù)認(rèn)為不收斂，對(duì)本問(wèn)題可以取值為3，觀察最終的ITERTIME（I），即每個(gè)嵌入節(jié)點(diǎn)最終的求解迭代次數(shù)，結(jié)果顯示均為1，由此可以說(shuō)明NR方法的2階收斂速度.

*DO，M，1，3

*MOPER，RES（1），JACOBI（1，1），SOLV，F(xiàn)FZERO（1）

S（I）=S（I）+RES（1）

T（I）=T（I）+RES（2）

R（I）=R（I）+RES（3）

*IF，RES（1）**2+RES（2）**2+RES（3）**2，LE，1E6，THEN

*EXIT

*ELSE

*ENDIF

ITERTIME（I）=ITERTIME（I）+1

*ENDDO

2算例分析

2.1計(jì)算模型和參數(shù)

為驗(yàn)證計(jì)算方法的合理性，將應(yīng)用Marc計(jì)算的結(jié)果與本文方法應(yīng)用ANSYS計(jì)算的結(jié)果進(jìn)行比較，通過(guò)設(shè)定具有相同的幾何模型和物理參數(shù)的懸臂柱，對(duì)柱頂進(jìn)行位移加載，比較柱頂?shù)姆戳ξ灰魄€.

柱截面取為300 mm×400 mm，保護(hù)層厚度取為20 mm，混凝土采用我國(guó)《混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范》[6]規(guī)定的混凝土軸心受壓應(yīng)力應(yīng)變關(guān)系，混凝土材料強(qiáng)度取C30，鋼筋為Q235鋼材，箍筋間距為100 mm，鋼筋混凝土構(gòu)件尺寸見(jiàn)圖2.有限元模型見(jiàn)圖3.

圖 2混凝土懸臂柱尺寸， mm

Fig.2Concrete cantilever column size， mm

（a） Marc模型（b） ANSYS模型圖 3有限元模型

Fig.3Finite element model

在有限元計(jì)算中，Marc需要輸入減去彈性應(yīng)變的等效塑性應(yīng)力應(yīng)變關(guān)系.[7]對(duì)于小應(yīng)變問(wèn)題，工程應(yīng)變與真實(shí)應(yīng)變的計(jì)算結(jié)果基本一致.本文輸入為工程應(yīng)力應(yīng)變關(guān)系，輸入的本構(gòu)關(guān)系見(jiàn)圖4.設(shè)定隨動(dòng)強(qiáng)化和等向強(qiáng)化材料屬性，對(duì)單向加載中單元無(wú)卸載的情況計(jì)算結(jié)果沒(méi)有影響.本文設(shè)定為隨動(dòng)強(qiáng)化模型，屈服準(zhǔn)則均設(shè)定為von Mises屈服準(zhǔn)則，鋼筋采用理想彈塑性模型.

圖 4混凝土單軸受壓本構(gòu)關(guān)系

Fig.4Concrete uniaxial compression constitutive relation

2.2計(jì)算結(jié)果分析

分3組模型進(jìn)行計(jì)算：第一組為不嵌入鋼筋的素混凝土，本構(gòu)關(guān)系設(shè)定為不考慮開(kāi)裂的彈塑性本構(gòu)模型；第二組為考慮嵌入鋼筋不設(shè)定開(kāi)裂影響的計(jì)算模型；第三組為考慮嵌入鋼筋和考慮開(kāi)裂的計(jì)算模型.ANSYS混凝土為SOLID65單元，鋼筋為L(zhǎng)INK8單元，Marc實(shí)體選用八節(jié)點(diǎn)7號(hào)單元，桿單元為二節(jié)點(diǎn)9號(hào)單元.

2.2.1第一、二組計(jì)算對(duì)比

第一、二組計(jì)算結(jié)果為不考慮或考慮鋼筋嵌入計(jì)算結(jié)果.通過(guò)將柱頂?shù)墓?jié)點(diǎn)進(jìn)行耦合，對(duì)保留節(jié)點(diǎn)施加水平方向40 mm位移進(jìn)行計(jì)算，結(jié)果見(jiàn)圖5.

圖 5第一組與第二組結(jié)果對(duì)比

Fig.5Result comparison of group 1 and group 2

對(duì)沒(méi)有嵌入鋼筋的素混凝土，2個(gè)程序的計(jì)算結(jié)果相同；對(duì)嵌入鋼筋的計(jì)算結(jié)果，ANSYS的峰值點(diǎn)略高于Marc的計(jì)算結(jié)果，在峰值之前兩者的計(jì)算結(jié)果相同.由此可以說(shuō)明，通過(guò)本文的方法可以實(shí)現(xiàn)在ANSYS程序中嵌入配筋的計(jì)算.

2.2.2第三組計(jì)算對(duì)比

為進(jìn)一步計(jì)算考慮混凝土加載過(guò)程受拉開(kāi)裂的結(jié)果，設(shè)定Marc開(kāi)裂的極限拉應(yīng)力為2 MPa，軟化模量為彈性模量的1/10，極限壓應(yīng)變?yōu)?.003 3，開(kāi)裂單元的剪力傳遞因數(shù)為0.3.ANSYS通過(guò)設(shè)定破壞準(zhǔn)則設(shè)定相同的開(kāi)裂拉應(yīng)力，張開(kāi)裂縫剪力傳遞因數(shù)取0.3，閉合裂縫剪力傳遞因數(shù)取0.95.關(guān)閉壓碎，同時(shí)設(shè)定SOLID65單元的單元選項(xiàng)為1和7，不考慮位移形函數(shù)附加項(xiàng)和考慮開(kāi)裂后的拉應(yīng)力釋放，并將收斂準(zhǔn)則CNVTOL設(shè)定為位移判定，保證計(jì)算收斂.[8]第三組結(jié)果對(duì)比見(jiàn)圖6.在沒(méi)有開(kāi)裂前的彈性段，2個(gè)程序的計(jì)算結(jié)果相同；在開(kāi)裂之后，計(jì)算的反力ANSYS要小于Marc，計(jì)算的結(jié)果趨勢(shì)相近.同時(shí)比較開(kāi)裂與不考慮開(kāi)裂的剛度和峰值，結(jié)果相差較大，開(kāi)裂峰值反力在20 kN左右，不開(kāi)裂峰值反力在80 kN左右.

圖 6第三組結(jié)果對(duì)比

Fig.6Result comparison of group 3

3結(jié)束語(yǔ)

通過(guò)有限元理論，應(yīng)用等參數(shù)單元的坐標(biāo)插值函數(shù)，利用NR方法求解自然坐標(biāo)系中的節(jié)點(diǎn)坐標(biāo)值，并將求解結(jié)果代入位移插值函數(shù)，應(yīng)用線性約束方程將嵌入單元節(jié)點(diǎn)的位移用被嵌入單元節(jié)點(diǎn)的位移表示，應(yīng)用ANSYS實(shí)現(xiàn)實(shí)體單元嵌入桿單元的功能.與Marc的計(jì)算結(jié)果進(jìn)行對(duì)比，結(jié)果表明：在不考慮混凝土開(kāi)裂情況下，兩者計(jì)算結(jié)果吻合良好，驗(yàn)證本文方法的合理性.2種軟件對(duì)開(kāi)裂的處理有一定的差別，考慮開(kāi)裂后計(jì)算結(jié)果雖然趨勢(shì)相近，但計(jì)算結(jié)果仍有一定差距.參考文獻(xiàn)：

[1]王新敏. ANSYS工程結(jié)構(gòu)數(shù)值分析[M]. 北京： 人民交通出版社， 2007： 479498.

[2]王新敏， 李義強(qiáng)， 許宏偉. ANSYS結(jié)構(gòu)分析單元與應(yīng)用[M]. 北京： 人民交通出版社， 2011： 187195.

[3]李人憲. 有限元法基礎(chǔ)[M]. 2版. 北京： 國(guó)防工業(yè)出版社， 2004.

[4]殷有泉. 非線性有限元基礎(chǔ)[M]. 北京： 北京大學(xué)出版社， 2007： 18.

[5]博弈創(chuàng)作室. APDL參數(shù)化有限元分析技術(shù)及其應(yīng)用實(shí)例[M]. 北京： 中國(guó)水利水電出版社， 2004： 9094.

[6]GB 50010—2010混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范[S].

[7]陸新征， 葉列平， 繆志偉， 等. 建筑抗震彈塑性分析[M]. 北京： 中國(guó)建筑工業(yè)出版社， 2009： 169202.

[8]陸新征， 江見(jiàn)鯨. 用ANSYS SOLID65單元分析混凝土組合構(gòu)件復(fù)雜應(yīng)力[J]. 建筑結(jié)構(gòu)， 2003， 33（6）： 2224.

LU Xinzheng， JIANG Jianjing. ANSYS SOLID65 element analysis in concrete composite components with complex stress state[J]. Building Struct， 2003， 33（6）： 2224.

（編輯武曉英）endprint