基于面积坐标和解析试函数的薄板元

王丽+赵玉明+龙志飞+魏佳+卢玉林

作者简介： 王丽（1982—），女，辽宁沈阳人，讲师，博士，研究方向为计算力学，(Email)wangl@cidp.edu.cn0引言

基于第二类四边形面积坐标的解析试函数法在文献［1］中被提出.本文把文献［1］中提出的10项试函数延伸到12项，再将薄板理论的基本解析解作为试函数［23］，构造薄板元AATFBQ4.算例计算结果表明所构造的单元具有很高的计算精度且收敛性可靠，是高性能优质单元.

1构造单元

薄板理论采用Kirchhoff直法线假设，法线转角Ψx和Ψy与中面倾角w/x和w/y相等,忽略板的横向剪应变γzx和γyz.薄板的物理控制方程包括平衡方程和本构方程.设载荷集度q, mx和my为0，则平衡方程为Qx=Mxx+Mxyy

Qy=Mxyx+Myy

Qxx+Qyy=0 (1)本构方程为 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)设单元的挠度场和转角场为w=Fλ

Ψx=xFλ

Ψy=yFλ (4)式中：λ含12个参数λ=［λ1λ2λ3…λ12］T利用第二类四边形面积坐标［47］并满足双调和方程（3），设F为含12个系数的多项式，

F=1Z1Z2Z21Z1Z2Z22Z31Z21Z2

Z1Z22Z32f3Z41－f1Z31Z2f3Z42－f2Z32Z1(5)

式中：Z1和Z2均为第二类四边形面积坐标的坐标分量；f1=b12+c12;f2=b22+c22;f3=b1b2+c1c2.其中，b1，b2，c1和c2为第二类四边形面积坐标与直角坐标的转换系数.［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

利用已推导的薄板理论基本解析解作为试函数，采用4个节点的位移w,Ψx和Ψy作为基本未知量，由点协调条件构造一个4节点薄板元AATFBQ4.

首先，选取由式（6)和(7)含12个参数的位移场和曲率场构成单元的物理场子空间.

其次，单元采用12个点协调条件［810］，即wi=wi

Ψxi=Ψxi

Ψyi=Ψyi（8）协调条件为q=Tλ（9）式中：q为节点位移向量，q=［w1Ψx1Ψy1w2Ψx2Ψy2w3Ψx3Ψy3w4Ψx4Ψy4］T （10）对式(9)求逆得λ=T－1q（11）曲率场式(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）应变场ε及其矩阵B已经求出，因此单元刚度矩阵为Ke=∫1－1∫1－1BTDBJdξdη （17）式中：J为坐标变换的雅可比行列式；D为弹性矩阵.单元刚度矩阵可用高斯积分求得.该单元记为AATFBQ4.

2数值算例

2.1算例1

Razzaque斜板［11］问题典型4×4网格见图1，斜板受均布载荷q=1 N/m2作用，边AB和CD简支（w=0，SS1），边CB和DA自由板的材料和几何参数为：E=1 000 Pa,μ=0.31,L=100 mm,h=0.1 m,O为板中心.板中心挠度wo和弯矩My的计算结果见表1，单元AATFBQ4的解精度较高.

图 1Razzaque斜板问题

Fig.1Razzaque skew plate issue

表 1在均布载荷作用下Razzaque斜板中心的挠度和弯矩

Tab.1Central deflection and central bending moment of Razzaque skew plate under uniform load网格数中心挠度

wo/mm中心挠度

相对误差/%中心弯矩

My/(N/m)中心弯矩

相对误差/%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

计算承受均布载荷的简支方板和固支方板，边长为1，泊松比为0.3，均布载荷q=1 N/m2.采用图2所示矩形单元网格计算中心挠度和中心弯矩.由表2和3可看出计算结果随网格加密而收敛于精确解,表中网格数为1/4板的网格划分,d表示板厚.

图 2矩形单元网格

Fig.2Rectangular element mesh表 2在均布载荷作用下四边固支方板中心挠度和中心弯矩

Tab.2Central deflection and central bending moment of square plate clamped at four edges under uniform load网格数中心挠度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在均布载荷作用下四边简支方板中心挠度和中心弯矩

Tab.3Central deflection and central moment of square plate simply supported at four edges under uniform load网格数中心挠度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结束语

由两个算例可以看出，采用面积坐标和基于解析试函数的薄板元AATFBQ4具有较高的精度和较好的收敛性，同时还有较强的稳定性：算例1中计算结果具有较高的精度和稳定性，算例2中随着网格的加密，计算结果很快向精确解收敛.由于单元的构造采用广义协调条件，因此具有良好性能并能通过分片检验，是可靠的薄板元.

参考文献：

［1］陈晓明, 岑松. 基于四边形面积坐标的平面单元解析试函数法［J］. 清华大学学报: 自然科学版, 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］龙驭球, 傅向荣. 基于解析试函数的广义协调四边形厚板元［J］. 工程力学, 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］傅向荣, 龙驭球. 基于解析试函数的广义协调四边形膜元［J］. 工程力学, 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］陈晓明, 岑松, 龙驭球, 等. 含两个分量的四边形单元面积坐标理论［J］. 工程力学, 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.(下转第77页)第23卷 第3期2014年6月计 算 机 辅 助 工 程Computer Aided EngineeringVol.23 No.3Jun. 2014

计算机辅助工程2014年第3期冯敏，等：人工髋骨接触状态仿真文章编号：1006-[KG*9〗0871(2014)03[KG*9〗0069[KG*9〗04

DOI：10.13340/j.cae.2014.03.014

人工髋骨接触状态仿真

冯敏a，从曙光b，郑百林b

(同济大学 a. 体育教学部；b. 航空航天与力学学院， 上海200092)

摘要：利用有限元法对髋臼内衬在一个步态周期载荷下的应力状况进行仿真，探讨髋臼内衬在不同股骨头髋臼材料组合和不同接触位置下的接触状态.结果表明：髋臼内衬的最大接触压力集中在髋臼内衬后1/4部分，其大小与髋骨所受合力一致.髋臼内衬应力受其接触位置影响很大.在某些状态下，不同的股骨头髋臼材料组合往往会对其接触状态产生根本性影响.

关键词：人工髋骨； 髋臼； 内衬； 接触压力； 有限元