Abaqus 空间管路 自动化建模与分析

目 录

• Blog links
• 一、前沿
• 二、刚体的空间定位
• • 2.1 欧拉定理
  • 2.2 旋转矩阵
  • 2.3 欧拉角坐标
• 三、几何模型的创建
• • 3.1 空间管路的自由度
  • 3.2 空间管路的旋转矩阵
  • 3.3 局部坐标系下各关键点的坐标
  • 3.4 全局坐标系下各关键点的坐标
  • 3.5 创建管路的几何轴线
  • 3.6 创建管路的几何扫掠面
  • 3.7 创建管路的几何扫掠体
  • 3.9 创建空间管路部件/计算用
• 四、有限元网格的划分
• 五、属性/分析步/加载
• • 5.1 截面属性的指定
  • 5.2 分析步的创建
  • • 5.2.1 模态分析步
    • 5.2.2 静力通用分析步
  • 5.3 定义固定边界条件
  • 5.4 施加内壁均布荷载
  • 5.5 创建作业
• 六、算例
• • 6.1 模态分析
  • 6.2 静力分析
• 七、致谢
• 八、尾声
• 九、参考文献

Blog links

• DalNur | 博客总目录

• ANSYS 参数化建模 实用教程

• Abaqus 二次开发/建模 教程

• nCodeDL 疲劳分析 实用教程

• Python 二次开发 SAP2000 教程

• Python 二次开发 AutoCAD 教程

• Python 二次开发 HyperMesh 教程

• Python 二次开发 Office 教程

• 多体动力学分析 实用教程

• 水动力分析 实用教程

一、前沿

  管路是指液压系统中传输工作流体的管道。相对于管道而言，管路是一种合理安排的管道系统。因为管路的灵活性，管路常被用于液压系统等靠液体驱动的机械设备。

  工程三维建模中，管系建模和电路走线建模都是非常重要的工作，在 AutoCAD、UG、ProE 等软件中，空间管路的创建相对容易，而在 Abaqus 中，因其前处理能力相对较弱，管路的空间定位是不容易的，这就为管路系统的有限元分析带来了困难。

二、刚体的空间定位

  刚体 (Rigid body) 是指在运动中和受力作用后，形状和大小不变，而且内部各点的相对位置不变的物体。绝对刚体实际上是不存在的，只是一种理想模型。

  空间中的刚体具有 6 个自由度，分别为 3 个平动自由度和 3 个转动自由度。3 个平动自由度用于确定刚体的 位置 ，3 个转动自由度用于确定刚体的 姿态 即刚体的朝向。

名词	英文	释义
连体基	body reference frame	固定在刚体上并随其运动的坐标系，用以确定刚体的运动。
参考基	reference frame	全局坐标系

  在三维空间中，确定一个点的位置很容易，只要给出它的三个定位坐标，在笛卡尔坐标系中，为 (x, y, z)。确定一个刚体的位置也很容易，只需要给出刚体上任意一点的坐标，但此时刚体的姿态并没有完全确定，为了确定刚体的姿态，几百年来，科学家们做了持续的努力，也创建了许多方式来描述刚体的姿态，如：旋转向量法、欧拉四元数法、欧拉角法、卡尔丹角法等等。无论采用何种方法，均可描述刚体绕定点做有限转动后的姿态。

2.1 欧拉定理

  我们先讨论刚体转动前后刚体上相应点坐标间的转换关系，刚体平移没必要大篇幅讨论，就涉及加减运算，有啥可说的。

  刚体绕某定点由某位置至另一位置的有限角位移称为刚体的 有限转动 。刚体绕定点的任意有限转动可以由绕过该点某根轴的一次有限转动实现，这就是欧拉定理。即刚体的某一姿态可通过绕空间某一轴 一次转动 有限角度实现。欧拉定理是关于刚体有限转动的重要定理。

2.2 旋转矩阵

  为了方便讨论而又不失一般性，假设刚体转动前的连体基（初态连体基/初态局部坐标系）与参考基（全局坐标系）完全重合。在全局坐标系中，对于绕定轴有限转动的刚体，刚体转动前，其上任意一点 P 的矢径为 r ‾ = ( x 1 , y 1 , z 1 ) T pmb{overline{r}}=(x_1,y_1,z_1)^{rm T} rrr=(x1​,y1​,z1​)T，刚体发生有限转动后，P 点的矢径变为 r = ( x 2 , y 2 , z 2 ) T pmb{r}=(x_2,y_2,z_2)^{rm T} rrr=(x2​,y2​,z2​)T。在此过程中，连体基跟随刚体一同转动。根据刚体运动学基本理论有：

r = A ⋅ r ‾ (1) pmb{r} = pmb{A} cdot pmb{overline{r}} tag{1} rrr=AAA⋅rrr(1)

[ x 2 y 2 z 2 ] = [ a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 ] ⋅ [ x 1 y 1 z 1 ] begin{bmatrix}x_2 \ y_2 \ z_2 \end{bmatrix}= begin{bmatrix}a_{11} & a_{12} & a_{13} \ a_{21} & a_{22} & a_{23} \ a_{31} & a_{32} & a_{33} end{bmatrix} cdot begin{bmatrix}x_1 \ y_1 \ z_1 \ end{bmatrix} ⎣⎡​x2​y2​z2​​⎦⎤​=⎣⎡​a11​a21​a31​​a12​a22​a32​​a13​a23​a33​​⎦⎤​⋅⎣⎡​x1​y1​z1​​⎦⎤​

式中， A pmb{A} AAA 为旋转矩阵 (Rotation Matrix)，实际上它是刚体转动后的连体基（末态连体基/末态局部坐标系）关于参考基的方向余弦阵 (Direction Cosine Matrix) 。即末态连体基在参考基中的姿态由矩阵 A pmb{A} AAA 确定。

  值得注意的是：旋转矩阵中的 9 个元素通常不独立，独立的参数只有3个，所以不是任何矩阵都可以作为旋转矩阵，旋转矩阵需要是正交矩阵 (即逆矩阵等于转置矩阵)。即旋转矩阵必是正交矩阵，但也不是所有正交矩阵都可以作为旋转矩阵。根据描述转动参数的选取不同，旋转矩阵有不同的表达形式，但殊途同归，最终描述的都是同一个刚体定点转动。

Orthogonality of the Rotation Matrix

  作为小白，无需了解式 (1) 的详细推导过程，也没有必要，采用拿来主义，先学会应用，如果特别感兴趣，再去阅读相应专著。这里只需要知道如下结论：在某一参考基（全局坐标系）中，刚体绕定点做有限转动时（刚体姿态发生变化），除转动中心外的（转动中心不一定在刚体上） 刚体上的任意一点 P ，在转动前后，坐标将发生变化，变化前的坐标与变化后的坐标并非毫无任何关系（这都是众所周知的废话），是可以通过末态连体基的方向余弦阵进行换算的，换算公式即为式 (1) 。各种参数并不一定直接关联旋转矩阵。

2.3 欧拉角坐标

  为了简单明了地描述刚体姿态，数百年来，科学大神们创建了各种各样的方法，主要包括：旋转向量法、欧拉四元数法、欧拉角法、卡尔丹角法等。当采用不同的方式描述刚体姿态时，旋转矩阵 A pmb{A} AAA 将有不同的表现形式，它们都可以用来描述同一个旋转。

  当采用欧拉角描述刚体姿态时，用欧拉角表示的旋转矩阵 A = R ( α , β , γ ) pmb{A}=pmb{R}(alpha, beta, gamma) AAA=RRR(α,β,γ) 如下所示：

  关于以欧拉角表示的旋转矩阵，你就先看看形式，知道它大概长啥样就行。接下来，我们将详细说明以欧拉角来描述刚体姿态。为啥详细说它？因为，对于人类来讲，此方法最好理解，方便直观，且应用广泛。

  欧拉角由三个独立参数组成，即 章动角 θ、旋进角(进动角) ψ 和 自转角 φ 。它将欧拉定理中的一次有限转动角度即一次旋转分解为 3 个绕不同轴的旋转。

  任何一个旋转可以表示为依次绕着三个旋转轴旋三个角度的组合，这三个角度称为欧拉角。三个轴可以指固定的世界坐标系轴，也可以指被旋转的物体坐标系的轴。三个旋转轴次序不同，会导致结果不同，即旋转次序影响刚体转动后的姿态，如下图所示。从数学角度解释，则是矩阵乘法不满足交换律。

三、几何模型的创建
3.1 空间管路的自由度

  若想唯一确定管路系统在空间中的走向，除了各关键点处的坐标外，还需要知道管路在弯折处的形态，即非直线段处是如何进行空间过度的。那么可用扭转角、折弯角和弯折半径三个参数来描述管路弯折处的空间形态。

  扭转角：从管路的起始端看过去，上一个平面与下一个平面的顺时针夹角（即上一平面顺时针绕到下一平面需要的角度），平面仅为相邻管线组成的平面，不包括延伸线区域，因此扭转角的范围为 0°~360°。

  折弯角：从管路的起始端看过去，在下一个平面上绕折弯点旋转的角度。由于下一个平面不包括延伸平面区域，因此折弯角无顺逆时针旋转之说，具备唯一性。因此折弯角的范围为 0°~180°。

  根据上述陈述，我们可以制作一折弯表格，用于管路的空间定位。本文以一外径为 10 mm，壁厚为 2 mm，共 8 个直段 7 个折弯的空间管路为例，来介绍其 Abaqus 模型的自动化创建，它的弯折表如下所示。

空间管路折弯表

管段编号	类别	直线段长度 mm	扭转角度 deg	折弯角度 deg	折弯半径 mm
1	直线段	75	—	—	—
2	折弯段	—	0	30	25
3	直线段	75	—	—	—
4	折弯段	—	90	60	25
5	直线段	75	—	—	—
6	折弯段	—	90	50	25
7	直线段	75	—	—	—
8	折弯段	—	90	90	25
9	直线段	75	—	—	—
10	折弯段	—	90	90	25
11	直线段	75	—	—	—
12	折弯段	—	90	90	25
13	直线段	75	—	—	—
14	折弯段	—	90	90	25
15	直线段	75	—	—	—

  基于空间管路的各段的几何信息，我们可在 Python 中定义字典 geoData 用于存贮管路的上述几何信息，geoData 的键名为管段编号，键值为该段管路所对应的几何信息字典，几何信息字典的键名依次为：“LineType” （线型），它的取值可以为 “Straight” 或 “Curve”，分别代表直线段和弯折段；“Length” 为直线段的长度，若为弯折段则此项设置为 “N/A”；“TorsionAngle” 为扭转角，若为直线段则此项设置为 “N/A”；“BendingAngle” 为折弯角，若为直线段则此项设置为 “N/A”；“Radius” 为弯折半径，若为直线段则此项设置为 “N/A”。则上表所对应的 geoData 如下所示。

geoData = {1: {"LineType": "Straight", "Length": 75, "TorsionAngle": "N/A", "BendingAngle": "N/A", "Radius": "N/A"},
           2: {"LineType": "Curve", "Length": "N/A", "TorsionAngle": 0, "BendingAngle": 30, "Radius": 25},
           3: {"LineType": "Straight", "Length": 75, "TorsionAngle": "N/A", "BendingAngle": "N/A", "Radius": "N/A"},
           4: {"LineType": "Curve", "Length": "N/A", "TorsionAngle": 90, "BendingAngle": 60, "Radius": 25},
           5: {"LineType": "Straight", "Length": 75, "TorsionAngle": "N/A", "BendingAngle": "N/A", "Radius": "N/A"},
           6: {"LineType": "Curve", "Length": "N/A", "TorsionAngle": 90, "BendingAngle": 50, "Radius": 25},
           7: {"LineType": "Straight", "Length": 75, "TorsionAngle": "N/A", "BendingAngle": "N/A", "Radius": "N/A"},
           8: {"LineType": "Curve", "Length": "N/A", "TorsionAngle": 90, "BendingAngle": 90, "Radius": 25},
           9: {"LineType": "Straight", "Length": 75, "TorsionAngle": "N/A", "BendingAngle": "N/A", "Radius": "N/A"},
           10: {"LineType": "Curve", "Length": "N/A", "TorsionAngle": 90, "BendingAngle": 90, "Radius": 25},
           11: {"LineType": "Straight", "Length": 75, "TorsionAngle": "N/A", "BendingAngle": "N/A", "Radius": "N/A"},
           12: {"LineType": "Curve", "Length": "N/A", "TorsionAngle": 90, "BendingAngle": 90, "Radius": 25},
           13: {"LineType": "Straight", "Length": 75, "TorsionAngle": "N/A", "BendingAngle": "N/A", "Radius": "N/A"},
           14: {"LineType": "Curve", "Length": "N/A", "TorsionAngle": 90, "BendingAngle": 90, "Radius": 25},
           15: {"LineType": "Straight", "Length": 75, "TorsionAngle": "N/A", "BendingAngle": "N/A", "Radius": "N/A"}, }

def get_information_from_geodata(data):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：根据管道几何信息字典创建空间管路建模所需的输入参数列表/字典信息列表化
    # 参数：data为管路几何信息字典/geoData
    # 返回：输入参数列表/inputdata
    inpdata = []  # 输入数据列表
    radius = []  # 折弯处半径列表
    number = len(data)
    rst1, rst2 = int(number // 2.0), int(number % 2.0)  # 求整数取余数
    for i in range(rst1):
        label1, label2 = int(i * 2.0 + 1), int(i * 2.0 + 2)
        line1type = data[label1]["LineType"]
        line2type = data[label2]["LineType"]
        length1 = data[label1]["Length"]
        # length2 = data[label2]["Length"]
        # tangle1 = data[label1]["TorsionAngle"]
        tangle2 = data[label2]["TorsionAngle"]
        # bangle1 = data[label1]["BendingAngle"]
        bangle2 = data[label2]["BendingAngle"]
        # radius1 = data[label1]["Radius"]
        radius2 = data[label2]["Radius"]
        if line1type == "Straight" and line2type == "Curve":
            inpdata.append([length1, tangle2, bangle2, radius2])
            radius.append(radius2)
    inpdata = np.array(inpdata)
    return inpdata, radius

3.2 空间管路的旋转矩阵

  根据上文的刚体空间定位基本原理，便很容易得到确定管道弯折处姿态的旋转矩阵，获得旋转矩阵的 Python 函数 get_rotation_matrix 如下所示，该函数需要两个参数即扭转角度和折弯角度，函数的返回值为旋转矩阵。

def get_rotation_matrix(tangle, bangle):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：获取旋转矩阵/用于坐标变换
    # 参数：长度、扭转角度、折弯角度、弯折半径
    # 返回：旋转矩阵
    # 说明：局部坐标系以折弯点为原点，下一直管段为局部x轴
    # 注意：角度的单位为弧度
    pi = np.pi  # 圆周率/3.14159265
    # tangle = tangle * pi / 180.0  # 扭转角单位换算/弧度制
    # bangle = bangle * pi / 180.0  # 折弯角度单位换算/弧度制
    axis = [1, 0, 0]  # x轴方向
    currang = tangle  # 当前角度
    currsin = np.sin(currang)  # 正弦值
    currcos = np.cos(currang)  # 余弦值
    x, y, z = axis[0], axis[1], axis[2]
    rotmatrix1 =  . . . （略）
    axis = [0, 0, 1]  # z轴方向
    currang = bangle  # 当前角度
    currsin = np.sin(currang)  # 正弦值
    currcos = np.cos(currang)  # 余弦值
    x, y, z = axis[0], axis[1], axis[2]
    rotmatrix2 = . . . （略）
    rotmatrix = rotmatrix2.dot(rotmatrix1)  # 旋转矩阵
    # invmatrix = np.linalg.inv(rotmatrix)  # 旋转矩阵的逆矩阵
    return rotmatrix

3.3 局部坐标系下各关键点的坐标

def get_local_coords_info_of_points(data):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：获取管路各关键点在局部坐标系下的坐标值
    # 参数：data为输入数据列表
    # 返回：各关键点在局部坐标系下的信息列表
    lcoordsinfo = []
    pi = np.pi  # 圆周率/3.14159265
    for i in range(data.shape[0]):
        length = data[i][0]  # 直线段长度
        tangle = data[i][1] * pi / 180.0  # 扭转角单位换算/弧度制
        bangle = data[i][2] * pi / 180.0  # 折弯角度单位换算/弧度制
        radius = data[i][3]  # 折弯段半径
        rotmatrix = get_rotation_matrix(tangle, bangle)  # 旋转矩阵的逆矩阵
        invmatrix = np.linalg.inv(rotmatrix)  # 旋转矩阵的逆矩阵
        if data.shape[0] == 1:
            blength = length
            location = np.array([blength, 0, 0])
            lcoordsinfo.append([invmatrix, location])
        else:
            if i == 0:
                blength = length + radius * np.tan(bangle / 2.0)
                location = np.array([blength, 0, 0])
                lcoordsinfo.append([invmatrix, location])
            elif i == data.shape[0] - 1:
                bangle0 = data[i - 1][2] * pi / 180.0
                radius0 = data[i - 1][3]
                blength = length + radius0 * math.tan(bangle0 / 2.0)
                location = np.array([blength, 0, 0])
                lcoordsinfo.append([invmatrix, location])
            else:
                bangle0 = data[i - 1][2] * pi / 180.0
                radius0 = data[i - 1][3]
                blength = length + radius * np.tan(bangle / 2.0) + radius0 * math.tan(bangle0 / 2.0)
                location = np.array([blength, 0, 0])
                lcoordsinfo.append([invmatrix, location])
    return lcoordsinfo

3.4 全局坐标系下各关键点的坐标

def get_global_coords_of_keypoints(data):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：获取各关键点在全局坐标系下的坐标值/用于创建wire
    # 参数：data为输入数据列表
    # 返回：关键点在全局坐标系下的信息列表
    # 说明：最后一个折弯点的坐标应排除圆弧段长度
    pointcoords = []
    lcoordsinfo = get_local_coords_info_of_points(data)  # 局部坐标系下各点信息列表
    for i in range(data.shape[0]):
        if i == 0:
            pointcoords.append(lcoordsinfo[i][1])
        else:
            j = i
            point = lcoordsinfo[j][1]
            while j > 0:
                point = np.array(lcoordsinfo[j - 1][0]).dot(point) + lcoordsinfo[j - 1][1]
                j = j - 1
            pointcoords.append(point)
    return pointcoords

3.5 创建管路的几何轴线

def create_axial_wires(mdbname, prtname, points, radius):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：在Abaqus中创建空间管路的轴线部件/创建两个/一个用于扫掠面/另一个用于扫掠体
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；prtname为部件名称
    # 参数：points为管路各点的坐标值；radius为弯折段半径列表
    # 返回：无
    currmdb = mdb.models[mdbname]
    prt1name = prtname + "-1"
    prtobj = currmdb.Part(name=prt1name, dimensionality=THREE_D, type=DEFORMABLE_BODY)
    dplane = prtobj.DatumPointByCoordinate(coords=(0, 0, 0))
    dplids = [dplane.id, ]
    for i in range(len(points)):
        dplane = prtobj.DatumPointByCoordinate(coords=tuple(points[i]))
        dplids.append(dplane.id)
    datums = prtobj.datums
    tmppoints = ()
    for i in range(len(dplids)):
        if i > 0:
            tmppoints = tmppoints + ((datums[dplids[i - 1]], datums[dplids[i]]),)
    prtobj.WirePolyLine(points=tmppoints, mergeType=IMPRINT, meshable=ON)
    vertices = prtobj.vertices
    for i in range(len(points) - 1):
        prtobj.Round(radius=radius[i], vertexList=(vertices.findAt(coordinates=tuple(points[i])),))
    # 创建线及硬点集合/按顺序排列
    edges = prtobj.edges
    setname = "Sweep-Path-Edges"  # 全部边集合/扫略线
    prtobj.Set(edges=edges, name=setname)
    sltes = edges.findAt(((0, 0, 0),))
    number = len(edges)  # 线的条数
    for i in range(number - 1):
        edges_1 = sltes[i].getAdjacentEdges()
        edges_2 = [edge for edge in edges_1 if edge not in sltes]
        sltes = sltes + edges[edges_2[0].index: edges_2[0].index + 1]
    for i, edge in enumerate(sltes):
        point = edge.pointOn
        slteg = edges.findAt(point, )  # 几何边
        setname = "Edge-%s" % (i + 1)
        prtobj.Set(edges=slteg, name=setname)
    vertices = prtobj.vertices
    sltvobjs = vertices.findAt(((0, 0, 0),))
    sltvtids = [sltvobjs[0].index, ]
    for j in range(len(edges)):
        indices = sltes[j].getVertices()
        for index in indices:
            if index not in sltvtids:
                sltvtids.append(index)
    for i, vtid in enumerate(sltvtids):
        vertex = vertices[vtid]
        point = vertex.pointOn
        sltvt = vertices.findAt(point, )  # 几何边
        setname = "Vertex-%s" % (i + 1)
        prtobj.Set(vertices=sltvt, name=setname)  # 全部硬点集合
    prt2name = prtname + "-2"
    currmdb.Part(name=prt2name, objectToCopy=currmdb.parts[prt1name])  # 复制部件
    prt3name = prtname + "-3"
    currmdb.Part(name=prt3name, objectToCopy=currmdb.parts[prt1name])  # 复制部件

空间管路的轴线/Pipe-1



3.6 创建管路的几何扫掠面

def create_sweep_faces(mdbname, prtname, d):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：在Abaqus中创建空间管路的扫掠面/用于切割
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；prtname为部件名称。
    # 参数：d为外径；t为壁厚。
    # 返回：无。
    currmdb = mdb.models[mdbname]
    prtname = prtname + "-2"
    prtobj = currmdb.parts[prtname]  # 扫掠面
    dples = prtobj.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=0.0)
    daxes = prtobj.DatumAxisByPrincipalAxis(principalAxis=ZAXIS)
    datms = prtobj.datums
    pripl = datms[dples.id]
    priax = datms[daxes.id]
    trs = prtobj.MakeSketchTransform(sketchPlane=pripl, sketchUpEdge=priax,
                                     sketchPlaneSide=SIDE1, sketchOrientation=RIGHT, origin=(0.0, 0.0, 0.0))
    x1, x2 = -0.8 * d, 0.8 * d
    skh = currmdb.ConstrainedSketch(name='__profile__', sheetSize=1000.0, gridSpacing=10.0, transform=trs)
    skh.setPrimaryObject(option=SUPERIMPOSE)
    skh.ConstructionLine(point1=(-500.0, 0.0), point2=(500.0, 0.0))  # 构造线
    skh.ConstructionLine(point1=(0.0, -500.0), point2=(0.0, 500.0))  # 构造线
    prtobj.projectReferencesOntoSketch(sketch=skh, filter=COPLANAR_EDGES)
    skh.Line(point1=(x1, 0.0), point2=(x2, 0.0))
    edges = prtobj.sets["Sweep-Path-Edges"].edges
    skh.unsetPrimaryObject()
    pathedges = edges
    prtobj.ShellSweep(path=pathedges, sketchPlane=pripl, sketchUpEdge=priax,
                      sketchOrientation=RIGHT, profile=skh, flipSweepDirection=ON)
    del mdb.models['Model-1'].sketches['__profile__']
    faces = prtobj.faces
    prtobj.Set(faces=faces, name="Faces")

空间管路的切割面/Pipe-2



3.7 创建管路的几何扫掠体

def create_sweep_solids(mdbname, prtname, d, t):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：在Abaqus中创建空间管路的几何体。
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；prtname为部件名称。
    # 参数：d为外径；t为壁厚。
    # 返回：无。
    r1, r2, r3 = 0.5 * d, 0.5 * d - t, 0.5 * (d - t)  # 外半径/内半径/中半径
    prtname = prtname + "-3"
    currmdb = mdb.models[mdbname]
    prtobj = currmdb.parts[prtname]  # 扫掠体
    dples = prtobj.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=0.0)
    daxes = prtobj.DatumAxisByPrincipalAxis(principalAxis=ZAXIS)
    datms = prtobj.datums
    pripl = datms[dples.id]
    priax = datms[daxes.id]
    trs = prtobj.MakeSketchTransform(sketchPlane=pripl, sketchUpEdge=priax,
                                     sketchPlaneSide=SIDE1, sketchOrientation=RIGHT, origin=(0.0, 0.0, 0.0))
    skh = currmdb.ConstrainedSketch(name='__profile__', sheetSize=1000.0, gridSpacing=10.0, transform=trs)
    skh.setPrimaryObject(option=SUPERIMPOSE)
    skh.ConstructionLine(point1=(-500.0, 0.0), point2=(500.0, 0.0))  # 构造线
    skh.ConstructionLine(point1=(0.0, -500.0), point2=(0.0, 500.0))  # 构造线
    prtobj.projectReferencesOntoSketch(sketch=skh, filter=COPLANAR_EDGES)
    skh.CircleByCenterPerimeter(center=(0.0, 0.0), point1=(r2, 0.0))
    skh.CircleByCenterPerimeter(center=(0.0, 0.0), point1=(r1, 0.0))
    skh.unsetPrimaryObject()
    edges = prtobj.sets["Sweep-Path-Edges"].edges
    skh.unsetPrimaryObject()
    pathedges = edges
    setnames = prtobj.sets.keys()
    if len(setnames) == 4:
        prtobj.SolidSweep(path=pathedges, sketchPlane=pripl, sketchUpEdge=priax,
                          sketchOrientation=RIGHT, profile=skh, flipSweepDirection=OFF)
    else:
        prtobj.SolidSweep(path=pathedges, sketchPlane=pripl, sketchUpEdge=priax,
                          sketchOrientation=RIGHT, profile=skh, flipSweepDirection=ON)
    del currmdb.sketches['__profile__']
    cells = prtobj.cells
    prtobj.Set(cells=cells, name="Cells")
    # 切分扫掠体
    radius = r1 + 0.01
    currmdb = mdb.models[mdbname]
    prtobj = currmdb.parts[prtname]  # 扫掠面
    setnames = prtobj.sets.keys()
    vs, fs, es, cs = prtobj.vertices, prtobj.faces, prtobj.edges, prtobj.cells
    number = int((len(setnames) - 3) / 2) + 1  # 硬点总数
    for i in range(number):
        vsetname = "Vertex-%s" % (i + 1)
        point = prtobj.sets[vsetname].vertices[0].pointOn[0]
        edges = es.getByBoundingBox(point[0] - radius, point[1] - radius,
                                    point[2] - radius, point[0] + radius,
                                    point[1] + radius, point[2] + radius)
        for edge in edges:
            if abs(edge.getRadius() - r1) < 0.001:
                i = i + 1
                setname = "Partition-Edge-%s" % (i)
                prtobj.Set(edges=es[edge.index:edge.index + 1], name=setname)
                try:
                    findex = edge.getFaces()
                    cindex = prtobj.faces[findex[0]].getCells()
                    pickededges = (edge,)
                    prtobj.PartitionCellByPatchNEdges(cell=cs[cindex[0]], edges=pickededges)
                except:
                    pass
    cells = prtobj.cells
    for i, cell in enumerate(cells):
        point = cell.pointOn[0]
        sltcl = cells.findAt((point,))
        setname = "Cell-%s" % (number - i - 1)
        prtobj.Set(cells=sltcl, name=setname)
    setnames = []
    for i in range(number):
        vsetname = "Vertex-%s" % (i + 1)
        point = prtobj.sets[vsetname].vertices[0].pointOn[0]
        edges = es.getByBoundingBox(point[0] - radius, point[1] - radius,
                                    point[2] - radius, point[0] + radius,
                                    point[1] + radius, point[2] + radius)
        setname = "Circle-Edges-%s" % (i + 1)
        prtobj.Set(edges=edges, name=setname)
        setnames.append(setname)
    currmdb = mdb.models[mdbname]
    prtobj = currmdb.parts[prtname]
    setobjs = [prtobj.sets[setname] for setname in setnames]
    prtobj.SetByBoolean(name="Circle-Edges", sets=setobjs)

空间管路的扫掠体/Pipe-3



3.9 创建空间管路部件/计算用

def create_pipe_part(mdbname, prtname, d, t):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：创建最终计算用的管道部件
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；prtname为部件名称。
    # 参数：d为外径；t为壁厚。
    # 返回：无。
    r1, r2, r3 = 0.5 * d, 0.5 * d - t, 0.5 * (d - t)  # 外半径/内半径/中半径
    prt1name, prt2name, prt3name = prtname + "-1", prtname + "-2", prtname + "-3"
    currmdb = mdb.models[mdbname]
    prt2obj = currmdb.parts[prt2name]
    prt3obj = currmdb.parts[prt3name]
    asm = currmdb.rootAssembly
    asm.Instance(name=prt2name, part=prt2obj, dependent=ON)
    asm.Instance(name=prt3name, part=prt3obj, dependent=ON)
    prt4name = prtname + "-4"
    asm.InstanceFromBooleanMerge(name=prt4name, instances=(asm.instances[prt2name], asm.instances[prt3name]),
                                 keepIntersections=ON, originalInstances=SUPPRESS, domain=GEOMETRY)
    currmdb = mdb.models[mdbname]
    prt4obj = currmdb.parts[prt4name]
    faces = prt4obj.sets["Faces"].faces
    fs1 = faces.findAt(((0.0, 1.6 * r1, 0.0),))
    sltfs = fs1[0].getFacesByFaceAngle(angle=30.0)
    prt4obj.Set(faces=sltfs, name='temp-1-Faces')
    fs2 = faces.findAt(((0.0, 0.0, 0.0),))
    sltfs = fs2[0].getFacesByFaceAngle(angle=30.0)
    prt4obj.Set(faces=sltfs, name='temp-2-Faces')
    fs3 = faces.findAt(((0.0, -1.6 * r1, 0.0,),))
    sltfs = fs3[0].getFacesByFaceAngle(angle=30.0)
    prt4obj.Set(faces=sltfs, name='temp-3-Faces')
    dfaces = prt4obj.sets["temp-1-Faces"].faces
    prt4obj.RemoveFaces(faceList=dfaces, deleteCells=False)
    dfaces = prt4obj.sets["temp-2-Faces"].faces
    prt4obj.RemoveFaces(faceList=dfaces, deleteCells=False)
    dfaces = prt4obj.sets["temp-3-Faces"].faces
    prt4obj.RemoveFaces(faceList=dfaces, deleteCells=False)
    # # 创建最终计算
    currmdb = mdb.models[mdbname]
    asm = currmdb.rootAssembly
    prt1obj = currmdb.parts[prt1name]
    asm.Instance(name=prt1name, part=prt1obj, dependent=ON)
    istname = prt4name + "-1"
    asm.InstanceFromBooleanMerge(name="DalNur-ZDW", instances=(asm.instances[prt1name], asm.instances[istname]),
                                 keepIntersections=ON, originalInstances=SUPPRESS, domain=GEOMETRY)
    currmdb = mdb.models[mdbname]
    asm = currmdb.rootAssembly
    istnames = asm.instances.keys()
    asm.deleteFeatures(istnames)
    currmdb.Part(name=prtname, objectToCopy=currmdb.parts["DalNur-ZDW"])  # 复制部件
    del currmdb.parts['DalNur-ZDW']
    # 创建集合/布种子边集合
    currmdb = mdb.models[mdbname]
    prtobj = currmdb.parts[prtname]
    edges = prtobj.edges
    slt1es, slt2es = edges[:1], edges[:1]
    for edge in edges:
        length = edge.getSize()
        if length > t - 0.01 and length < t + 0.01:
            point = edge.pointOn[0]
            tempe = edges.findAt((point,))
            slt1es = slt1es + tempe
        try:
            radius = edge.getRadius()
            print(radius)
        except:
            point = edge.pointOn[0]
            tempe = edges.findAt((point,))
            slt2es = slt2es + tempe
    slt1es = slt1es[1:]
    setname = "Radius-Edges"
    prtobj.Set(edges=slt1es, name=setname)
    slt2es = slt2es[1:]
    setname = "Straight-Edges"
    prtobj.Set(edges=slt2es, name=setname)
    prtobj.SetByBoolean(name="Length-Edges", operation=DIFFERENCE,
                        sets=(prtobj.sets["Straight-Edges"], prtobj.sets["Radius-Edges"],))
    # 创建集合/加载/边界条件
    prtobj = currmdb.parts[prtname]
    setnames = prtobj.sets.keys()
    vs, fs = prtobj.vertices, prtobj.faces
    number = int((len(setnames) - 10) / 5) + 1  # 硬点总数
    vobj1, vobj2 = prtobj.sets["Vertex-1"].vertices[0], prtobj.sets["Vertex-2"].vertices[0]
    center1 = vobj1.pointOn[0]
    center2 = vobj2.pointOn[0]
    faces = fs.getByBoundingCylinder(center1=(center1[0] - 1, center1[1], center1[2]),
                                     center2=(center2[0] + 1, center2[1], center2[2]), radius=r3)
    side1Faces = faces[0].getFacesByFaceAngle(angle=30.0)
    prtobj.Surface(side1Faces=side1Faces, name='Surf-Load')
    face1 = fs.getByBoundingSphere(center=(0, 0, 0), radius=r1 + 0.001)
    setname = "Vertex-%s" % number
    vobj = prtobj.sets[setname].vertices[0]
    center = vobj.pointOn[0]
    x, y, z = center
    face2 = fs.getByBoundingSphere(center=(x, y, z), radius=r1 + 0.001)
    prtobj.Set(faces=face1, name='Face-BC-1')
    prtobj.Set(faces=face2, name='Face-BC-2')
    prtobj.Surface(side1Faces=face1, name='Surf-BC-1')
    prtobj.Surface(side1Faces=face2, name='Surf-BC-2')
    prtobj.Set(faces=side1Faces, name='Face-Load')
    currmdb = mdb.models[mdbname]
    prtobj = currmdb.parts[prtname]
    vertices = prtobj.vertices
    xMin, xMax = -0.001, 0.001
    yMin, yMax = -10241024.0, 10241024.0
    zMin, zMax = -10241024.0, 10241024.0
    sltvts = vertices.getByBoundingBox(xMin=xMin, xMax=xMax,
                                       yMin=yMin, yMax=yMax,
                                       zMin=zMin, zMax=zMax)
    prtobj.Set(vertices=sltvts, name='temp-vertices')
    currmdb = mdb.models[mdbname]
    prtobj = currmdb.parts[prtname]
    edge = prtobj.sets["Sweep-Path-Edges"].edges
    prtobj.RemoveWireEdges(wireEdgeList=edge)

四、有限元网格的划分

def create_mesh(mdbname, prtname, lsize=0, csize=None, cnumber=None, rnumber=None, bnumber=None):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：在Abaqus中创建空间管路的体网格
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；prtname为部件名称。
    # 参数：lsize,为直线段沿长度方向的网格尺寸
    # 参数：csize, rsize分别为长度方向、圆周方向、半径方向的网格尺寸
    # 参数：cnumber为圆周方向的网格尺寸
    # 参数：bnumber为弯折段的网格数目/各段统一数目
    # 返回：无。
    currmdb = mdb.models[mdbname]
    prtobj = currmdb.parts[prtname]
    setobj = prtobj.sets["temp-vertices"]
    radius, i = 0, 0
    for vt in setobj.vertices:
        x, y, z = vt.pointOn[0]
        tmpr = (x ** 2 + y ** 2 + z ** 2) ** 0.5
        if tmpr >= radius:
            radius = tmpr
    cumber1, cnumber2 = 5, 5  # 圆周方向的网格不等能小于5个
    if cnumber is not None:
        cumber1 = cnumber
    if csize is not None:
        if float(csize) > 0:
            cnumber2 = int(2 * 3.14159265 * radius / csize + 1)  # 圆周方向上的网格数目
    cnumber = int(max(cumber1, cnumber2))  # 最终建模用网格
    if bnumber < 5:
        bnumber = 5
    currmdb = mdb.models[mdbname]
    prtobj = currmdb.parts[prtname]
    edges, cells = prtobj.edges, prtobj.cells
    prtobj.setMeshControls(regions=cells, technique=STRUCTURED)  # 结构网格
    prtobj.seedEdgeByNumber(edges=edges, number=bnumber, constraint=FINER) # 弯折段网格尺寸
    edges = prtobj.sets["Circle-Edges"].edges
    prtobj.seedEdgeByNumber(edges=edges, number=cnumber, constraint=FINER)
    edges = prtobj.sets["Straight-Edges"].edges
    prtobj.seedEdgeBySize(edges=edges, size=lsize, constraint=FINER)
    edges = prtobj.sets["Radius-Edges"].edges
    prtobj.seedEdgeByNumber(edges=edges, number=rnumber, constraint=FINER)
    prtobj.generateMesh()
    # 删除多余集合
    currmdb = mdb.models[mdbname]
    prtobj = currmdb.parts[prtname]
    setnames = prtobj.sets.keys()
    setname = "Sweep-Path-Edges"
    del currmdb.parts[prtname].sets[setname]
    setname = "temp-vertices"
    del currmdb.parts[prtname].sets[setname]
    for setname in setnames:
        if setname[0] == "E" and "Edge-" in setname:
            del currmdb.parts[prtname].sets[setname]
        if "Partition-Edge-" in setname:
            del currmdb.parts[prtname].sets[setname]
        if "Vertex-" in setname:
            del currmdb.parts[prtname].sets[setname]
    # 装配部件
    prtobj = currmdb.parts[prtname]
    asm = currmdb.rootAssembly
    asm.Instance(name=prtname, part=prtobj, dependent=ON)

五、属性/分析步/加载
5.1 截面属性的指定

def assign_section_props(mdbname="Model-1", prtname="Pipe", matname="Mat-ZDW", rho=1024.1024, e=1024.1024,
                         nu=0.314159265):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：为管道部件指定截面属性
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；prtname为部件名称
    # 参数：matname为材料名称；rho、e、nu分别为质量密度、弹性模量和泊松比。
    # 返回：无。
    currmdb = mdb.models[mdbname]
    currmdb.Material(name=matname)
    rho, e, nu = float(rho), float(e), float(nu)
    currmdb.materials[matname].Density(table=((rho,),))
    currmdb.materials[matname].Elastic(table=((e, nu),))
    secname = "Section-Pipe"
    currmdb.HomogeneousSolidSection(name=secname, material=matname, thickness=None)
    try:
        currmdb.materials["Mat-ZDW"].Density(table=((rho,),))
        currmdb.materials["Mat-ZDW"].Elastic(table=((e, nu),))
        secname = "Section-DalNur"
        currmdb.HomogeneousSolidSection(name=secname, material=matname, thickness=None)
    except:
        pass
    prtobj = currmdb.parts[prtname]
    region = prtobj.sets['Cells']
    prtobj.SectionAssignment(region=region, sectionName=secname, offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='',
                             thicknessAssignment=FROM_SECTION)

5.2 分析步的创建
5.2.1 模态分析步

def create_freq_step(mdbname, freqnum=10):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：创建模态分析步
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；
    # 参数：stepname为模态分析步的名称；freqnum为模态阶数
    # 返回：无。
    currmdb = mdb.models[mdbname]
    stepname = "Step-1"
    currmdb.FrequencyStep(name=stepname, previous='Initial', numEigen=freqnum)

5.2.2 静力通用分析步

def create_static_step(mdbname, incresize=0.1, nlgem=True):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：创建通用静力分析步
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；
    # 参数：incresize为增量步步长，需要小于1.0。
    # 参数：nlgem为几何非线性选项，若为True，则开启几何非线性；若为False则关闭几何非线性。
    # 返回：无。
    currmdb = mdb.models[mdbname]
    if incresize >= 1.0:
        incresize = 0.5
    if nlgem:
        currmdb.StaticStep(name='Step-2', previous='Step-1', maxNumInc=10000, timeIncrementationMethod=FIXED,
                           initialInc=incresize, noStop=OFF, nlgeom=ON)
    else:
        currmdb.StaticStep(name='Step-2', previous='Step-1', maxNumInc=10000, timeIncrementationMethod=FIXED,
                           initialInc=incresize, noStop=OFF)

5.3 定义固定边界条件

def create_boundary_condition(mdbname, prtname):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：定义左右两端固定边界条件
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；prtname为部件名称
    # 返回：无。
    currmdb = mdb.models[mdbname]
    asm = currmdb.rootAssembly
    region = asm.instances[prtname].sets['Face-BC-1']
    currmdb.DisplacementBC(name='BC-1', createStepName='Initial',
                           region=region, u1=SET, u2=SET, u3=SET, ur1=SET, ur2=SET, ur3=SET,
                           amplitude=UNSET, distributionType=UNIFORM, fieldName='', localCsys=None)
    region = asm.instances[prtname].sets['Face-BC-2']
    currmdb.DisplacementBC(name='BC-2', createStepName='Initial',
                           region=region, u1=SET, u2=SET, u3=SET, ur1=SET, ur2=SET, ur3=SET,
                           amplitude=UNSET, distributionType=UNIFORM, fieldName='', localCsys=None)

5.4 施加内壁均布荷载

def apply_pressure_load(mdbname, prtname, pressure=21929.0):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：时间管道内壁的均布荷载
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；prtname为部件名称
    # 参数：pressure为均布荷载的取值
    # 返回：无。
    currmdb = mdb.models[mdbname]
    asm = currmdb.rootAssembly
    region = asm.instances[prtname].surfaces['Surf-Load']
    pressure = float(pressure)
    currmdb.Pressure(name='Load-1', createStepName='Step-2', region=region, distributionType=UNIFORM, field='',
                     magnitude=pressure, amplitude=UNSET)

5.5 创建作业

def create_job(mdbname, jobname=None, submit=False):
    # 作者：Zhao Dongwei、DalNur
    # 邮箱：18700910380@163.com
    # 功能：创建分析作业
    # 参数：mdbname模型数据库的名称，默认为"Model-1"；job为作业名称。
    # 参数：submit为作业提交选项，若为True，则提交作业，自动计算。
    # 返回：无
    if jobname is None:
        jobname = "Pipe-20210929-133523"
        mdb.Job(name=jobname, model=mdbname)
    try:
        mdb.Job(name=jobname, model=mdbname)
    except:
        pass
    if submit:
        try:
            mdb.jobs[jobname].submit(consistencyChecking=OFF)
        except:
            pass

六、算例
6.1 模态分析

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
    ============================
    Author: Zhao Dongwei、DalNur
    Email: 18700910380@163.com
    Email: liyang@alu.hit.edu.cn
    Date: 2021/09/29 14:18:22
    ============================
"""

from abaqus import *
from abaqusConstants import *
from caeModules import *
import numpy as np


data = geoData
[inpdata, radius] = get_information_from_geodata(data)
mdbname = "Model-1"  # 模型数据库名称
prtname = "Pipe"  # 管道部件名称
d, t = 12, 3  # 外径壁厚
points = get_global_coords_of_keypoints(data=inpdata)
create_axial_wires(mdbname, prtname, points, radius)
create_sweep_faces(mdbname, prtname, d)
create_sweep_solids(mdbname, prtname, d, t)
create_pipe_part(mdbname, prtname, d, t)
lsize = 10  # 直线段的网格近似尺寸
csize = None  # 圆周方向上的网格近似尺寸
cnumber = 30  # 圆周方向上的网格近似数目
rnumber = 5  # 半径方向上的网格近似数目
bnumber = 10  # 弯折处网格尺寸/需反复调整
create_mesh(mdbname, prtname, lsize, csize, cnumber, rnumber, bnumber)
matname = "Mat-Pipe"
rho, e, nu = 7850e-9, 2.06e5, 0.3
assign_section_props(mdbname, prtname, matname, rho, e, nu)
create_boundary_condition(mdbname, prtname)
create_freq_step(mdbname, freqnum=10)
# create_static_step(mdbname, incresize=0.1, nlgem=False)
# apply_pressure_load(mdbname, prtname, pressure=2021.0929)
jobname = "Test-Pipe-Job-1"
create_job(mdbname, jobname, submit=True)

第8阶模态



6.2 静力分析

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
    ============================
    Author: Zhao Dongwei、DalNur
    Email: 18700910380@163.com
    Email: liyang@alu.hit.edu.cn
    Date: 2021/09/29 14:18:22
    ============================
"""

from abaqus import *
from abaqusConstants import *
from caeModules import *
import numpy as np


data = geoData
[inpdata, radius] = get_information_from_geodata(data)
mdbname = "Model-1"  # 模型数据库名称
prtname = "Pipe"  # 管道部件名称
d, t = 12, 3  # 外径壁厚
points = get_global_coords_of_keypoints(data=inpdata)
create_axial_wires(mdbname, prtname, points, radius)
create_sweep_faces(mdbname, prtname, d)
create_sweep_solids(mdbname, prtname, d, t)
create_pipe_part(mdbname, prtname, d, t)
lsize = 10  # 直线段的网格近似尺寸
csize = None  # 圆周方向上的网格近似尺寸
cnumber = 30  # 圆周方向上的网格近似数目
rnumber = 5  # 半径方向上的网格近似数目
bnumber = 10  # 弯折处网格尺寸/需反复调整
create_mesh(mdbname, prtname, lsize, csize, cnumber, rnumber, bnumber)
matname = "Mat-Pipe"
rho, e, nu = 7850e-9, 2.06e5, 0.3
assign_section_props(mdbname, prtname, matname, rho, e, nu)
create_boundary_condition(mdbname, prtname)
# create_freq_step(mdbname, freqnum=10)
create_static_step(mdbname, incresize=0.1, nlgem=False)
apply_pressure_load(mdbname, prtname, pressure=2021.0929)
jobname = "Test-Pipe-Job-2"
create_job(mdbname, jobname, submit=True)

Mises应力云图/全局

Mises应力云图/管壁内表面



七、致谢

  衷心地感谢 CALT 的 Fang Hongrong、Wu Yuanhao 和 Wang Ruwen！

  让我有机会能接触并了解到这个领域的相关工作。

  同时，没能加入你们实属遗憾，祝你们工作顺利，天天开心！

  本文最重要的建模架构来自 CALT 的 Zhao Dongwei，祝你前程似锦，万里鹏程！

  代码的总成及部分修改工作由阿阳完成。

八、尾声

  以上便是自动化创建空间管路 Abaqus 模型的基本内容。

  自 2020-12-10 00:17:00 至 2021-09-29 13:35:23，

  经过努力，已基本实现空间管路在 Abaqus 中的自动建模；

  特此记录，借以总结，同时也能方便后学者。

  本文仅用于个人学习，除此之外，无其他任何用途。

  因个人水平有限，文中难免有所疏漏，还请各位大神不吝批评指正。

  如有疑问，请联系邮件联系，Email： 18700910380@163.com 。

  胸藏文墨怀若谷，腹有诗书气自华，希望各位都能在知识的 pāo 子里快乐徜徉。

  本文首次发表于 2021-09-29 14:42:09，Beijing 。

  欢迎大家点赞、评论及转载，转载请注明出处！

  为我打call，不如为我打款！

  最后，祝各位攻城狮们，珍爱生命，保护发际线！

九、参考文献

[1]. 刚体. 百度百科.

[2]. 计算多体系统动力学. 洪嘉振.

[3]. 刚体在三维空间的旋转(关于旋转矩阵、DCM、旋转向量、四元数、欧拉角). MulinB.

[4]. 旋转矩阵、欧拉角、四元数理论及其转换关系. AI人工智能科学.

[5]. 各种旋转方式总结. Ladd7.

[6]. 旋转变换（二）欧拉角. csxiaoshui.

.