止于至善

0%

全局坐标系与局部坐标系下弹性矩阵的变换

发表于 更新于 分类于

变换矩阵\(~\boldsymbol{Q}\)

将全局坐标系\(~O-\boldsymbol{e}_1\boldsymbol{e}_2\boldsymbol{e}_3~\)与局部坐标系\(~o-\boldsymbol{n}_1\boldsymbol{n}_2\boldsymbol{n}_3\)的单位正交基矢量均记为列阵形式，即 \[ \tag{1a} \boldsymbol{e} = [\boldsymbol{e}_1~,~\boldsymbol{e}_2~,~\boldsymbol{e}_3]^\mathrm{T} \] \[ \tag{1b} \boldsymbol{n} = [\boldsymbol{n}_1~,~\boldsymbol{n}_2~,~\boldsymbol{n}_3]^\mathrm{T} \]

于是，可以通过两个坐标系的正交基矢量来定义变换矩阵\(~\boldsymbol{Q}\)，即 \[ \tag{2a} \boldsymbol{e} = \boldsymbol{Q} \cdot \boldsymbol{n} \] \[ \tag{2b} \boldsymbol{e}_i = Q_{ij} \boldsymbol{n}_{j} \]

根据变换矩阵\(~\boldsymbol{Q}\)，可以定义固体力学中常用的一阶、二阶以及四阶张量的变换格式。

一阶张量的变换

一阶张量\(~\boldsymbol{a}~\)在全局坐标系、局部坐标系下的分量分别为\(~a_i\)、\(\tilde{a}_i\)，将两者构成的列阵分别记为\(~\{a\}\)、\(\{\tilde{a}\}\)，则有 \[ \tag{3a} \{a\} = \boldsymbol{Q} \cdot \{\tilde{a}\} \] \[ \tag{3b} a_i = \tilde{a}_j Q_{ij} \]

二阶张量的变换

二阶张量\(~\boldsymbol{F}~\)在全局坐标系、局部坐标系下的分量分别为\(~F_{ij}\)、\(\tilde{F}_{ij}\)，将两者构成的矩阵分别记为\(~\big[F\big]\)、\(\big[\tilde{F}\big]\)，则有 \[ \tag{4a} \big[F\big] = \boldsymbol{Q} \cdot \big[\tilde{F}\big] \cdot \boldsymbol{Q}^\mathrm{T} \] \[ \tag{4b} F_{ij} = \tilde{F}_{mn} Q_{im} Q_{jn} \]

四阶张量的变换

四阶张量\(~\boldsymbol{\mathcal{C}}~\)在全局坐标系、局部坐标系下的分量分别为\(~\mathcal{C}_{ijkl}\)、\(\tilde{\mathcal{C}}_{ijkl}\)，则有 \[ \tag{5} \mathcal{C}_{ijkl} = \tilde{\mathcal{C}}_{mnpq} Q_{im} Q_{jn} Q_{kp} Q_{lq} \]

弹性矩阵的变换

线弹性体的广义胡克定律如下 \[ \tag{6} \boldsymbol{\sigma} = \boldsymbol{\mathcal{C}} : \boldsymbol{\varepsilon} \]

式中，\(\boldsymbol{\mathcal{C}}~\)为四阶弹性张量，\(\boldsymbol{\sigma}\)、\(\boldsymbol{\varepsilon}~\)分别为二阶应力张量、二阶应变张量。

考虑到应力、应变张量的对称性，可将式\((6)\)简化为 \[ \tag{7} \boldsymbol{\sigma} = \boldsymbol{C} \cdot \boldsymbol{\varepsilon} \]

式中，\(\boldsymbol{C}~\)为二阶弹性张量，\(\boldsymbol{\sigma}\)、\(\boldsymbol{\varepsilon}~\)分别为一阶应力张量、一阶应变张量。

按照ABAQUS中的存储顺序，将一阶应力张量\(~\boldsymbol{\sigma}\)、一阶应变张量\(~\boldsymbol{\varepsilon}~\)记为列阵形式，即 \[ \tag{8a} \{\sigma\} = \big[\sigma_{11}~,~\sigma_{22}~,~\sigma_{33}~,~\sigma_{12}~,~\sigma_{13}~,~\sigma_{23}\big]^\mathrm{T} \] \[ \tag{8b} \{\varepsilon\} = \big[\varepsilon_{11}~,~\varepsilon_{22}~,~\varepsilon_{33}~,~2\varepsilon_{12}~,~2\varepsilon_{13}~,~2\varepsilon_{23}\big]^\mathrm{T} \]

需要注意，列阵\(~\{\varepsilon\}~\)中采用的是工程剪应变，而非张量剪应变

二阶应力张量\(~\boldsymbol{\sigma}~\)在全局坐标系、局部坐标系下的分量分别为\(~\sigma_{ij}\)、\(\tilde{\sigma}_{ij}\)，将两者构成的矩阵分别记为\(~\big[\sigma\big]\)、\(\big[\tilde{\sigma}\big]\)，则有 \[ \tag{9} \big[\sigma\big] = \boldsymbol{Q} \cdot \big[\tilde{\sigma}\big] \cdot \boldsymbol{Q}^\mathrm{T} \]

将正交变换矩阵\(~\boldsymbol{Q}~\)的分量代入式\((9)\)，按照式\((8)\)的简记形式，整理可得 \[ \tag{10a} \{\sigma\} = \boldsymbol{T}_{\sigma} \cdot \{\tilde{\sigma}\} \] \[ \tag{10b} \boldsymbol{T}_{\sigma} = \left[ \begin{array}{cc} ~\boldsymbol{T}_1~ & 2\boldsymbol{T}_2 \\ ~\boldsymbol{T}_3~ & \boldsymbol{T}_4 \end{array}\right]_{6\times6} \]

式中，\(\{\tilde{\sigma}\}~\)与\(~\{\sigma\}~\)的形式一致，且 \[ \begin{aligned} \boldsymbol{T}_1 & = \begin{bmatrix} Q_{11}^2 & Q_{12}^2 & Q_{13}^2 \\ Q_{21}^2 & Q_{22}^2 & Q_{23}^2 \\ Q_{31}^2 & Q_{32}^2 & Q_{33}^2 \end{bmatrix} \\[3pt] \boldsymbol{T}_2 & = \begin{bmatrix} Q_{11}Q_{12} & Q_{11}Q_{13} & Q_{12}Q_{13} \\ Q_{21}Q_{22} & Q_{21}Q_{23} & Q_{22}Q_{23} \\ Q_{31}Q_{32} & Q_{31}Q_{33} & Q_{32}Q_{33} \end{bmatrix} \\[3pt] \boldsymbol{T}_3 & = \begin{bmatrix} Q_{11}Q_{21} & Q_{12}Q_{22} & Q_{13}Q_{23} \\ Q_{11}Q_{31} & Q_{12}Q_{32} & Q_{13}Q_{33} \\ Q_{21}Q_{31} & Q_{22}Q_{32} & Q_{23}Q_{33} \end{bmatrix} \\[3pt] \boldsymbol{T}_4 & = \begin{bmatrix} Q_{11}Q_{22}+Q_{12}Q_{21} & Q_{11}Q_{23}+Q_{13}Q_{21} & Q_{12}Q_{23}+Q_{13}Q_{22} \\ Q_{11}Q_{32}+Q_{12}Q_{31} & Q_{11}Q_{33}+Q_{13}Q_{31} & Q_{12}Q_{33}+Q_{13}Q_{32} \\ Q_{21}Q_{32}+Q_{22}Q_{31} & Q_{21}Q_{33}+Q_{23}Q_{31} & Q_{22}Q_{33}+Q_{23}Q_{32} \end{bmatrix} \end{aligned} \]

  • 《复合材料力学（第2版）》266页也给出了列阵\(~\{\sigma\}~\)、\(\{\tilde{\sigma}\}~\)之间的变换矩阵，不过其并不等于\(~\boldsymbol{T}_{\sigma}\)，因为该书中列阵采用的存储顺序为\(~\big[\sigma_{11}~,~\sigma_{22}~,~\sigma_{33}~,~\sigma_{23}~,~\sigma_{13}~,~\sigma_{12}\big]^\mathrm{T}\)。
  • 将式\((9)\)展开推导\(~\boldsymbol{T}_{\sigma}~\)的具体表达式较为繁琐，建议采用支持符号运算的MATLAB、Mathematicsa等软件来辅助推导。

同样地，对于二阶应变张量\(~\boldsymbol{\varepsilon}\)，考虑其不同坐标系下的分量形式，可得 \[ \tag{11a} \{\varepsilon\} = \boldsymbol{T}_{\varepsilon} \cdot \{\tilde{\varepsilon}\} \] \[ \tag{11b} \boldsymbol{T}_{\varepsilon} = \left[\begin{array}{cc} \boldsymbol{T}_1 & ~\boldsymbol{T}_2~ \\ 2\boldsymbol{T}_3 & ~\boldsymbol{T}_4~ \end{array}\right]_{6\times6} \]

式中，\(\{\tilde{\varepsilon}\}~\)与\(~\{\varepsilon\}~\)的形式一致，且\(~\boldsymbol{T}_1\)、\(\boldsymbol{T}_2\)、\(\boldsymbol{T}_2\)、\(\boldsymbol{T}_4~\)与式\((10\mathrm{b})\)中定义一致。

正是由于式\((8\mathrm{b})\)中列阵\(~\{\varepsilon\}~\)采用了工程剪应变，故\(~\boldsymbol{T}_{\varepsilon} \neq \boldsymbol{T}_{\sigma}\)。

考虑不同坐标系下的胡克定律形式，由式\((7)\)可得 \[ \tag{12a} \{\sigma\} = \boldsymbol{C} \cdot \{\varepsilon\} \] \[ \tag{12b} \{\tilde{\sigma}\} = \tilde{\boldsymbol{C}} \cdot \{\tilde{\varepsilon}\} \]

将式\((10)\)、\((11)\)代入式\((12)\)可得 \[ \tag{13a} \boldsymbol{C} = \boldsymbol{T}_{\sigma} \cdot \tilde{\boldsymbol{C}} \cdot \boldsymbol{T}_{\varepsilon}^{-1} \] \[ \tag{13b} \tilde{\boldsymbol{C}} = \boldsymbol{T}_{\sigma}^{-1} \cdot \boldsymbol{C} \cdot \boldsymbol{T}_{\varepsilon} \]

关于矩阵\(~\boldsymbol{T}_{\sigma}\)、\(\boldsymbol{T}_{\varepsilon}\)，具有如下性质 \[ \tag{14a} \boldsymbol{T}_{\sigma}^\mathrm{T} = \boldsymbol{T}_{\varepsilon}^{-1} \] \[ \tag{14b} \boldsymbol{T}_{\varepsilon}^\mathrm{T} = \boldsymbol{T}_{\sigma}^{-1} \] \[ \tag{14c} \mathrm{det}(\boldsymbol{T}_{\sigma}) = \mathrm{det}(\boldsymbol{T}_{\varepsilon}) = 1 \]

令\(~\boldsymbol{T}=\boldsymbol{T}_{\sigma}\)，则\(~\boldsymbol{T}^\mathrm{T}=\boldsymbol{T}_{\varepsilon}^{-1}\)，由式\((13)\)可得 \[ \tag{15a} \boldsymbol{C} = \boldsymbol{T} \cdot \tilde{\boldsymbol{C}} \cdot \boldsymbol{T}^\mathrm{T} \] \[ \tag{15b} \tilde{\boldsymbol{C}} = \boldsymbol{T}^{-1} \cdot \boldsymbol{C} \cdot \boldsymbol{T}^{-\mathrm{T}} \]

至此，推导完成。

各向同性弹性矩阵的验证

对于各向同性的材料，可用两个\(~\mathrm{Lam} \acute {\mathrm{e}}~\)常数（\(\lambda~\)和\(~\mu\)）来表示其固体力学特性，进一步将式\((7)\)中的弹性矩阵写作 \[ \tag{16} \boldsymbol{C} = \left[ \begin{array}{cc} \begin{matrix} \lambda+2\mu & \lambda & \lambda \\ ~ & \lambda+2\mu & \lambda \\ ~ & ~ & \lambda+2\mu \end{matrix} & \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \\ \begin{matrix} ~ & ~ & ~ \\ ~ & \mathrm{symm} & ~ \\ ~ & ~ & ~ \end{matrix} & \begin{matrix} \mu & 0 & 0 \\ ~ & \mu & 0 \\ ~ & ~ & \mu \end{matrix} \end{array}\right] \]

任意由单位正交基矢量构成的直角坐标系（螺旋方向一致，一般取右手系）下，各向同性弹性矩阵应是一致的，可将式\((16)\)代入式\((15)\)进行证明。

下面给出数值验证的代码，其中mat_local为局部坐标系下的弹性矩阵，Q为变换矩阵，mat_global为全局坐标系下的弹性矩阵，并由函数trans_EM计算得到。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
% version - R2018b
clc;clear;

% Lame constant
lambda=115;
mu=77;

% isotropic elastic matrix in local coordinate system
mat_local=zeros(6);
mat_local(1:3,1:3)=lambda*ones(3)+diag(2*mu*ones(3,1));
mat_local(4:6,4:6)=diag(mu*ones(3,1));

% orthogonal matrix Q
Q=orth(rands(3,3));

mat_global=trans_EM(mat_local,Q);
disp([' err = ',num2str(norm(mat_global-mat_local))]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stress_array=[s11,s22,s33,  s12,  s13,  s23]
% strain_array=[e11,e22,e33,2*e12,2*e13,2*e23] (engineering strain)
% INPUT
%   x    - elastic matrix in local coordinate system,  matrix 6x6
%   rot  - rotation matrix, matrix 3x3
%
% OUTPUT
%   y    - elastic matrix in global coordinate system, matrix 6x6
function y=trans_EM(x,rot)

    mat1=[rot(:,1),rot(:,1:2)];mat2=[rot(:,2:3),rot(:,3)];
    mat3=[rot(1,:);rot(1:2,:)];mat4=[rot(2:3,:);rot(3,:)];

    T1=rot.*rot;
    T2=mat1.*mat2;
    T3=mat3.*mat4;

    % elementary matrix
    e1=[1,0,0;0,0,1;0,1,0];e2=[0,1,0;1,0,0;0,0,1];

    % elementary column transformation
    T4=(mat3*e1).*(mat4*e2)+(mat3*e2).*(mat4*e1);

    % transformation matrix
    T=[T1,2*T2;T3,T4];
    
    % elastic matrix in global coordinate system
    y=T*x*T';

end