首页/文章/ 详情

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

1月前浏览1083

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

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

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

主要内容有:

  • 常应变四面体单元 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(11);
y1 = elemNodes(12);
z1 = elemNodes(13);
x2 = elemNodes(21);
y2 = elemNodes(22);
z2 = elemNodes(23);
x3 = elemNodes(31);
y3 = elemNodes(32);
z3 = elemNodes(33);
x4 = elemNodes(41);
y4 = elemNodes(42);
z4 = elemNodes(43);
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) / 20,    0],
        [0,      0,      0,      0,          (1 - 2 * nu) / 20],
        [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

支反力的计算

边界条件采用罚函数的方法,支反力可以极为方便的一步求出,具体细节见:
罚函数处理边界条件(理论推导+数值实现)

数值案例

以一个单轴拉伸的模型为例,分别使用 C3D4 和 C3D10 单元进行网格划分并计算,下一期再此基础上,使用 C3D8、C3D20 单元进行网格划分并计算。网格信息以及边界条件设置如下:

计算结果

C3D4 单元

自研程序 UmagAbaqus Umag
自研程序 RFmagAbaqus RFmag

C3D10 单元

自研程序 UmagAbaqus Umag
自研程序 RFmagAbaqus RF mag

精度与 Abaqus 一致。

很高兴认识你

以文字为媒介,分享有限元数值思想的点点滴滴

“积跬步,至千里”



来源:易木木响叮当
ACTAbaqus通用MATLABpythonUM理论
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-09-15
最近编辑:1月前
易木木响叮当
硕士 有限元爱好者
获赞 212粉丝 244文章 344课程 2
点赞
收藏
作者推荐
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈