基于APDL的ANSYS嵌入式配筋實現

杜永峰+張玉星+朱前坤

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

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

中圖分類號： TU311； P315.9文獻標志碼： B

0引言

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

1計算方法

1.1等參數單元形函數

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

單元坐標插值函數為f1=ni=1Nixi-x=0

f2=ni=1Niyi-y=0

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

v=ni=1Nivi

w=ni=1Niwi（3）通過求解式（2）中被約束節點在自然坐標系中的s，t和r坐標值，應用式（3）求解形函數在約束節點的值，將所求值作為所要建立約束方程中各節點位移的分項系數，實現嵌入節點與被嵌入單元位移的協調.

1.2NewtonRaphson方法

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

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

1.4APDL實現步驟

實現過程應用到ANSYS程序內部的函數命令，其中ENEARN是獲得節點臨近單元號的內部函數命令，*MOPER為求解線性方程組的命令.

基本步驟如下：

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

（2）應用*GET命令獲取嵌入節點的數量，并建立嵌入節點號數組.

（3）應用*DO循環和內部函數命令ENEARN得到每個嵌入節點臨近的嵌入單元號，并存入數組.

（4）獲得嵌入單元節點與被嵌入單元節點在物理坐標系下的坐標值，應用*IF命令判斷單元形函數所屬類型，選擇相應的形函數.

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

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

實現流程見圖1.

圖 1實現流程

Fig.1Implementation process

NR迭代APDL過程，應用*MOPER計算線性方程組，需要將雅可比矩陣和不平衡項FFZERO的第一項寫入命令相應位置.[5]可以設定循環求解的次數，超出求解迭代次數認為不收斂，對本問題可以取值為3，觀察最終的ITERTIME（I），即每個嵌入節點最終的求解迭代次數，結果顯示均為1，由此可以說明NR方法的2階收斂速度.

*DO，M，1，3

*MOPER，RES（1），JACOBI（1，1），SOLV，FFZERO（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計算模型和參數

為驗證計算方法的合理性，將應用Marc計算的結果與本文方法應用ANSYS計算的結果進行比較，通過設定具有相同的幾何模型和物理參數的懸臂柱，對柱頂進行位移加載，比較柱頂的反力位移曲線.

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

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

Fig.2Concrete cantilever column size， mm

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

Fig.3Finite element model

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

圖 4混凝土單軸受壓本構關系

Fig.4Concrete uniaxial compression constitutive relation

2.2計算結果分析

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

2.2.1第一、二組計算對比

第一、二組計算結果為不考慮或考慮鋼筋嵌入計算結果.通過將柱頂的節點進行耦合，對保留節點施加水平方向40 mm位移進行計算，結果見圖5.

圖 5第一組與第二組結果對比

Fig.5Result comparison of group 1 and group 2

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

2.2.2第三組計算對比

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

圖 6第三組結果對比

Fig.6Result comparison of group 3

3結束語

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

[1]王新敏. ANSYS工程結構數值分析[M]. 北京： 人民交通出版社， 2007： 479498.

[2]王新敏， 李義強， 許宏偉. ANSYS結構分析單元與應用[M]. 北京： 人民交通出版社， 2011： 187195.

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

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

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

[6]GB 50010—2010混凝土結構設計規范[S].

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

[8]陸新征， 江見鯨. 用ANSYS SOLID65單元分析混凝土組合構件復雜應力[J]. 建筑結構， 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）設定初值0，求解s，t和r值代入位移插值函數形成各位移的分項系數，應用CE命令建立節點約束方程.

實現流程見圖1.

圖 1實現流程

Fig.1Implementation process

NR迭代APDL過程，應用*MOPER計算線性方程組，需要將雅可比矩陣和不平衡項FFZERO的第一項寫入命令相應位置.[5]可以設定循環求解的次數，超出求解迭代次數認為不收斂，對本問題可以取值為3，觀察最終的ITERTIME（I），即每個嵌入節點最終的求解迭代次數，結果顯示均為1，由此可以說明NR方法的2階收斂速度.

*DO，M，1，3

*MOPER，RES（1），JACOBI（1，1），SOLV，FFZERO（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計算模型和參數

為驗證計算方法的合理性，將應用Marc計算的結果與本文方法應用ANSYS計算的結果進行比較，通過設定具有相同的幾何模型和物理參數的懸臂柱，對柱頂進行位移加載，比較柱頂的反力位移曲線.

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

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

Fig.2Concrete cantilever column size， mm

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

Fig.3Finite element model

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

圖 4混凝土單軸受壓本構關系

Fig.4Concrete uniaxial compression constitutive relation

2.2計算結果分析

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

2.2.1第一、二組計算對比

第一、二組計算結果為不考慮或考慮鋼筋嵌入計算結果.通過將柱頂的節點進行耦合，對保留節點施加水平方向40 mm位移進行計算，結果見圖5.

圖 5第一組與第二組結果對比

Fig.5Result comparison of group 1 and group 2

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

2.2.2第三組計算對比

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

圖 6第三組結果對比

Fig.6Result comparison of group 3

3結束語

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

[1]王新敏. ANSYS工程結構數值分析[M]. 北京： 人民交通出版社， 2007： 479498.

[2]王新敏， 李義強， 許宏偉. ANSYS結構分析單元與應用[M]. 北京： 人民交通出版社， 2011： 187195.

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

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

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

[6]GB 50010—2010混凝土結構設計規范[S].

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

[8]陸新征， 江見鯨. 用ANSYS SOLID65單元分析混凝土組合構件復雜應力[J]. 建筑結構， 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）設定初值0，求解s，t和r值代入位移插值函數形成各位移的分項系數，應用CE命令建立節點約束方程.

實現流程見圖1.

圖 1實現流程

Fig.1Implementation process

NR迭代APDL過程，應用*MOPER計算線性方程組，需要將雅可比矩陣和不平衡項FFZERO的第一項寫入命令相應位置.[5]可以設定循環求解的次數，超出求解迭代次數認為不收斂，對本問題可以取值為3，觀察最終的ITERTIME（I），即每個嵌入節點最終的求解迭代次數，結果顯示均為1，由此可以說明NR方法的2階收斂速度.

*DO，M，1，3

*MOPER，RES（1），JACOBI（1，1），SOLV，FFZERO（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計算模型和參數

為驗證計算方法的合理性，將應用Marc計算的結果與本文方法應用ANSYS計算的結果進行比較，通過設定具有相同的幾何模型和物理參數的懸臂柱，對柱頂進行位移加載，比較柱頂的反力位移曲線.

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

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

Fig.2Concrete cantilever column size， mm

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

Fig.3Finite element model

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

圖 4混凝土單軸受壓本構關系

Fig.4Concrete uniaxial compression constitutive relation

2.2計算結果分析

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

2.2.1第一、二組計算對比

第一、二組計算結果為不考慮或考慮鋼筋嵌入計算結果.通過將柱頂的節點進行耦合，對保留節點施加水平方向40 mm位移進行計算，結果見圖5.

圖 5第一組與第二組結果對比

Fig.5Result comparison of group 1 and group 2

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

2.2.2第三組計算對比

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

圖 6第三組結果對比

Fig.6Result comparison of group 3

3結束語

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

[1]王新敏. ANSYS工程結構數值分析[M]. 北京： 人民交通出版社， 2007： 479498.

[2]王新敏， 李義強， 許宏偉. ANSYS結構分析單元與應用[M]. 北京： 人民交通出版社， 2011： 187195.

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

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

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

[6]GB 50010—2010混凝土結構設計規范[S].

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

[8]陸新征， 江見鯨. 用ANSYS SOLID65單元分析混凝土組合構件復雜應力[J]. 建筑結構， 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