首页/文章/ 详情

从理论到实践：空间四面体单元的有限元编程与数值实现

易木木响叮当

2月前浏览1328

本期分享的空间四面体单元的有限元编程格式，从理论到数值实现，步步紧跟，面向有限元初级学习者，若对本节理论内容感到疑惑，可先阅览《有限元基础编程百科全书》前面章节内容。

空间四面体的理论知识与平面等参单元理论极为相似，如果上一章节的内容把握到位的话，空间实体单元对于你来说，简直就是洒洒水。

虽然单元类型看起来不一样、空间维度不同，但是在应用有限元理论时、在代码移植过程中，可谓是极其相似。只需要在原始代码中修改关键地方即可，但是本章节还是需要事无巨细地将空间实体单元部分呈现给大家，望读者莫要烦躁，耐心看完，也许会有另一番收获！

主要内容有：

常应变四面体单元 C3D4

直接坐标法
等参变换法
单元介绍
有限元格式
单元刚度矩阵计算

高阶四面体单元 C3D10

单元介绍
有限元格式
单元刚度矩阵计算

数值案例

常应变四面体单元 C3D4

单元介绍

如下图所示的 C3D4 单元，又常常被称为常应变四面体单元，这个叫法可参考常应变三角形单元介绍。


物理坐标系	自然坐标系

每个节点的位移有 3 个分量，、、，可以写为矩阵形式：

单元共有 12 个节点位移分量：

有限元格式

接下来会带着大家再次熟悉有限元基本格式，在以后的有限元问题描述上我都会不断的串讲有限元基本格式，加深印象。因为格式真的很固定，特别适合计算机处理，当然这不代表有限元简单，只不过很多内容不需要我们去推倒，直接拿来用就行（应用型教学）。以下的理论读者不要觉得枯燥，是和编程直接相关的理论基础，本书将与编程无关的理论均剔除，在后面的程序中，以下的理论公式均有所体现。

形函数

摘抄 Abaqus 对于 C3D4 单元位移模式的定义:

即

位移模式

单元内某一点的位移，可通过形函数与节点位移描述。

形函数也可以反映母单元和子单元的映射关系，设母单元内部一点的坐标为，可由形函数插值得到：

雅可比矩阵可表示为：

不要被上面的表达式所吓倒哈！其实在数值编程的过程中，可以把它换成另一种形式：

如此以来是不是简单许多，左侧项为形函数的偏导矩阵，右侧项为节点坐标。

function [JacobianMatrix,XYDerivatives] = Jacobian(nodeCoordinates,naturalDerivatives)
    JacobianMatrix = nodeCoordinates'*naturalDerivatives;
    XYDerivatives = naturalDerivatives/JacobianMatrix;
end

几何方程

单元内某一点的应变，可通过几何矩阵B与节点位移描述。

其中：

本构矩阵

单元刚度矩阵计算

直接坐标法

function k =C3D4Stiffness(prop,elemNodes)
E=prop(1);
NU=prop(2);
x1 = elemNodes(1, 1);
y1 = elemNodes(1, 2);
z1 = elemNodes(1, 3);
x2 = elemNodes(2, 1);
y2 = elemNodes(2, 2);
z2 = elemNodes(2, 3);
x3 = elemNodes(3, 1);
y3 = elemNodes(3, 2);
z3 = elemNodes(3, 3);
x4 = elemNodes(4, 1);
y4 = elemNodes(4, 2);
z4 = elemNodes(4, 3);
xyz = [1 x1 y1 z1 ; 1 x2 y2 z2 ; 1 x3 y3 z3 ; 1 x4 y4 z4];
V = det(xyz)/6;
mbeta1 = [1 y2 z2 ; 1 y3 z3 ; 1 y4 z4];
mbeta2 = [1 y1 z1 ; 1 y3 z3 ; 1 y4 z4];
mbeta3 = [1 y1 z1 ; 1 y2 z2 ; 1 y4 z4];
mbeta4 = [1 y1 z1 ; 1 y2 z2 ; 1 y3 z3];
mgamma1 = [1 x2 z2 ; 1 x3 z3 ; 1 x4 z4];
mgamma2 = [1 x1 z1 ; 1 x3 z3 ; 1 x4 z4];
mgamma3 = [1 x1 z1 ; 1 x2 z2 ; 1 x4 z4];
mgamma4 = [1 x1 z1 ; 1 x2 z2 ; 1 x3 z3];
mdelta1 = [1 x2 y2 ; 1 x3 y3 ; 1 x4 y4];
mdelta2 = [1 x1 y1 ; 1 x3 y3 ; 1 x4 y4];
mdelta3 = [1 x1 y1 ; 1 x2 y2 ; 1 x4 y4];
mdelta4 = [1 x1 y1 ; 1 x2 y2 ; 1 x3 y3];
beta1 = -1*det(mbeta1);
beta2 = det(mbeta2);
beta3 = -1*det(mbeta3);
beta4 = det(mbeta4);
gamma1 = det(mgamma1);
gamma2 = -1*det(mgamma2);
gamma3 = det(mgamma3);
gamma4 = -1*det(mgamma4);
delta1 = -1*det(mdelta1);
delta2 = det(mdelta2);
delta3 = -1*det(mdelta3);
delta4 = det(mdelta4);
B1 = [beta1 0 0 ; 0 gamma1 0 ; 0 0 delta1 ; 
   gamma1 beta1 0 ; 0 delta1 gamma1 ; delta1 0 beta1];
B2 = [beta2 0 0 ; 0 gamma2 0 ; 0 0 delta2 ; 
   gamma2 beta2 0 ; 0 delta2 gamma2 ; delta2 0 beta2];
B3 = [beta3 0 0 ; 0 gamma3 0 ; 0 0 delta3 ; 
   gamma3 beta3 0 ; 0 delta3 gamma3 ; delta3 0 beta3];
B4 = [beta4 0 0 ; 0 gamma4 0 ; 0 0 delta4 ; 
   gamma4 beta4 0 ; 0 delta4 gamma4 ; delta4 0 beta4];
B = [B1 B2 B3 B4]/(6*V);
D = (E/((1+NU)*(1-2*NU)))*[1-NU NU NU 0 0 0 ; NU 1-NU NU 0 0 0 ; NU NU 1-NU 0 0 0 ;
   0 0 0 (1-2*NU)/2 0 0 ; 0 0 0 0 (1-2*NU)/2 0 ; 0 0 0 0 0 (1-2*NU)/2];
k = abs(V)*B'*D*B;
end

等参变换法

刚度矩阵的计算也可以使用更为通用的方法—等参单元的有限元格式，和之前平面单元章节的内容极为类似，格式相同，唯一不同的就是积分点位置和权重。

对于三维四面体单元，刚度矩阵的计算可以表示为：

使用 Hammer 积分点进行数值积分，可以表示为：

其中：为 Hammer 积分点的数量，为积分点的权重系数，是在积分点处的形函数导数矩阵。是雅可比矩阵的行列式值，用于将参考坐标系下的积分转换到物理坐标系。通过对每个 Hammer 积分点进行求和，可以得到单元的刚度矩阵。

四面体单元中，hammer积分点个数与权重系数见下表：

高阶四面体单元 C3D10

单元介绍

C3D10单元是指三维四节点等参四面体单元的二次元素。与C3D4单元相比，C3D10单元有10个节点，包括四个顶点节点和六个边界中点节点，能够更好地捕捉弯曲和变形。单元节点分布见下图：

每个节点的位移同样包含三个分量，分别是 , , ，对应物理坐标系下 , , 的位移。对于任意节点，其位移向量可以表示为矩阵形式：

有限元格式

格式与平面等参单元一致，位移模式：

单元刚度矩阵计算

采用四个 hammer 积分点进行计算

def C3D10Stiffness(prop, elemNodeCoordinate, elemType):
    """
    Calculate the stiffness matrix for a C3D10 element.

    Parameters:
    prop (list): A list containing the material properties (E and nu).
    elemNodeCoordinate (numpy.ndarray): The coordinates of the nodes in the element.
    elemType (str): The type of the element (e.g., 'C3D10').

    Returns:
    k (numpy.ndarray): The stiffness matrix of the element.

    Raises:
    ValueError: If the Jacobian matrix is not square, this exception is raised.
    """
    E = prop[0][0]
    nu = prop[0][1]
    DOF = 3
    elenode = len(elemNodeCoordinate)

    k = np.zeros((DOF * elenode, DOF * elenode))

    hammerWeights, hammerLocations_cols = hammer(4)

    factor = E / ((1 + nu) * (1 - 2 * nu))

    D = factor * np.array([
        [1 - nu,  nu,     nu,     0,          0,          0],
        [nu,     1 - nu, nu,     0,          0,          0],
        [nu,     nu,     1 - nu, 0,          0,          0],
        [0,      0,      0,      (1 - 2 * nu) / 2, 0,    0],
        [0,      0,      0,      0,          (1 - 2 * nu) / 2, 0],
        [0,      0,      0,      0,          0,          (1 - 2 * nu) / 2]
    ])

    for n in range(len(hammerWeights)):
        xi_Hammer = hammerLocations_cols[n][0]
        eta_Hammer = hammerLocations_cols[n][1]
        gama_Hammer = hammerLocations_cols[n][2]

        B, Jacob = Bmatrix(elenode, elemNodeCoordinate,
                           elemType, DOF, xi_Hammer, eta_Hammer, gama_Hammer)

        k += np.dot(np.dot(B.T, D), B) * \
            hammerWeights[n] * np.linalg.det(Jacob)

    return k


def Bmatrix(elenode, elemNodeCoordinate, elemType, DOF, xi_Gauss, eta_Gauss, gama):
    B = np.zeros((6, DOF * elenode))
    _, naturalDerivatives = shapeFunction(xi_Gauss, eta_Gauss, gama, elemType)
    Jacob, XYderivatives = Jacobian(elemNodeCoordinate, naturalDerivatives)
    B[0, np.arange(0, DOF * elenode, 3)] = XYderivatives[:, 0].T
    B[1, np.arange(1, DOF * elenode, 3)] = XYderivatives[:, 1].T
    B[2, np.arange(2, DOF * elenode, 3)] = XYderivatives[:, 2].T
    B[3, np.arange(0, DOF * elenode, 3)] = XYderivatives[:, 1].T
    B[3, np.arange(1, DOF * elenode, 3)] = XYderivatives[:, 0].T
    B[4, np.arange(1, DOF * elenode, 3)] = XYderivatives[:, 2].T
    B[4, np.arange(2, DOF * elenode, 3)] = XYderivatives[:, 1].T
    B[5, np.arange(0, DOF * elenode, 3)] = XYderivatives[:, 2].T
    B[5, np.arange(2, DOF * elenode, 3)] = XYderivatives[:, 0].T
    return B, Jacob