[力学]关于张量的旋转变换
在弹性力学入门时候, 我们接触的主要是各向同性的本构模型. 各向同性就是说, 材料的力学性质在任意方向都相同. 我们可以证明, 对于各向同性材料, 其四阶张量的模量, 只依赖于两个独立参数(三维情形, 下同). 这组参数可以是杨氏模量和泊松比, 也可以是两个拉梅常数等.
然而, 对于更普遍的弹性材料, 各向同性的模型是远远不够的, 我们需要各向异性的模型, 也就是力学性质依赖于材料坐标架. 不仅仅是弹性中存在各向异性, 塑性, 损伤等非线性各向异性本构模型也越来越重要, 尤其在复合材料的力学分析中尤为普遍. 常见的各向异性有横观各向同性 (5个独立参数), 正交各向异性(9个独立参数), 以及各向异性(21个独立参数).
在开发本构模型中, 我们会接触到不同坐标架之间张量的旋转变换问题. 普遍来说, 有限元问题是定义在全局坐标架中, 而针对每个材料积分点, 其材料坐标架可以各不相同. 因此, 对于应变张量, 我们需要从全局坐标架旋转到材料坐标架; 而对于本构模型得到的应力和模量, 我们则需要反其道而行之.
下面就简单整理一下关于二阶和四阶张量旋转的一些问题.
首先, 我们需要定义旋转矩阵
对于旋转矩阵, 我们要求这是一个正交变换, 也就是得满足
其中I 是单位矩阵. 如果
是全局坐标架, 
是局部坐标架的话, 那么
---
对于二阶张量的旋转, 我们有
或者
这里
和
分别是全局和局部坐标架下的一个二阶张量.
从Fortran代码实现的角度说, 只需要
mp = MATMUL(rot,MATMUL(m0,TRANSPOSE(rot)))其中m0和mp分别是全局和局部坐标下的张量, 而rot是旋转矩阵. 这个算法对于三维问题来说, 需要54次浮点乘法计算.
需要提醒的是, 在力学范畴里面, 小变形下的应力应变, Cauchy应力, 应变率, PK2应力, Green应变等都是可以证明是对称矩阵, 因此, 其存储是采用所谓的Voigt记号, 例如:
或者
因此, 在处理旋转应变或者应力的时候, 如果要采用矩阵形式, 不要忘记应变的剪向是工程应变, 对应于两倍的张量应变. 而事实上, 不少有限元程序并不采用之前的矩阵相乘形式, 而是从Voigt记号出发, 推导
具体的每一项的表达式, 然后用Voigt记号改写成矩阵和向量的点乘的形式. 例如对于应变, 我们有
这里:
对于应力, 我们有
这里:
尽管这中算法只需要36次浮点乘法计算, 但计算矩阵
或者
时候, 是需要额外不少浮点乘法计算的. 因此, 后面这个算法适合大量积分点的材料坐标架相同的情况, 因为此时矩阵或者
只需要计算一次.
---
对于四阶张量的旋转, 我们有
或者
这里
和
分别是全局和局部坐标架下的一个二阶张量.
我们需要注意的是, 在本构模型的程序的编写中, 我们并不会将模量用一个四阶数组
REAL(KIND=RKIND) :: D(3,3,3,3)来存储材料的模量, 而是采用Voigt记号, 将这个张量用
REAL(KIND=RKIND) :: D(6,6)来表示.
在说程序之前, 我们先回顾一些力学上的知识点. 我们知道, 三维四阶张量拥有81个独立变量, 而显然, 在Voigt记号下, 这个张量仅仅只有36个独立变量. 所以, 我们得先弄清楚为什么可以这么改写. 原因其实很朴素, 对于四阶张量, 我们称满足
为major symmetry, 而满足
为minor symmetry.
前者来自于应变能密度对应变导数的可交换性, 后者来自于角动量守恒(注: 并非所有材料都满足major symmetry, 但对于小变形下的应力应变, Cauchy应力, 应变率, PK2应力, Green应变, minor symmetry是一定满足的). 对于满足minor symmetry的模量, 我们就可以将其改写成D(6,6)的Voigt记号, 也就是可以从81个独立变量削减成36个. 进而, 若major symmetry也满足, 那么D(6,6)即是对称的, 从而可以进一步削减成21个独立变量.
我们不妨考虑最通用的代码实现, 也就是模量矩阵D(6,6)并不是对称的. 我们有如下两种代码实现
mp = ZERO DO i = 1, 3DO j = i, 3 II = tensor_index_get_symCondense(i,j, 3) DO k = 1, 3 DO l = k, 3 JJ = tensor_index_get_symCondense(k,l, 3) DO m = 1, 3 DO n = 1, 3 KK = tensor_index_get_symCondense(m,n, 3) DO p = 1, 3 DO q = 1, 3 LL = tensor_index_get_symCondense(p,q, 3) mp(II,JJ) = mp(II,JJ) + rot(i,m)*rot(j,n)*rot(k,p)*rot(l,q)*m0(KK,LL) ENDDO ENDDO ENDDO ENDDO ENDDO ENDDOENDDOENDDO 以及
vp = MATMUL(vrot,MATMUL(v0,TRANSPOSE(vrot)))其中vrot是之前提及的, 函数tensor_index_get_symCondense用来计算下标对{i,j}对应的Voigt下标值II:
II = i, 如果 i == j
II = 9-i-j, 如果 i /= j.
我们可以估计计算复杂度, 第一个算法大约需要11664次浮点乘法运算以及大约2916次的赋值, 而第二个算法只需要432次浮点乘法运算以及36次(存疑)赋值. 显然后者要经济许多. 从实践角度来看, 后者的确要快不少. 在我这里的环境下, 单核, 各重复100000次旋转操作, 前者需要2.28s的CPU时间, 后者是0.031s.
注意: 1) 这两个算法并不要求D(6,6)是对称的; 2) 第一个算法中, 内层4个循环需要从1遍历到3, 这个可以从张量旋转的定义中取得; 3) 除去Voigt记号外, 还有Mandel记号, 也是常用的记法; 4) 尽管不常见, 但是要记得, 柔度张量不是简单的求D(6,6)的逆.
相关文章，在仿真秀官网搜索：
[力学]从张量应变和工程应变开始谈起
[力学]关于刚度和模量
无穷魅力的力学系列（四）——任重而道远
无穷魅力的力学系列（三）——魅力四射，五彩缤纷
无穷魅力的力学系列（二）——众说纷纭的力学
无穷魅力的力学系列（一）无处不在的力学
