等 参 单 元

前言：本文是作者对于等参数元学习的一些理解和看法，有不对的地方欢迎批评指正，对于不同的见解也欢迎交流。

等参变换

等参元主要涉及到坐标变换：将规则的几何形状看做为母单元，之后可由母单元变换成所需要的各种子单元。我们可以假定这种变换为$(x,y)=f(\xi ,\eta )$ （不需要知道变换关系的具体表达，只是为了便于理解假定存在这种变换。)

​ 在母单元构造形函数$N_i^{\prime }(\xi ,\eta )$，由坐标变换可知：在任意子单元中存在形函数$N_i(x,y)$与之相对应。在局部坐标中形函数满足：
$$\xi =\sum _{i=1}^m{N}_i^{\prime }{\xi }_i\eta =\sum _{i=1}^m{N}_i^{\prime }{\eta }_i\text{\hspace{1em}(1)}$$
相应子单元的形函数也满足：
$$x=\sum _{i=1}^m{N}_i{x}_{i}y=\sum _{i=1}^m{N}_i{y}_i\text{\hspace{1em}(2)}$$
  由坐标变换可知：自变量一一对应，因变量值并不发生改变，只是函数关系发生改变而已。

所以式（2）可改写为
$$x=\sum _{i=1}^m{N}_i^{\prime }(\xi ,\eta )x_{i}y=\sum _{i=1}^m{N}_i^{\prime }(\xi ,\eta )y_i\text{\hspace{1em}(3)}$$
（3）其实也就是之前假设变换关系的插值函数的形式。

​   经过以上的推论，得到了坐标变换的插值函数形式，由于母单元形函数是比较容易求得的，所以也就比较轻松的得到了坐标变换关系，便于之后在积分时进行换元和求雅克比矩阵，利于程序化，这也是等参元的优势。

   上述坐标变换的插值形式中所利用到的节点数目为$m$，在表达场函数时，也是利用节点进行插值表达，假设此时利用的节点数目为$n$。如果两者采用相同的节点，且插值函数相同，那么称为等参变换，此时$m=n$。若$m>n$，为超参；$m<n$,为次(亚)参。

应用

平面应变问题中，在整体坐标下，单元矩阵的表达为：
$$k={\int }_{\Omega }B(N_{i,x},N_{i,y})DBdxdy$$
利用换元到局部坐标：
$$k={\iint }_{{\Omega }^{\prime }}B(N_{i,x}^{\prime },N_{i,y}^{\prime })DB|J|d\xi d\eta$$
雅克比矩阵：
$$\begin{gathered} J=\frac{\partial (x,y)}{\partial (\xi ,\eta )}=\left[\begin{array}{cc}{\sum }_{i=1}^m\frac{\partial N_i^{\prime }}{\partial \xi }{x}_i& {\sum }_{i=1}^m\frac{\partial N_i^{\prime }}{\partial \xi }yi\\ {\sum }_{i=1}^m\frac{\partial N_i^{\prime }}{\partial \xi }{x}_i& {\sum }_{i=1}^m\frac{\partial N_i^{\prime }}{\partial \xi }yi\end{array}\right] \\ =\left[\begin{array}{cccc}\frac{\partial N_1^{\prime }}{\partial \xi }& \frac{\partial N_2^{\prime }}{\partial \xi }& \dots & \frac{\partial N_m^{\prime }}{\partial \xi }\\ \frac{\partial N_1^{\prime }}{\partial \eta }& \frac{\partial N_2^{\prime }}{\partial \eta }& \dots & \frac{\partial N_m^{\prime }}{\partial \eta }\end{array}\right]\left[\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\\ ⋮& ⋮\\ x_m& y_m\end{array}\right] \end{gathered}$$
由
$$\left[\begin{array}{c}\frac{\partial N_i^{\prime }}{\partial \xi }\\ \frac{\partial N_i^{\prime }}{\partial \eta }\end{array}\right]=\left[\begin{array}{cc}\frac{\partial x}{\partial \xi }& \frac{\partial y}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& \frac{\partial x}{\partial \eta }\end{array}\right]\left[\begin{array}{c}\frac{\partial N_i^{\prime }}{\partial x}\\ \frac{\partial N_i^{\prime }}{\partial y}\end{array}\right]=J\left[\begin{array}{c}\frac{\partial N_i^{\prime }}{\partial x}\\ \frac{\partial N_i^{\prime }}{\partial y}\end{array}\right]$$
得
$$\left[\begin{array}{c}\frac{\partial N_i^{\prime }}{\partial x}\\ \frac{\partial N_i^{\prime }}{\partial y}\end{array}\right]=J^{-1}\left[\begin{array}{c}\frac{\partial N_i^{\prime }}{\partial \xi }\\ \frac{\partial N_i^{\prime }}{\partial \eta }\end{array}\right]$$

采用Lagrange单元，所以四节点单元的形函数表达为：
$$N^{\prime }=\frac{1}{4}(1+{\xi }_i\xi )(1+{\eta }_i\eta )$$
形函数的偏导
$$\left[\begin{array}{c}\frac{\partial N^{\prime }}{\partial \xi }\\ \frac{\partial N^{\prime }}{\partial \eta }\end{array}\right]=\frac{1}{4}\left[\begin{array}{c}{\xi }_i(1+{\eta }_i\eta )\\ (1+{\xi }_i\xi ){\eta }_i\end{array}\right]$$

• 四边形四节点单元应用等参元计算刚度矩阵的Matlab程序

function K = Stiffness(gElement,gNode,E,nu)

% 本例每个方向上取4个高斯点
Guass = [-0.861136311594052575224  0.347854845137453857373
         -0.339981043584856264803  0.652145154862546142627
         0.339981043584856264803  0.652145154862546142627
         0.861136311594052575224  0.347854845137453857373];

D = E/(1-nu^2)*[1 nu 0; nu 1 0; 0 0 (1-nu)/2]; % 弹性矩阵

K = zeros(2*size(gNode,1));               % 刚度矩阵预定义
for iel = 1 : size(gElement,1)      % 对单元进行循环
    Coor_SU = gNode(gElement(iel,:),:);    % 整体坐标       [x,y]
    Coor_PU = [-1,-1;1,-1;1,1;-1,1];            %  局部坐标     【xi,eta】
    
    k = zeros(8);
    for i = 1:4
        for  j = 1 : 4
            % 雅克比矩阵
            dNdxe = zeros(2,4);   % 用于储存高斯点处，形函数的偏导值   第一行： 对xi的偏导   第二行： 对eta的偏导
            for m = 1 : 4
                dNdxe(1,m) = Coor_PU(m,1)*(1+Coor_PU(m,2)*Guass(j,1)) /4;      
                dNdxe(2,m) = (1+Coor_PU(m,1)*Guass(i,1))*Coor_PU(m,2) /4;
            end
            Jacob = dNdxe*Coor_SU;
            
            dNdxy = Jacob\dNdxe;  % 插值函数对x y的偏导矩阵
            
            % 应变矩阵
            B = zeros(3,8);
            for  k = 1 : 4
                B(1,2*k-1) = dNdxy(1,k);B(2,2*k) = dNdxy(2,k);
                B(3,2*k-1) = dNdxy(2,k);B(3,2*k) = dNdxy(1,k);
            end
            
            k = k + B'*D*B*det(Jacob)*Guass(i,2)*Guass(j,2);
        end
    end
    
    % 组装刚度矩阵
    for m = 1 : 4
        for n = 1 : 4
            kmn = k(2*m-1:2*m,2*n-1:2*n);
            I = gElement(iel,m); J = gElement(iel,n);
            K(2*I-1:2*I,2*J-1:2*J) = K(2*I-1:2*I,2*J-1:2*J)+ kmn;
        end
    end
end

参考文献：
[王勖成. 有限单元法[M]. 清华大学出版社, 2003.]
[徐荣桥. 结构分析的有限元法与MATLAB程序设计[M]. 人民交通出版社, 2006.]