杜永峰+張玉星+朱前坤

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

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

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

0引言

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

1計算方法

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

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

1.2NewtonRaphson方法

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

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

1.4APDL實現(xiàn)步驟

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

基本步驟如下：

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

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

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

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

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

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

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

圖 1實現(xiàn)流程

Fig.1Implementation process

NR迭代APDL過程，應(yīng)用*MOPER計算線性方程組，需要將雅可比矩陣和不平衡項FFZERO的第一項寫入命令相應(yīng)位置.[5]可以設(shè)定循環(huán)求解的次數(shù)，超出求解迭代次數(shù)認(rèn)為不收斂，對本問題可以取值為3，觀察最終的ITERTIME（I），即每個嵌入節(jié)點最終的求解迭代次數(shù)，結(jié)果顯示均為1，由此可以說明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計算模型和參數(shù)

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

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

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

Fig.2Concrete cantilever column size， mm

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

Fig.3Finite element model

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

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

Fig.4Concrete uniaxial compression constitutive relation

2.2計算結(jié)果分析

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

2.2.1第一、二組計算對比

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

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

Fig.5Result comparison of group 1 and group 2

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

2.2.2第三組計算對比

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

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

Fig.6Result comparison of group 3

3結(jié)束語

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

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

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

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

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

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

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

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

[8]陸新征， 江見鯨. 用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ù)形成各位移的分項系數(shù)，應(yīng)用CE命令建立節(jié)點約束方程.

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

圖 1實現(xiàn)流程

Fig.1Implementation process

NR迭代APDL過程，應(yīng)用*MOPER計算線性方程組，需要將雅可比矩陣和不平衡項FFZERO的第一項寫入命令相應(yīng)位置.[5]可以設(shè)定循環(huán)求解的次數(shù)，超出求解迭代次數(shù)認(rèn)為不收斂，對本問題可以取值為3，觀察最終的ITERTIME（I），即每個嵌入節(jié)點最終的求解迭代次數(shù)，結(jié)果顯示均為1，由此可以說明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計算模型和參數(shù)

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

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

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

Fig.2Concrete cantilever column size， mm

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

Fig.3Finite element model

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

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

Fig.4Concrete uniaxial compression constitutive relation

2.2計算結(jié)果分析

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

2.2.1第一、二組計算對比

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

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

Fig.5Result comparison of group 1 and group 2

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

2.2.2第三組計算對比

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

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

Fig.6Result comparison of group 3

3結(jié)束語

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

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

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

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

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

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

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

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

[8]陸新征， 江見鯨. 用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ù)形成各位移的分項系數(shù)，應(yīng)用CE命令建立節(jié)點約束方程.

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

圖 1實現(xiàn)流程

Fig.1Implementation process

NR迭代APDL過程，應(yīng)用*MOPER計算線性方程組，需要將雅可比矩陣和不平衡項FFZERO的第一項寫入命令相應(yīng)位置.[5]可以設(shè)定循環(huán)求解的次數(shù)，超出求解迭代次數(shù)認(rèn)為不收斂，對本問題可以取值為3，觀察最終的ITERTIME（I），即每個嵌入節(jié)點最終的求解迭代次數(shù)，結(jié)果顯示均為1，由此可以說明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計算模型和參數(shù)

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

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

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

Fig.2Concrete cantilever column size， mm

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

Fig.3Finite element model

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

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

Fig.4Concrete uniaxial compression constitutive relation

2.2計算結(jié)果分析

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

2.2.1第一、二組計算對比

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

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

Fig.5Result comparison of group 1 and group 2

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

2.2.2第三組計算對比

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

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

Fig.6Result comparison of group 3

3結(jié)束語

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

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

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

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

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

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

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

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

[8]陸新征， 江見鯨. 用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