河 北 水 利 电 力 学 院 学 报
JournalofHebeiUniversityof WaterResourcesandElectricEngineering
2021 年3 月第31卷第1期
Mar2021Vol31 No1
文章编号：2096 — 5680(2021)01 — 0075 — 06
孙晓勇1 2 ，宋兴海2，侯娜娜12，付建航2，刘立悦1,2
(1.河北省数据中心相变热管理技术创新中心，河北省沧州市重庆路1号061001；
2.河北水利电力学院土木工程学院，河北省沧州市重庆路1号061001)
摘要：经典弹性力学理论用位移梯度表示无限小变形，不考虑旋转变形，把微元体的旋转视为刚体旋转。含偶应力弹性力学
理论将旋转变形以旋转张量表示，微元体旋转和微元体平动位移同量级，而旋转张量和应变张量同量级，旋转张量与旋转矢
量一一对应，用旋转矢量的梯度表示旋转变形。含偶应力弹性力学理论本构关系包括应力-应变关系和偶应力-曲率张量关 系，用等参变换方法离散单元位移到节点上，从虚功原理出发，增加罚函数项以降低有限元方程对高阶单元的需求，推导了拟 解决三维及二维问题的含偶应力弹性线力学有限元理论，可得三维及二维问题中位移、应力、应变等分布情况，对结构进行力 学评价。
关键词：偶应力；旋转变形；旋转张量；张量分析中图分类号：O343
文献标识码:A DOI ： 10. 16046/j. cnki. issn2096-5680. 2021. 01. 015
1经典线弹性理论与考虑偶应力线弹 性理论
在经典弹塑性力学理论中，物体内任意一点的 应力状态只和应变或应变的历史有关，其基本变量
为位移，对位移求梯度得到应变张量，用位移梯度描
述无限小的变形，然后再由一点的应变张量分析得 到应力张量［1］。含偶应力的线弹性力学理论认为, 物体内任意一点的微元体，除有各个方向的位移外,
还有本身的旋转变形，而这种旋转变形并非单纯的 以旋转角表达，而是用和应变张量一个量级的旋转
张量来表示［］。经典弹性力学理论中微元体也有旋 转，其将旋转视作微小的体积单元的刚体转动，即只
有旋转位移，没有旋转变形。这样处理将带来其他
问题：如果以刚体微小旋转表示单元体旋转，那么相 邻的两个微元体将会控制不同的微小旋转，从而无
法满足旋转自由度的连续性。对位移求梯度，会得
到9个不同的分量，而在经典弹性力学理论中，由位
移求梯度得到的应变张量却只有6个独立的分量,
剩下的3个分量去哪里了？含偶应力的弹性理论认 为，另外的3个分量应由旋转张量提供［］。含偶应 力弹性理论中，基本变量仍然是位移，对位移求梯
度，将位移梯度分解为应变张量与旋转张量之和，其
中应变张量是对称张量，而旋转张量是反对称张 量4。旋转张量与旋转矢量一一对应，对旋转矢量
求梯度得到曲率张量，用曲率张量表示微元体的旋 转变形［］。很显然，曲率张量比应变张量低一个量
级。这就解决了相邻单元旋转自由度不连续的矛 盾，因为相邻微元体的旋转矢量是连续的，所以用旋 转矢量或旋转张量表示的旋转自由度也是连续的。 因此，含偶应力弹性理论是经典弹性理论的一个补
充，微元体变形不再只有平动变形，增加了旋转变 形，平动变形以应变张量来表示，而旋转变形以曲率
张量来表示。 含偶应力弹性理论中增加了偶应力与 曲率张量的本构关系，引入了旋转模量作为材料的
基本参数［］。因为偶应力相对于应力小一个量级，因
此在宏观尺寸下，经典弹性理论不考虑偶应力不会影
响到结果的精确性，而在微尺度下，偶应力所占比重
增加，不考虑偶应力将带来结果上的很大误差，这也 就是为什么微尺度受力构件会有尺寸效应的原因。
2含偶应力线弹性理论
含偶应力弹性理论有悠久的发展历史。早在
com
收稿日期20200325 修回日期20200911
基金项目：河北省教育厅教学改革课题(2018GJJG386)
第一作者简介：孙晓勇(1989-)，男，天津市人，河北水利电力学院助教，研究方向为固体力学。E-mail ： 893941282@qq.
76河北水利电力学院学报2021
1887年,Voigt7]就提出，将物体分解为很多个体积
接近于零但不为零的微小单元，微小单元上作用有
体力偶和面力偶，相邻的微小单元接触面上的力偶
连续。1909年，Cosserat兄弟[]提出了Cosserat理
论，他们把物体分解成刚体颗粒，每个颗粒可以平移
还可以转动。物体的变形就以每个颗粒的平移和转
动来表示。在1960年后，Toupin[],Mindlin和
Tiersten[10'11],Koiter[12]等建立了较完善的偶应力
理论[13]，在Cosserat理论基础上，他们认为，微元体
的变形仍由平移和转动组成，但表示转动的矢量是
由位移矢量得到的，相互不独立。2002年,Yang[14]
在Mindlin基础上，又对偶应力理论进行了改进，从
应变能的角度，建立了由应变张量和旋转矢量梯度
共同组成的应变能密度表达式。Yang的理论简化
了偶应力理论的本构关系，将材料参数定为3个，除
传统的杨氏模量和泊松比外，增加了一个内禀特征
尺寸参量，用以表示旋转变形本构关系，即后来的旋
转模量。
2.1含偶应力弹性力学理论有限元推导
2.1.1基于偶应力理论的几何描述
对位移求梯度，可以表示材料内部的无限小变
形，得到非对称的二阶张量，可以分解为2部分：
u犻=）犻+Q犻
其中，应变张量£和旋转张量Q与位移的关系分别
为
）犻—1（u犻+u犻）
旋转矢量a与旋转张量£2——对应，关系式为
3犻——^Q犻：)犻2—22▽x u
Q犻——3犻：）k
其中，）是置换张量。写为矩阵形式:
~dth 狓1du
狓
d u1
d x3
u 狓1d u2d u2 d x2d x3
狌d u3d u3
狓1d x2d x3_
)11)12)13-0Q12Q13-)12)22)23+—Q120。23
)13)23)33_—Q13Q230_
)100燄0—30燄
0)20+300
00)3_000_
旋转张量Q犻是斜对称张量，其物理意义为微元体体积不变时的旋转变形。引入曲率张量x表示材料中任一点转动变形的程度：
X犻—3犻
曲率张量由旋转矢量求梯度得到，而旋转矢量又是位移旋度的一半，可以得到曲率张量和位移矢量的关系表达式：
x—2▽▽x u
而曲率张量的迹为零，即tr X=0o
根据£：£=0，可知应变张量）和旋转张量Q犻相互独立。平动变形和转动变形以图形表示，如图
Fig.1Infinitesimal deformation is decomposed into
translational deformation and rotational deformation 2.1.2基于偶应力理论的Hamilton原理
对一个纯力学过程建立场方程，变形体外表面受外力珋和外力偶亦作用，体内受体力b和体力偶c 作用，动量定理方程为
J u犞=冷十犞
动量矩定理：
¥[p x X u d犞=[狓X i d S+[z n d S
d t p J s J s
+J p x X b d犞+J p c d犞
动能定理：
J p（2…+犝）犞
+J pb•u d犞+-2J p C•▽X u d犞
式中：x为微元体位置，为速度，为单位质量机械能。
建立虚功原理方程
：
第1期孙晓勇等：用张量分析方法推导含偶应力弹性力学有限元理论77
yp(.b•5u+犮•8co)d V+J(t•8u+m•&o)d S—J p•Sw d V=J(:£+m：S%)d V
在含偶应力弹性理论中，位移是惟一的独立变量，旋转角根据位移求得。但是在进行有限元分析时，将节点位移视作未知量可能会引起节点旋转角的不连续。为了避免因未知旋转角引起的需要建立C1连续高阶元，增加虚功原理的约束条件，引入罚函数项：
J u p u•d V+J(8犼6犼+^x^m■犼)d V+
J8(<p，一)<p犻一CO t)d V一
J(狋犻+5«?m,)d S—J(flo S b j+p S«j!c!)d V=0
式中®是未知的旋转角，a是罚函数因子。
把弹性体离散为六面体单元，每个单元有8个节点，每个节点有6个自由度。将单元内任一点位移写成矩阵列向量形式：
u=[x u狔u狕<f，x(p y甲狕]
每个节点的位移也写成矩阵列向量形式：
d i=Lu-ic u狔u*<p x<p狔g*]T(犻=1〜8)将单元内任一点位移离散：
u=Nd=[犖1N2…N8][d1d2…d8]T 其中，矩阵[N i]=N i l6X6，形函数
N i=t(1+黨))1+)1+ZZ))犻=1〜8)应变表示为矩阵的形式：
£=[X]T
其中
£IZ£xx£狔狔£狕狕£xy£狔狕£xz]
x=X xx X m x*X x x*X x狕X x x”狔X x]因此，应变矩阵由节点位移矩阵表示为
£=L u=LNd=Bd
3x000003
3y
03
3”L=030
3
000
30 L<?3；y3x3”
00
3033000
3；y3x
由q—
«l—q)i+2Q k
i—p i+2u)k犲ijk,，得
p一(Mi=L二L Nd=B d
其中
B a=[B“B a2B a B«4B5B6B a7B8]3X48且
-13N i_13N
0N0
0-
23z23狔
B a=一13N i013N0N0
2dz23狓
13N i13N
000N 23狔一23狓
213基于偶应力理论的本构关系
含偶应力理论中应力表示为中国加入世贸组织
£=[m]T
其中
G^(^xx G狔狔O”G x狔O狔”O x]
o=[m x m”m”m x狔m”m x”m”m”
同样，应变表示为
£=[£x]T
其中
]
£IZ£xz£yy£””£xy£yz£””
x=】x x x狔狔x”x x*x狔x x x*x x”
由本构关系G i=2£+A犐和m i=4巧犻得:
O=犇£
m”]
x”=]
其中
D u
D-
D=
其中"A+2^
LN=[B1B2…B8]A+2^
B i=LN i，L=「300 x
030
003
L u0
_0L一
303-
330
3x3”
033
33x
D q=4犐9X9
A+2“
0
持阳伞的女人2.1.4
将体力和体力偶合写为矩阵列向量的形式：二
基于偶应力理论的有限元求解格式
b犮]T。将面力和面力偶合写为矩阵列向量的形
式：二[m]T。将单元体有限元控制方程写为矩
D
卩
00
78
河 北 水 利 电 力 学 院 学 报
2021
阵形式：M e d + (K e +a K 犲d — P e
其中
M = J p N T N d 犞
K — [ B T DB d V
J 犞
K — J B T B a d V
J 犞
P — J p Nb d V + J Nt d S
对坐标进行等参变换：
x —〉： N i
X [,=〉：犖狔犻，狕—〉：N i Z
i
犻=1 i =1 i =1
其中
N i — 1(1 + ee ) (1 + n ) (1 + ZZ )( — 1 〜8)1 1 1 1 1
由 J 」】丄p Nb j i d e d nZ +J J n tA dd 求得
单元应力后，节点应力由高斯点上的应力外推得到。
2. 1. 5平面问题有限元求解格式
对于平面问题，将六面体单元简化为四边形单 元， 每个单元的4 个节点上均存在3 个自由度。
将单元位移列矩阵表示为u = [u 狏°]T，节
点位移列矩阵表示为d = [u 狏 p ]T ( = 1
〜4),
单元上任意一点位移离散：
u = Nd = [N ] N 2 N N 4][d 1 d 2 d 3
d 4]T
其中形函数矩阵：
[Ni ]
— NJ 3X 3
而
好旺角房屋中介
N i = 1(1 + ee ) )1 + nn') )i = 1 〜4)
平面问题应变矩阵：
将N i 对局部坐标求导得:
应变由节点位移表示：
)=Lu = LNd = Bd
其中
B = LN = LNd = [B 1 B 2 B 3 B 4]
B i = LN i
其中，雅阁比矩阵丿可改写为
[3N i
3x 3y 3z 燄3N i 烌
3N i 烌
3e 3e 3e 3Z
3x 3x
3N i
3x 3y
3z 3N i T
3N i
| 3y =
3tj 3rj
3tj 烅3y
＞犑3y 烍3N i 3x
3y 3z 3N i 3N i 3Z
3Z
3Z _3z ,3富
[-L2 0 -L =
[0 L j
L u —33x
3
3狔
3
33x
8
2
i =1
8
2
i =1
8
2
i =1
3N i
i x i 3e x i 8
2
i =1
3N i
3e y i 8
2
i =1
3N i -----z 3e 3N i
i x i 3rj 8
2
i =1
3N i 3狔8
2
i =1
3N i
3n z 3N i
i x i 3Z i
8
2
i =1
3N i 3Z y i
8
2
i =1
3N i
i z i 3Z
根据
于是，质量矩阵M 、刚度矩阵K 、罚函数矩阵K 以
及外力矩阵P 分别写为
111
M ‘ = ] J P N T N | J l ded n dZ
—1—1—1
111
K = J J B T DB | J | ded n dZ
111
K = J
J
b t b ” | j | ded n dz
33x 33y
申
i
3i p i + £ Q 犼犲
丄丄
^pi
电容式压力传感器I 2 u^2 犲 ij2
所以有
p z — 3z
1 3
2 3；y 1
2 3x
p i
3i — L a u — L a N d — B a d
B a — [B a
1 B a
2
B a 3 B 詔]
其中
L 甲
—
J —
即
u
3 1讥
P
B 'ai
1 3N i 1 3N i
2 37 ― 2 討
M]犻=1 〜4)平面问题应力列矩阵为
珋—[狓 O yy
O xy
犿xz 犿狔z ] T
平面问题应力与应变矩阵关系为
O = D )
其中
第1期孙晓勇等：用张量分析方法推导含偶应力弹性力学有限元理论79
D= A+-
D u0一0Dq_
纽A0-
D u=A A+2“0
002“
D q=4犐2X2
平面问题的体力表示为
b=b x b犮狕]
面力表示为
t=t x t y m”]T
于是，单元有限元控制方程可写为
M e d+（K犲+aK：犲）d=犘犲
经等参变换，其中质量矩阵、刚度矩阵、罚函数矩阵和外力矩阵分别可表示为
11
犕=J p N T N|J|d^dy
—
1—1
11
K=J B T DB|J|d e d,
张庭玉
—1—1
11
K a=J b b”|j|ded n
—1—1
111
犘=J P N T b|j|d e d n+J Nt d n
—1—1—1
3结束语
含偶应力弹性理论是经典弹性理论的进一步发展。文中从几何描述、能量原理、本构方程等方面阐述了考虑偶应力对弹性力学理论产生的影响；按照以位移为基本变量，旋转变形由位移导出，总应变分解为对称和反对称部分，对称部分代表经典弹性理论中的应变，即平动变形，反对称部分代表旋转变形，对虚功原理方程增加了罚函数因子，从而引入约束条件，避免了需要构造C1连续单元求解未知旋转角带来的高阶元求解的困难。运用等参元离散有限元方程，对三维及二维问题的有限元求解过程进行了详细推导，其结果可用于计算平板上偶应力分布及旋转变形分布等问题，对三维及二维问题中结构变形、裂纹的产生以及材料的破坏等均可进行力学评价。
参考文献
[1]王敏中，王炜，武际可.弹性力学教程[M].北京：北京
大学出版社，2002.[2]Tang Q G,Meng J P,Liang J S,et al.Effects of
copperbaseda l oysonthe nucleation and growth of calcium carbonate scale[J].Journal of Alloys&Com-pounds，2010，491(1-2)：242-247
[]M,N,Krivosheina,et al.The relationship between e-lasticconstantsandstructureofshock wavesinazinc singlecrystal[J]AIPConferenceProceedings，2017，
1909(1)：20107-20107
[4]Zhanfang Liu,Shijun Yan,Zhi Fu.Dynamic analysis
on generalizedlinear elastic body subjected tolarge scalerigid rotations[J]Applied Mathematics and
Mechanics(EnglishEdition)，2013，34(8)：1001-1016 []颜世军，刘占芳.修正的偶应力线弹性理论及广义线弹性体的有限元方法固体力学学报,2012(3):52-60
[6]Zhanfang Liuand ZhiFu Scalee f ectsofthestress
symmetryin generalized elasticity[J]International
Journal of Aerospace and Lightweight Structures，
2012，2(4)509-521
[7]Brunswik H.Neuere untersuchungen Uber die sexua-
litatsverhaltnisse bei den pilzen[J].MGG-Molecular andGeneralGenetics，1924，34(1)214-228
[8]Eugene et Fran ois Cosserat.Theorie des corps defor-
mables[J]Nature，1909，81(2072)67
[9]Toupin R A Elastic materials withcouple-stresses
[J].Archive for Rational Mechanics&Analysis,
1962，11(1)385-414
[10]MindlinR D，Tiersten H F E f ectsofcouple-stres-
sesinlinearelasticity[J]ArchiveforRational Me缝隙腐蚀
chanics and Analysis，1962，11(1)415-448
[11]R D Mindlin Influenceofcouple-stressesonstress
concentrations[J]ExperimentalMechanics，1963，3
(1)1-7
[12]Ivan Hlavacek,Miroslav Hlavacek.On the existence
anduniquenessofsolutionandsomevariationalprin-
ciplesinlineartheoriesofelasticitywithcouple-stres-
ses.I：Cosserat continuum[J].Clinical&Experi-
mentalA l ergy，1969，19(2)197-202
[13]Hadjesfandiari A R，Dargush G F Couplestress
theoryforsolids[J]internationaljournalofsolids
&structures,2011,48(18)24962510.
[14]F Yang，A C M Chong，D C C Lam，etal
Couplestressbasedstraingradienttheoryforelastic-
ty[J].International Journal of Solids&Structures,
2002，39(10)：2731-2743
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议