基于面積坐標(biāo)和解析試函數(shù)的薄板元

王麗+趙玉明+龍志飛+魏佳+盧玉林

作者簡(jiǎn)介： 王麗（1982—），女，遼寧沈陽(yáng)人，講師，博士，研究方向?yàn)橛?jì)算力學(xué)，(Email)wangl@cidp.edu.cn0引言

基于第二類四邊形面積坐標(biāo)的解析試函數(shù)法在文獻(xiàn)［1］中被提出.本文把文獻(xiàn)［1］中提出的10項(xiàng)試函數(shù)延伸到12項(xiàng)，再將薄板理論的基本解析解作為試函數(shù)［23］，構(gòu)造薄板元AATFBQ4.算例計(jì)算結(jié)果表明所構(gòu)造的單元具有很高的計(jì)算精度且收斂性可靠，是高性能優(yōu)質(zhì)單元.

1構(gòu)造單元

薄板理論采用Kirchhoff直法線假設(shè)，法線轉(zhuǎn)角Ψx和Ψy與中面傾角w/x和w/y相等,忽略板的橫向剪應(yīng)變?chǔ)脄x和γyz.薄板的物理控制方程包括平衡方程和本構(gòu)方程.設(shè)載荷集度q, mx和my為0，則平衡方程為Qx=Mxx+Mxyy

Qy=Mxyx+Myy

Qxx+Qyy=0 (1)本構(gòu)方程為 Mx=－D2wx2+μ2wy2

My=－D2wy2+μ2wx2

Mxy=－(1－μ)D2wxy (2)式中：D=Eh312(1－μ2)，其中h為板厚，μ為泊松比.

將式（2）代入式（1）得Δ4w=4wx4+24wx2y2+4wy4=0 (3)設(shè)單元的撓度場(chǎng)和轉(zhuǎn)角場(chǎng)為w=Fλ

Ψx=xFλ

Ψy=yFλ (4)式中：λ含12個(gè)參數(shù)λ=［λ1λ2λ3…λ12］T利用第二類四邊形面積坐標(biāo)［47］并滿足雙調(diào)和方程（3），設(shè)F為含12個(gè)系數(shù)的多項(xiàng)式，

F=1Z1Z2Z21Z1Z2Z22Z31Z21Z2

Z1Z22Z32f3Z41－f1Z31Z2f3Z42－f2Z32Z1(5)

式中：Z1和Z2均為第二類四邊形面積坐標(biāo)的坐標(biāo)分量；f1=b12+c12;f2=b22+c22;f3=b1b2+c1c2.其中，b1，b2，c1和c2為第二類四邊形面積坐標(biāo)與直角坐標(biāo)的轉(zhuǎn)換系數(shù).［4］

因此w,Ψx和Ψy可表示為w

Ψx

Ψy= 1Z1Z2Z21Z1Z2Z22Z31Z21Z2Z1Z22

0b1Ab2A2b1AZ1b1Z2+b2Z1A2b2AZ23b1AZ21Λ1Λ3

0c1Ac2A2c1AZ1c1Z2+c2Z1A2c2AZ23c1AZ21Λ2Λ4Z32Λ5Λ8

3b2AZ23Λ6Λ9

3c2AZ22Λ7Λ10λ(6)式中：Λ1=2b1Z1Z2+b2Z21A

Λ2=2c1Z1Z2+c2Z21A

Λ3=2b2Z1Z2+b1Z22A

Λ4=2c2Z1Z2+c1Z22A

Λ5=f3Z41－f1Z31Z2

Λ6=4b1f3Z31－3b1f1Z21Z2－b2f1Z31A2

Λ7=4c1f3Z31－3c1f1Z21Z2－c2f1Z31A2

Λ8=f3Z42－f2Z32Z1

Λ9=4b2f3Z32－3b2f2Z22Z1－b1f2Z32A2

Λ10=4c2f3Z32－3c2f2Z22Z1－c1f2Z32A2

曲率k=kx

ky

2kxy=κ1κ2κ3κ4λ(7)式中：κ1=000

000

000

κ2=2b12A22b1b2A22b22A2

2c12A22c1c2A22c22A2

4b1c1A22b1c2+b2c1A24b2c2A2κ3=6b12Z1A22b12Z2+4b1b2Z1A22b22Z1+4b1b2Z2A26b22Z2A2

6c12Z1A22c12Z2+4c1c2Z1A22c22Z1+4c1c2Z2A26c22Z2A2

12b1c1Z1A24b1c1Z2+b2c1Z1+b1c2Z1A24b2c2Z1+b2c1Z2+b1c2Z2A212b2c2Z2A2

κ4=12b12f3Z21－6b1b2f1Z21－6b12f1Z1Z2A212b22f3Z22－6b1b2f2Z22－6b22f2Z1Z2A2

12c12f3Z21－6c1c2f1Z21－6c12f1Z1Z2A212c22f3Z22－6c1c2f2Z22－6c22f2Z1Z2A2

24b1c1f3Z21－6c1b2f1Z21－6b1c2f1Z21－12b1c1f1Z1Z2A224b2c2f3Z22－6c1b2f2Z22－6b2c2f2Z22－12b2c2f2Z1Z2A2

利用已推導(dǎo)的薄板理論基本解析解作為試函數(shù)，采用4個(gè)節(jié)點(diǎn)的位移w,Ψx和Ψy作為基本未知量，由點(diǎn)協(xié)調(diào)條件構(gòu)造一個(gè)4節(jié)點(diǎn)薄板元AATFBQ4.

首先，選取由式（6)和(7)含12個(gè)參數(shù)的位移場(chǎng)和曲率場(chǎng)構(gòu)成單元的物理場(chǎng)子空間.

其次，單元采用12個(gè)點(diǎn)協(xié)調(diào)條件［810］，即wi=wi

Ψxi=Ψxi

Ψyi=Ψyi（8）協(xié)調(diào)條件為q=Tλ（9）式中：q為節(jié)點(diǎn)位移向量，q=［w1Ψx1Ψy1w2Ψx2Ψy2w3Ψx3Ψy3w4Ψx4Ψy4］T （10）對(duì)式(9)求逆得λ=T－1q（11）曲率場(chǎng)式(7)可記為ε=κ=Sλ=ST－1q=Bq（12）由式(7)給出S.又B=ST－1（13）令=［wΨxΨy］T（14）則式(6)可記為=Fλ=FT－1q=Nq（15）式中：F由式(5)給出.又N=FT－1 （16）應(yīng)變場(chǎng)ε及其矩陣B已經(jīng)求出，因此單元?jiǎng)偠染仃嚍镵e=∫1－1∫1－1BTDBJdξdη （17）式中：J為坐標(biāo)變換的雅可比行列式；D為彈性矩陣.單元?jiǎng)偠染仃嚳捎酶咚狗e分求得.該單元記為AATFBQ4.

2數(shù)值算例

2.1算例1

Razzaque斜板［11］問(wèn)題典型4×4網(wǎng)格見(jiàn)圖1，斜板受均布載荷q=1 N/m2作用，邊AB和CD簡(jiǎn)支（w=0，SS1），邊CB和DA自由板的材料和幾何參數(shù)為：E=1 000 Pa,μ=0.31,L=100 mm,h=0.1 m,O為板中心.板中心撓度wo和彎矩My的計(jì)算結(jié)果見(jiàn)表1，單元AATFBQ4的解精度較高.

圖 1Razzaque斜板問(wèn)題

Fig.1Razzaque skew plate issue

表 1在均布載荷作用下Razzaque斜板中心的撓度和彎矩

Tab.1Central deflection and central bending moment of Razzaque skew plate under uniform load網(wǎng)格數(shù)中心撓度

wo/mm中心撓度

相對(duì)誤差/%中心彎矩

My/(N/m)中心彎矩

相對(duì)誤差/%2×20.759 14.460.112 116.94×40.777 62.140.101 15.438×80.788 00.820.097 51.6716×160.791 90.330.096 40.52注：wo精確解為0.794 5 mm，My精確解為0.095 89 N/m.2.2算例2

計(jì)算承受均布載荷的簡(jiǎn)支方板和固支方板，邊長(zhǎng)為1，泊松比為0.3，均布載荷q=1 N/m2.采用圖2所示矩形單元網(wǎng)格計(jì)算中心撓度和中心彎矩.由表2和3可看出計(jì)算結(jié)果隨網(wǎng)格加密而收斂于精確解,表中網(wǎng)格數(shù)為1/4板的網(wǎng)格劃分,d表示板厚.

圖 2矩形單元網(wǎng)格

Fig.2Rectangular element mesh表 2在均布載荷作用下四邊固支方板中心撓度和中心彎矩

Tab.2Central deflection and central bending moment of square plate clamped at four edges under uniform load網(wǎng)格數(shù)中心撓度wc/mm中心彎矩Mc/(N/m)d=10-30 mmd=10-3 mmd=10-2 mmd=10-2 mmd=10-30 mmd=10-3 mmd=10-2 mmd=10-2 mm2×20.140 30.140 30.140 30.140 30.027 80.027 80.027 80.027 84×40.130 40.130 50.130 40.130 30.024 00.024 10.024 00.023 98×80.127 60.127 50.127 40.127 50.023 10.023 20.023 30.023 216×160.126 80.126 90.126 80.126 70.023 00.023 10.023 00.023 032×320.126 60.126 60.126 60.126 60.022 90.022 90.022 90.023 0精確解［12］0.126 50.126 50.126 50.149 90.022 9050.022 9050.0229 050.023 1

表 3在均布載荷作用下四邊簡(jiǎn)支方板中心撓度和中心彎矩

Tab.3Central deflection and central moment of square plate simply supported at four edges under uniform load網(wǎng)格數(shù)中心撓度wc/mm中心彎矩Mc/(N/m)d=10-30 mmd=10-3 mmd=10-2 mmd=10-1 mmd=10-30 mmd=10-3 mmd=10-2 mmd=10-1 mm2×20.432 80.432 80.432 80.432 80.052 20.052 20.052 20.052 24×40.412 90.412 90.412 90.412 90.048 90.048 90.048 90.048 98×80.407 90.407 90.407 90.407 90.048 10.048 10.048 10.048 116×160.406 70.406 70.406 70.406 70.048 00.048 00.048 00.048 032×320.406 30.406 30.406 30.406 30.047 90.047 90.047 90.047 9精確解［12］0.406 20.406 20.406 40.427 3 0.047 89 0.047 89 0.047 89 0.047 89

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

由兩個(gè)算例可以看出，采用面積坐標(biāo)和基于解析試函數(shù)的薄板元AATFBQ4具有較高的精度和較好的收斂性，同時(shí)還有較強(qiáng)的穩(wěn)定性：算例1中計(jì)算結(jié)果具有較高的精度和穩(wěn)定性，算例2中隨著網(wǎng)格的加密，計(jì)算結(jié)果很快向精確解收斂.由于單元的構(gòu)造采用廣義協(xié)調(diào)條件，因此具有良好性能并能通過(guò)分片檢驗(yàn)，是可靠的薄板元.

參考文獻(xiàn)：

［1］陳曉明, 岑松. 基于四邊形面積坐標(biāo)的平面單元解析試函數(shù)法［J］. 清華大學(xué)學(xué)報(bào): 自然科學(xué)版, 2008, 48(2): 289293.

CHEN Xiaoming, CEN Song. Analytical trial function method for plane elements based on quadrilateral area coordinates theory［J］. Tsinghua Univ: Sci & Tech, 2008, 48(2): 289293.

［2］龍馭球, 傅向榮. 基于解析試函數(shù)的廣義協(xié)調(diào)四邊形厚板元［J］. 工程力學(xué), 2002, 19(3): 1015.

LONG Yuqiu, FU Xiangrong. Two generalized conforming quadrilateral thick plate elements based on analytical trial functions［J］. Eng Mech, 2002, 19(3):1015.

［3］傅向榮, 龍馭球. 基于解析試函數(shù)的廣義協(xié)調(diào)四邊形膜元［J］. 工程力學(xué), 2002, 19(4): 1216.

FU Xiangrong, LONG Yuqiu. Generalized conforming quadrilateral plane elements based on analytical trial functions ［J］. Eng Mech, 2002, 19(4): 1216.

［4］陳曉明, 岑松, 龍馭球, 等. 含兩個(gè)分量的四邊形單元面積坐標(biāo)理論［J］. 工程力學(xué), 2007, 24(S1): 3235.

CHEN Xiaoming, CEN Song, LONG Yuqiu, et al. A twocomponent area coordinate method for quadrilateral elements［J］. Eng Mech, 2007, 24(S1): 3235.(下轉(zhuǎn)第77頁(yè))第23卷 第3期2014年6月計(jì) 算 機(jī) 輔 助 工 程Computer Aided EngineeringVol.23 No.3Jun. 2014

計(jì)算機(jī)輔助工程2014年第3期馮敏，等：人工髖骨接觸狀態(tài)仿真文章編號(hào)：1006-[KG*9〗0871(2014)03[KG*9〗0069[KG*9〗04

DOI：10.13340/j.cae.2014.03.014

人工髖骨接觸狀態(tài)仿真

馮敏a，從曙光b，鄭百林b

(同濟(jì)大學(xué) a. 體育教學(xué)部；b. 航空航天與力學(xué)學(xué)院， 上海200092)

摘要：利用有限元法對(duì)髖臼內(nèi)襯在一個(gè)步態(tài)周期載荷下的應(yīng)力狀況進(jìn)行仿真，探討髖臼內(nèi)襯在不同股骨頭髖臼材料組合和不同接觸位置下的接觸狀態(tài).結(jié)果表明：髖臼內(nèi)襯的最大接觸壓力集中在髖臼內(nèi)襯后1/4部分，其大小與髖骨所受合力一致.髖臼內(nèi)襯應(yīng)力受其接觸位置影響很大.在某些狀態(tài)下，不同的股骨頭髖臼材料組合往往會(huì)對(duì)其接觸狀態(tài)產(chǎn)生根本性影響.

關(guān)鍵詞：人工髖骨； 髖臼； 內(nèi)襯； 接觸壓力； 有限元