CN107944135B - 三维长管的支撑位置优化方法 - Google Patents

Abstract

三维长管的支撑位置优化方法，涉及航空航天以及工业生产领域，特别涉及管道支撑位置的优化方法。为了解决现有的人为选取或等间距选取支撑点支撑长管的方法存在长管容易变形的问题。本发明利用相邻的节点表示一个直管单元，计算在非支撑点的节点和支撑点的节点处的总体坐标下对直管单元负载向量，根据局部坐标系下直管单元的刚度矩阵计算总体坐标系下的刚度矩阵，并根据总体坐标系下刚度矩阵计算三维长管由n个直管单元组成时弯管在总坐标系下的总刚度矩阵；根据直管单元负载向量计算总体坐标系下总载荷向量，根据式K_tD_total＝F_total计算得到总变形矢量；最后通过PSO算法对支撑位置进行优化得到最优的支撑位置。本发明适用于长管的支撑位置选择。

Description

三维长管的支撑位置优化方法

技术领域

本发明涉及航空航天以及工业生产领域，特别涉及管道支撑位置的优化方法。

背景技术

工业及航空航天领域经常会使用到管道，管道可能承受许多种外力的作用，包括本身的重量等。为了保证管道的强度和刚度，必须设置各种支(吊)架、支座等支撑装置。支架的设置根据管道的直径、材质、管子壁厚和载荷等条件决定。尤其是针对航空航天用的管道，由于航空航天的设备要求，虽然其本身一般会使用高强度金属，但是所用到的管道的管道壁较薄，在管道较长的情况下，薄壁长管存在由于自重而产生薄壁长管变形过大的问题，尤其是薄壁长管以后还有可能会在一些特殊条件下工作，例如自身重量等条件以及工作压力等工况条件，极有可能造成薄壁长管产生变形。但是目前为止，针对管道的支撑都是人为选取或等间距选取支撑点进行支撑。很容易造成薄壁长管支撑中产生较大变形，尤其是针对航空航天领域使用的薄壁长管，变形问题更加突出。

发明内容

本发明为了解决现有的人为选取或等间距选取支撑点支撑长管的方法存在长管容易变形的问题。

三维长管的支撑位置优化方法，包括下述步骤：

步骤1、将三维长管模型输入二次开发工具，提取三维长管模型的轴线；

步骤2、基于有限元法的三维长管刚度矩阵的建立：

在二次开发工具中三维长管的总体坐标系为XYZ，将三维长管分解成若干个直管单元，直管单元对应的局部坐标系为xyz；

局部坐标系下直管单元的刚度矩阵

h₁至h₁₀为中间变量，

h₁＝EA/l，h₄＝GJ/l，

式中，E为材料的杨氏模量，A为弯管截面面积，l为单元杆长度，b为刚度修正系数，G为材料的剪切模量，I_y表示截面对y轴的惯性矩、I_z表示截面对z轴的惯性矩，J表示截面的极惯性矩；

利用相邻的节点表示一个直管单元，在总体坐标下的直管单元负载向量

F⁽ⁱ⁾＝(F₁ ⁽ⁱ⁾F₂ ⁽ⁱ⁾)

在非支撑点的节点处，直管单元只受重力产生均布载荷q_i，对应直管单元长度为l_i，非支撑点处的载荷如下

在支撑点的节点处载荷为(000000)^T；

根据局部坐标系下直管单元的刚度矩阵k计算总体坐标系下的刚度矩阵：

K＝RkR^T

其中：R^T表示R的转置矩阵； cos＜·,·＞表示局部坐标系与总体坐标系的夹角；

当三维长管由n个直管单元组成时，弯管在总坐标系XYZ下的总刚度矩阵K_t：

其中上角标的i代表直管单元i；表示直管单元1对应的K的分块矩阵形式，表示直管单元i对应的K的分块矩阵形式；

根据非支撑点处和支撑点的节点处直管单元负载向量，计算总体坐标系下总载荷向量F_total：

根据式K_tD_total＝F_total计算得到总变形矢量D_total，总变形矢量D_total即为弯管的最大变形量；

步骤3、以每种弯管支撑位置组合下的弯管最大总变形矢量为适应度函数，通过PSO算法对支撑位置进行优化寻找适应度函数即最大总变形矢量最小的支撑位置，得到最优的支撑位置。

进一步地，步骤3的具体过程如下：

将弯管的每一种支撑位姿看作一个粒子，包含弯管的支撑位置和弯管的旋转角度；适应度是指每种弯管支撑位置组合下的弯管最大变形量；

首先初始化一群随机粒子，所有粒子组成一个群体，再通过迭代更新每个粒子的速度和位置；粒子个体自身找到的最优解，称为个体最优pbest；整个群体到目前为止找到的最优解，称为群体最优gbest；在每一次更新中，粒子通过跟踪个体最优pbest和群体最优gbest来更新自己的状态；

对于第k次迭代，PSO中的第i个粒子按照以下公式进行更新：

其中，i＝1,2,3,...,M；M为粒子群中粒子的总数；是粒子i第k-1次迭代的飞行速度矢量；是粒子i第k-1次迭代后的位置；是第k-1次迭代后，粒子i历史上最优位置；gbest^k是第k-1次迭代后，所有粒子中的最优位置；c1，c2为权重因子；rand()为0到1内的随机数，w为惯性权重。

进一步地，步骤3的粒子X＝(n₁,n₂,...,n_i,a)，n_i表示第i个支撑选取的位置，a表示弯管的旋转角度。

进一步地，步骤2中的总体坐标系下的刚度矩阵K＝RkR^T是基于三维长管的总体坐标系和直管单元对应的局部坐标系的得到的，具体过程如下：

在二次开发工具中，在总体坐标系下的总体刚度矩阵K、负载向量F和变形矢量D，关系如下：

KD＝F (1)

在局部坐标系下直管单元的刚度矩阵k、负载向量f和变形矢量d，关系如下：

kd＝f (2)

利用相邻的节点表示一个直管单元；节点直管单元1包括节点1、节点2，直管单元2包括节点2、节点3；依次类推，直管单元i包括节点i、节点(i+1)；

局部坐标系下变形矢量d为：

d＝(u₁ V₁ W₁ a_x1 v_y1 v_z1 u₂ V₂ W₂ a_x2 v_y2 v_z2)^T (4)

其中，参数下标中的“1”和“2”表示相邻节点中的节点1和节点2；u₁为节点1的x方向的位移，V₁为节点1的y方向的位移，W₁为节点1的z方向的位移，a_x1为节点1的x方向的转角，v_y1为节点1的y方向的转角，v_z1为节点1的z方向的转角；u₂为节点2的x方向的位移，V₂为节点2的y方向的位移，W₂为节点2的z方向的位移，a_x2为节点2的x方向的转角，v_y2为节点2的y方向的转角，v_z2为节点2的z方向的转角；

局部坐标系下负载向量f为

f＝(f_x1 F_y1 F_z1 M_x1 M_y1 M_z1 f_x2 F_y2 F_z2 M_x2 M_y2 M_z2)^T (5)

其中，f_x1、F_y1、F_z1、M_x1、M_y1、M_z1分别表示节点1在xyz三个方向的力与力矩；f_x2、F_y2、F_z2、M_x2、M_y2、M_z2分别表示节点2在xyz三个方向的力与力矩；

在总体坐标下的直管单元负载向量F为

F＝(f_x1' F_y1' F_z1' M_x1' M_y1' M_z1' f_x2' F_y2' F_z2' M_x2' M_y2' M_z2')^T (6)

其中，f_x1'、F_y1'、F_z1'、M_x1'、M_y1'、M_z1'分别表示节点1在XYZ三个方向的力与力矩；f_x2'、F_y2'、F_z2'、M_x2'、M_y2'、M_z2'分别表示节点2在XYZ三个方向的力与力矩；

根据总体坐标系XYZ与局部坐标系xyz间的空间位置关系可得到总体坐标系下变形矢量D和负载向量F与直管单元局部坐标系下变形矢量d和负载向量f的关系如式：

d＝R^TD (7)

f＝R^TF (8)

根据公式(1)、(2)、(7)、(8)得到总体坐标系下的刚度矩阵K＝RkR^T。

本发明具有以下有益效果：

本发明提供了一种快速准确的找到最优支撑位置的方法，使管件在不同支撑状态下由自重产生的变形量最小，从而提高管件的形状精度。尤其是针对航空航天用的薄壁长管，能够有效的防止长管变形，极大的提高管件的形状精度，从而保证安装精度。

基于实施例的仿真条件下，利用本发明计算出的管件最大变形量与ANSYS分析得到的最大变形量的平均偏差为0.07mm，最大偏差为0.13mm。在整个方法单纯从运算速度上相较传统的有限元分析软件也大大提高，以本发明的实施例来说，从开始运行到得到最优运算结果及生成迭代曲线等后处理结果，平均只需要30秒以内，这包括计算了500种支撑情况并从中自动选出最优的支撑位置，而传统有限元软件计算一次支撑情况就需要近1分钟。

附图说明

图1为三维长管模型在ProE中的示意图；

图2为ProE对轴线进行提取的示意图；

图3为任意直管单元在总体坐标系下的受力情况示意图；

图4为管件节点位置及旋转角度示意图；

图5(a)至图5(f)为ANSYS管件变形量分析结果；

图6为PSO迭代曲线；

图7为PSO粒子群算法得到的最优支撑位置；

图8(a)至图8(c)为最优支撑位置在x-y-z方向的分量沿弯管轴线变化曲线。

具体实施方式

具体实施方式一：

本发明所称的三维长管可以是空间圆管，可以是其他形状的长管，如截面是方形的长管。本发明可以适用于各种壁厚的长管，由于厚壁的长管不易产生变形，而薄壁长管很容易产生变形问题，所以本发明尤其适用于三维薄壁长管。

三维长管的支撑位置优化方法，包括下述步骤：

步骤1、将三维长管模型输入二次开发工具(ProE等)，提取三维长管模型的轴线；

二次开发工具(ProE)中三维长管如图1所示，二次开发技术对轴线进行提取如图2所示；

步骤2、基于有限元法的三维长管刚度矩阵的建立：

在二次开发工具中三维长管的总体坐标系为XYZ，将三维长管分解成若干个直管单元，直管单元对应的局部坐标系为xyz；

根据材料力学，可以得到局部坐标系下直管单元的刚度矩阵

h₁至h₁₀为中间变量，

h₁＝EA/l，h₄＝GJ/l，

式中，E为材料的杨氏模量，A为弯管截面面积，l为单元杆长度，b为刚度修正系数，G为材料的剪切模量，I_y表示截面对y轴的惯性矩、I_z表示截面对z轴的惯性矩，J表示截面的极惯性矩；

利用相邻的节点表示一个直管单元，在总体坐标下的直管单元负载向量

F⁽ⁱ⁾＝(F₁ ⁽ⁱ⁾F₂ ⁽ⁱ⁾)

在非支撑点的节点处，直管单元只受重力产生均布载荷q_i，对应直管单元长度为l_i，非支撑点处的载荷如下

在支撑点的节点处载荷为(000000)^T；

非支撑点处和支撑点处对应的载荷与均布载荷并不是上下位概念，是材料力学上的术语，非支撑点处和支撑点处对应的载荷包括三个方向上的力和三个方向上的转矩；

根据局部坐标系下直管单元的刚度矩阵k计算总体坐标系下的刚度矩阵：

K＝RkR^T

其中：R^T表示R的转置矩阵(实际应该是R^-1，这里R^T和R^-1是一样的)；cos＜·,·＞表示局部坐标系与总体坐标系的夹角，cos＜x,X＞表示局部坐标系x轴与总体坐标系X轴的夹角；

当三维长管由n个直管单元组成时，弯管在总坐标系XYZ下的总刚度矩阵K_t：

其中上角标的i代表直管单元i；表示直管单元1对应的K的分块矩阵形式，表示直管单元i对应的K的分块矩阵形式；之所以将K写成分块矩阵，是由于K为12维的矩阵，为了表示在K_t的矩阵中i+1行对应的位置为直管单元i对应的和直管单元i+1对应的和的形式；

根据非支撑点处和支撑点的节点处直管单元负载向量，计算总体坐标系下总载荷向量F_total：

根据式K_tD_total＝F_total计算得到总变形矢量D_total，总变形矢量即为弯管的最大变形量；

步骤3、以每种弯管支撑位置组合下的弯管最大总变形矢量为适应度函数，通过PSO算法对支撑位置进行优化寻找适应度函数即最大总变形矢量最小的支撑位置，得到最优的支撑位置。

具体实施方式二：

本实施方式所述步骤3的具体过程如下：

PSO粒子群算法是一种模拟鸟类觅食的多目标优化算法；每个可能解都被看做为一个粒子，群体所有粒子信息共享；在本发明中，将弯管的每一种支撑位姿看作一个粒子，包含弯管的支撑位置和弯管的旋转角度；所有粒子都有一个由适应函数fitness()决定的适应度，在本发明中适应度是指每种弯管支撑位置组合下的弯管最大变形量；

首先初始化一群随机粒子，所有粒子组成一个群体，再通过迭代更新每个粒子的速度和位置；粒子个体自身找到的最优解，称为个体最优pbest，可以看做粒子自身的经验；整个群体到目前为止找到的最优解，称为群体最优gbest，可以看做粒子群体的经验；在每一次更新中，粒子通过跟踪个体最优pbest和群体最优gbest来更新自己的状态；

对于第k次迭代，PSO中的第i个粒子按照以下公式进行更新：

其中，i＝1,2,3,...,M；M为粒子群中粒子的总数；是粒子i第k-1次迭代的飞行速度矢量；是粒子i第k-1次迭代后的位置；是第k-1次迭代后，粒子i历史上最优位置；gbest^k是第k-1次迭代后，所有粒子中的最优位置；c1，c2为权重因子；rand()为0到1内的随机数，w为惯性权重。

其他步骤和参数与具体实施方式一相同。

具体实施方式三：

本实施方式所述步骤3的粒子X＝(n₁,n₂,...,n_i,a)，n_i表示第i个支撑选取的位置，a表示弯管的旋转角度。

本方法可以进行任意支撑数量的优化。以一个实际工程为例，当有三个支撑时，并且三维长管可绕一端法兰轴线360度旋转，因此若把弯管的一个支撑位置看作粒子群中的一个粒子如下式：

X＝(n₁，n₂，n₃，a)

支撑位置的节点在计算中处理成固定约束，把相应节点处的载荷向量相应的力和变形矢量相应的位移置零，旋转角度a处理为局部坐标系与总体坐标系的夹角；

通过PSO算法的不断迭代，在一定次数的迭代之后，支撑位置的数据将会收敛，总变形量的值便不再减小，认为到达最优解。

其他步骤和参数与具体实施方式一或二相同。

具体实施方式四：

本实施方式所述步骤2中的总体坐标系下的刚度矩阵K＝RkR^T是基于三维长管的总体坐标系和直管单元对应的局部坐标系的得到的，具体过程如下：

在二次开发工具中，以提取的轴线特征为基础，利用有限元思想推导出在自重状态并有支撑的情况下弯管的有限元模型，在总体坐标系下的总体刚度矩阵K、负载向量F和变形矢量D，关系如下：

KD＝F (1)

在局部坐标系下直管单元的刚度矩阵k、负载向量f和变形矢量d，关系如下：

kd＝f (2)

利用相邻的节点表示一个直管单元；节点直管单元1包括节点1、节点2，直管单元2包括节点2、节点3；依次类推，直管单元i包括节点i、节点(i+1)；

根据材料力学，还可以得到刚度矩阵

局部坐标系下变形矢量d为：

d＝(u₁ V₁ W₁ a_x1 v_y1 v_z1 u₂ V₂ W₂ a_x2 v_y2 v_z2)^T (4)

其中，参数下标中的“1”和“2”表示相邻节点中的节点1和节点2；u₁为节点1的x方向的位移，V₁为节点1的y方向的位移，W₁为节点1的z方向的位移，a_x1为节点1的x方向的转角，v_y1为节点1的y方向的转角，v_z1为节点1的z方向的转角；u₂为节点2的x方向的位移，V₂为节点2的y方向的位移，W₂为节点2的z方向的位移，a_x2为节点2的x方向的转角，v_y2为节点2的y方向的转角，v_z2为节点2的z方向的转角；

局部坐标系下负载向量f为

f＝(F_x1 F_y1 F_z1 M_x1 M_y1 M_z1 F_x2 F_y2 F_z2 M_x2 M_y2 M_z2)^T (5)

其中，F_x1、F_y1、F_z1、M_x1、M_y1、M_z1分别表示节点1在xyz三个方向的力与力矩；F_x2、F_y2、F_z2、M_x2、M_y2、M_z2分别表示节点2在xyz三个方向的力与力矩；

在总体坐标下的直管单元负载向量F为

F＝(F_x1'F_y1'F_z1'M_x1'M_y1'M_z1'F_x2'F_y2'F_z2'M_x2'M_y2'M_z2')^T(6)

其中，F_x1'、F_y1'、F_z1'、M_x1'、M_y1'、M_z1'分别表示节点1在XYZ三个方向的力与力矩；F_x2'、F_y2'、F_z2'、M_x2'、M_y2'、M_z2'分别表示节点2在XYZ三个方向的力与力矩；

在实际情况中，直管单元会与总体坐标系下各个坐标轴有一定的角度；任意直管单元在总体坐标系下的受力情况如图3所示；

根据总体坐标系XYZ与局部坐标系xyz间的空间位置关系可得到总体坐标系下变形矢量D和负载向量F与直管单元局部坐标系下变形矢量d和负载向量f的关系如式：

d＝R^TD (7)

f＝R^TF (8)

其中，公式(7)是说明总体坐标系下变形矢量D与直管单元局部坐标系下变形矢量d的关系，并不能直接利用变形矢量D求解总变形矢量D_total；

根据公式(1)、(2)、(7)、(8)得到总体坐标系下的刚度矩阵K＝RkR^T。

其他步骤和参数与具体实施方式一至三之一相同。

实施例

为了验证利用本发明所建立的弯管有限元模型计算管件变形的准确性，用ANSYS软件对同一管件、同一种支撑位置时的变形量进行仿真，与发明所建立有限元模型的计算结果进行对比。以图1中的管件模型为研究对象，在本发明建立的有限元模型中将管件分为35段进行分析，材料参数如表1，单元杆参数如表2。

表1空间管件的材料属性

Tab.1 Material attribution of space poles

表2单元杆参数表

Tab.2 Unit bar parameters

按照表1设置材料参数，分别选取三种支撑条件下的两个角度，共6种情况，在ANSYS中对管件进行受力分析。管件节点位置及旋转角度示意图如图4所示。

仿真结果和仿真得到的最大变形量如图5(a)～图5(f)所示。其中，图5(a)表示支撑1-15-35、旋转0度的管件变形量分析结果；图5(b)表示支撑1-15-35、旋转90度的管件变形量分析结果；图5(c)为支撑6-15-30、旋转0度的管件变形量分析结果；图5(d)为支撑6-15-30、旋转90度的管件变形量分析结果；图5(e)为支撑10-20-30、旋转0度的管件变形量分析结果；图5(f)为支撑10-20-30、旋转90度的管件变形量分析结果。

将各支撑情况下，ANSYS分析得到的管件的最大变形量与本发明所建立的有限元模型计算结果进行对比，如表3。可以看出，对于用于验证的6种支撑情况，本发明有限元模型计算出的管件最大变形量与ANSYS分析得到的最大变形量的平均偏差为0.07mm，最大偏差为0.13mm。

表3有限元结果与ANSYS仿真结果对比

在整个方法单纯从运算速度上相较传统的有限元分析软件也大大提高，以本发明的实例来说，从开始运行到得到最优运算结果及生成迭代曲线等后处理结果，平均只需要30秒以内，这包括计算了500种支撑情况并从中自动选出最优的支撑位置，而传统有限元软件计算一次支撑情况就需要近1分钟。

以上图1所示的薄壁长管为例，提取轴线，有限元模型的建立已经由上文详述，此处不再赘述。进行优化的过程如下：

粒子群大小取为10，迭代50次。PSO运算迭代过程如图6，得到的最优粒为(8，15，30，60)，即选择8，15，30节点作为支撑点，旋转60度为最优支撑位姿。

从图6中可以看出，迭代次数到达30次左右，总变形量的值的值便不再减小，认为到达最优解。根据有限元算法，计算出在最优支撑位姿下的弯管变形情况如图7。其中，细实线是旋转前的管子轴线；粗虚线是旋转后未变形轴线；粗实线表示转动管子且变形的轴线；实心圆表示支撑点。

变形量在x-y-z方向的分量沿弯管轴线变化曲线如图8(a)至图8(c)所示，图8(a)为Z轴变形分量，图8(b)为X轴变形分量，图8(c)为Y轴变形分量。

Claims (3)

1.三维长管的支撑位置优化方法，其特征在于，包括下述步骤：

步骤1、将三维长管模型输入二次开发工具，提取三维长管模型的轴线；

步骤2、基于有限元法的三维长管刚度矩阵的建立：

在二次开发工具中三维长管的总体坐标系为XYZ，将三维长管分解成若干个直管单元，直管单元对应的局部坐标系为xyz；

局部坐标系下直管单元的刚度矩阵

h₁至h₁₀为中间变量，

式中，E为材料的杨氏模量，A为弯管截面面积，l为单元杆长度，b为刚度修正系数，G为材料的剪切模量，I_y表示截面对y轴的惯性矩、I_z表示截面对z轴的惯性矩，J表示截面的极惯性矩；

利用相邻的节点表示一个直管单元，在总体坐标下的直管单元负载向量

F⁽ⁱ⁾＝(F₁ ⁽ⁱ⁾F₂ ⁽ⁱ⁾)

在非支撑点的节点处，直管单元只受重力产生均布载荷q_i，对应直管单元长度为l_i，非支撑点处的载荷如下

在支撑点的节点处载荷为(0 0 0 0 0 0)^T；

根据局部坐标系下直管单元的刚度矩阵k计算总体坐标系下的刚度矩阵：

K＝RkR^T

其中：R^T表示R的转置矩阵； cos＜·,·＞表示局部坐标系与总体坐标系的夹角；

当三维长管由n个直管单元组成时，弯管在总坐标系XYZ下的总刚度矩阵K_t：

其中上角标的i代表直管单元i；表示直管单元1对应的K的分块矩阵形式，表示直管单元i对应的K的分块矩阵形式；

根据非支撑点处和支撑点的节点处直管单元负载向量，计算总体坐标系下总载荷向量F_total：

根据式K_tD_total＝F_total计算得到总变形矢量D_total，总变形矢量D_total即为弯管的最大变形量；

步骤3、以每种弯管支撑位置组合下的弯管最大总变形矢量为适应度函数，通过PSO算法对支撑位置进行优化寻找适应度函数即最大总变形矢量最小的支撑位置，得到最优的支撑位置；

步骤3的具体过程如下：

将弯管的每一种支撑位姿看作一个粒子，包含弯管的支撑位置和弯管的旋转角度；适应度是指每种弯管支撑位置组合下的弯管最大变形量；

首先初始化一群随机粒子，所有粒子组成一个群体，再通过迭代更新每个粒子的速度和位置；粒子个体自身找到的最优解，称为个体最优pbest；整个群体到目前为止找到的最优解，称为群体最优gbest；在每一次更新中，粒子通过跟踪个体最优pbest和群体最优gbest来更新自己的状态；

对于第k次迭代，PSO中的第i个粒子按照以下公式进行更新：

其中，i＝1,2,3,...,M；M为粒子群中粒子的总数；是粒子i第k-1次迭代的飞行速度矢量；是粒子i第k-1次迭代后的位置；是第k-1次迭代后，粒子i历史上最优位置；gbest^k是第k-1次迭代后，所有粒子中的最优位置；c1,c2为权重因子；rand()为0到1内的随机数，w为惯性权重。

2.根据权利要求1所述的三维长管的支撑位置优化方法，其特征在于，步骤3的粒子X＝(n₁,n₂,...,n_i,a)，n_i表示第i个支撑选取的位置，a表示弯管的旋转角度。

3.根据权利要求1或2所述的三维长管的支撑位置优化方法，其特征在于，步骤2中的总体坐标系下的刚度矩阵K＝RkR^T是基于三维长管的总体坐标系和直管单元对应的局部坐标系的得到的，具体过程如下：

在二次开发工具中，在总体坐标系下的总体刚度矩阵K、负载向量F和变形矢量D，关系如下：

KD＝F (1)

在局部坐标系下直管单元的刚度矩阵k、负载向量f和变形矢量d，关系如下：

kd＝f (2)

利用相邻的节点表示一个直管单元；节点直管单元1包括节点1、节点2，直管单元2包括节点2、节点3；依次类推，直管单元i包括节点i、节点(i+1)；

局部坐标系下变形矢量d为：

d＝(u₁ V₁ W₁ a_x1 v_y1 v_z1 u₂ V₂ W₂ a_x2 v_y2 v_z2)^T (4)

其中，参数下标中的“1”和“2”表示相邻节点中的节点1和节点2；u₁为节点1的x方向的位移，V₁为节点1的y方向的位移，W₁为节点1的z方向的位移，a_x1为节点1的x方向的转角，v_y1为节点1的y方向的转角，v_z1为节点1的z方向的转角；u₂为节点2的x方向的位移，V₂为节点2的y方向的位移，W₂为节点2的z方向的位移，a_x2为节点2的x方向的转角，v_y2为节点2的y方向的转角，v_z2为节点2的z方向的转角；

局部坐标系下负载向量f为

f＝(F_x1 F_y1 F_z1 M_x1 M_y1 M_z1 F_x2 F_y2 F_z2 M_x2 M_y2 M_z2)^T (5)

其中，F_x1、F_y1、F_z1、M_x1、M_y1、M_z1分别表示节点1在xyz三个方向的力与力矩；F_x2、F_y2、F_z2、M_x2、M_y2、M_z2分别表示节点2在xyz三个方向的力与力矩；

在总体坐标下的直管单元负载向量F为

F＝(F_x1'F_y1'F_z1'M_x1'M_y1'M_z1'F_x2'F_y2'F_z2'M_x2'M_y2'M_z2')^T (6)

其中，F_x1'、F_y1'、F_z1'、M_x1'、M_y1'、M_z1'分别表示节点1在XYZ三个方向的力与力矩；F_x2'、F_y2'、F_z2'、M_x2'、M_y2'、M_z2'分别表示节点2在XYZ三个方向的力与力矩；

根据总体坐标系XYZ与局部坐标系xyz间的空间位置关系可得到总体坐标系下变形矢量D和负载向量F与直管单元局部坐标系下变形矢量d和负载向量f的关系如式：

d＝R^TD (7)

f＝R^TF (8)

根据公式(1)、(2)、(7)、(8)得到总体坐标系下的刚度矩阵K＝RkR^T。