CN117057161B - 一种基于数字孪生的风机安装力学分析方法及系统 - Google Patents

Abstract

本发明公开了一种基于数字孪生的风机安装力学分析方法及系统，包括如下步骤：获取杆件结构的基本信息；根据所述杆件结构的基本信息生成杆件单元集合和荷载集合；根据所述杆件单元集合和所述荷载集合集成压缩稀疏行矩阵格式的整体刚度矩阵与外荷载向量；根据所述整体刚度矩阵和所述外荷载向量，采用共轭梯度法获得整体位移向量；根据所述整体位移向量得到变形后的杆件单元集合；根据所述变形后的杆件单元集合获得杆件结构应力分布结果。本发明加速了矩阵的存储和求解效率，且不损失有效信息，从而兼顾计算精度的同时可有效提升计算效率，实现复杂大型结构的快速建模与计算。

Description

一种基于数字孪生的风机安装力学分析方法及系统

技术领域

本发明涉及风机安装领域，特别是涉及一种基于数字孪生的风机安装力学分析方法及系统。

背景技术

风能作为一种清洁的可再生能源越来越受到世界各国的重视。海上风电节约土地、资源量大，使得近海风力发电成为能源发展的重点领域之一，海上风力发电机的相关研究也随之成为海上工程领域的研究热点。为保证风机正常运行，顺利发电，须在设计期间对其构件，主要是杆件的变形与应力情况进行严密地模拟仿真，以评估并完善设计。随着计算机技术的发展，结构分析的计算机方法已日益成为主流。现有的结构力学求解器等计算软件可以完成结构的准确求解，但在面临复杂结构时，其求解效率往往不尽如人意；而当前复杂结构已数见不鲜。但现有技术无法分析风机支架复杂结构的变形与应力。

发明内容

本发明为了解决现有技术无法分析风机支架复杂结构的变形与应力的技术问题，本发明提出了一种基于数字孪生的风机安装力学分析方法及系统。

本发明的技术问题通过以下的技术方案予以解决：

一种基于数字孪生的风机安装力学分析方法，包括如下步骤：

S1：获取杆件结构的基本信息；

S2：根据所述杆件结构的基本信息生成杆件单元集合和荷载集合；

S3：根据所述杆件单元集合和所述荷载集合集成压缩稀疏行矩阵格式的整体刚度矩阵与外荷载向量；

S4：根据所述整体刚度矩阵和所述外荷载向量，采用共轭梯度法获得整体位移向量；

S5：根据所述整体位移向量得到变形后的杆件单元集合；

S6：根据所述变形后的杆件单元集合获得杆件结构应力分布可视化结果。

在一些实施例中，所述杆件结构的基本信息包括：杆件结构中各个杆件单元的物理性能参数，杆件结构中各个杆件单元的摆放位置，杆件单元与杆件单元之间的连接方式以及杆件单元的约束方式，施加在杆件结构节点上的集中荷载信息。

在一些实施例中，所述杆件单元与杆件单元之间的连接方式和所述杆件单元的约束方式采用节点六自由度编码的方式表示，杆端自由度以节点六自由度的形式表示。

在一些实施例中，所述节点六自由度具有X、Y、Z以及RX、RY、RZ六个不同的编号，分别用于表示节点沿X、Y、Z方向的3个平动以及绕X、Y、Z轴的3个转动；

刚性连接为两个节点六自由度的6个编号值完全相同，表明对应的杆端在此处刚性连接；

铰链连接为两个节点六自由度的代表平动的前3个编号值完全相同，代表转动的后3个编号值完全不同，表明对应的杆端在此处铰链连接；

所述节点六自由度满足以下规则：

当所述节点六自由度的某个自由度编号为预设的特殊值时，则表示该自由度编号的自由度被约束；

当两个节点六自由度在某个自由度上具有相同编号时，表明两个节点在此方向上必须具有相同位移或转角。

在一些实施例中，步骤S2中所述生成杆件单元集合的步骤如下：采用结构体系的数值化方法，对杆件结构的所有杆件构造对应的杆件单元，形成杆件单元集合，完成所有所述杆件单元的自由度编号后，所述杆件单元集合保存了杆件结构的所有信息。

在一些实施例中，步骤S3中所述外荷载向量的集成步骤如下：

A1：将所有节点上施加的外荷载分解成沿各个自由度方向的荷载分量；

A2：构造n维零向量F；

A3：将所有所述荷载分量按自由度对应的方式添加到所述n维向量F对应分量上，得到外荷载向量F。

在一些实施例中，步骤S3中所述整体刚度矩阵的集成步骤如下：

B1：计算局部单元刚度矩阵与坐标转换矩阵；

B2：计算全局单元刚度矩阵；

B3：按照自由度编号对应的规则将全局单元刚度矩阵的元素添加到元素全为0的n阶矩阵，得到整体刚度矩阵。

在一些实施例中，步骤S4具体包括：

根据压缩稀疏行矩阵格式的整体刚度矩阵和共轭梯度法求解整体位移线性方程组，得到整体位移向量，所述整体位移向量的表达式如下：

；

其中F为所有外加荷载组成的向量，为待求的所有位移或转角，K为整体刚度矩阵。

在一些实施例中，所述共轭梯度法包括以下步骤：

C1：初始化未知向量为迭代向量；

C2：计算剩余向量和前进方向；

C3：若所述置前进方向为零向量或所述剩余向量的模小于误差限值，则结束迭代，输出所述迭代向量；反之执行步骤C4；

C4：计算前进步长，更新所述迭代向量，计算更新所述迭代向量后的剩余向量；计算共轭系数，更新所述前进方向，之后执行步骤C3。

在一些实施例中，所述压缩稀疏行矩阵格式包括值列表V，列号列表C以及行起始索引列表RS；

其中，值列表V存储了矩阵中所有非0元素的值，满足行号小的元素在前，行号相同的元素则可任意放置；列号列表C与值列表V对应，存储了每个元素所在列的列号；行起始索引列表则按顺序存储了每一行元素在V中的最小索引。

在一些实施例中，步骤S5包括以下步骤：

S51：对每个杆件单元，按照自由度对应的方式从整体位移向量中提取出所述杆件单元的杆端全局位移向量；

S52：计算所述杆件单元的局部位移向量；

S53：将所述杆件单元沿轴向切分为有限多个杆件截面；

S54：计算局部坐标系下每个杆件截面中心的位移值并进行相应移动；

S55：将所述杆件截面重新连接，得到变形后的杆件单元。

在一些实施例中，步骤S6包括以下步骤：

S61：计算每个杆件单元局部坐标系下的杆端力；

S62：将所示杆件截面切分成有限多个扇形；

S63：计算所述扇形每个边缘角点处的应力值；

S64：根据所述应力值对所述杆件结构进行颜色标识，得到所述应力分布可视化结果；

S65：展示所述应力分布可视化结果。

本发明还提出了一种采用上述的基于数字孪生的风机安装力学分析方法的基于数字孪生的风机安装力学分析系统，其特征在于，包括：

结构生成模块：用于获取杆件结构的基本信息，生成杆件单元集合；

荷载管理模块：用于根据所述基本信息生成荷载集合；

结构求解模块：用于根据所述杆件单元集合和所述荷载集合集成整体刚度矩阵与外荷载向量，求解整体位移向量；

变形应力分析模块：用于根据所述整体位移向量计算结构的变形，获得变形后的杆件单元集合，计算结构应力分布，杆件端部的位移及受力，杆件内部指定位置处的位移及应力值均作为分析结果传递给结果展示模块；

结果展示模块：用于将所述分析结果可视化生成应力分布可视化结果，并展示所述应力分布可视化结果。

本发明与现有技术对比的有益效果包括：

本发明通过根据压缩稀疏行矩阵格式的所述整体刚度矩阵和所述外荷载向量，之后采用共轭梯度法求解整体位移向量，加速了矩阵的存储和求解效率，且不损失有效信息，从而兼顾计算精度的同时可有效提升计算效率，实现复杂大型结构的快速建模与计算。同时，本发明对结构体系的表示方法及结构求解方法，可推广应用至风机安装结构之外的任意杆系结构，具有良好的实用性和可扩展性。

本发明实施例中的其他有益效果将在下文中进一步述及。

附图说明

图1是本发明实施例中基于数字孪生的风机安装力学分析方法的流程图。

图2是本发明实施例中杆件单元的局部坐标系示意图。

图3a是本发明实施例中杆件单元在仅有一端侧向位移情况下的轴线变形曲线。

图3b是本发明实施例中杆件单元在仅有一端侧向转角情况下的轴线变形曲线。

图4是本发明实施例中杆件单元外表面应力计算示意图。

图5是本发明实施例中离散应力点分布示意图。

图6是本发明实施例中基于数字孪生的风机安装力学分析系统示意图。

图7是本发明实施例中杆件结构示意图。

图8a是本发明实施例1中第一角度可视化结果示意图。

图8b是本发明实施例1中第二角度可视化结果示意图。

图8c是本发明实施例1中第三角度可视化结果示意图。

图9a是本发明实施例2中第一角度整体可视化结果示意图。

图9b是本发明实施例2中第二角度整体可视化结果示意图。

图9c是本发明实施例2中第二角度局部可视化结果示意图。

具体实施方式

下面对照附图并结合优选的实施方式对本发明作进一步说明。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

需要说明的是，本实施例中的左、右、上、下、顶、底等方位用语，仅是互为相对概念，或是以产品的正常使用状态为参考的，而不应该认为是具有限制性的。

本发明能够实现对杆件组成的复杂结构，包括风机运输过程中支架中变形与应力的模拟仿真，从而解决风机支架在外力荷载下的变形与应力的分析问题。

本发明实施例提出的基于数字孪生的风机安装力学分析方法的原理如下：

风机支架一般属于超静定的杆系结构，其分析方法主要有两类：力法和位移法。位移法以结构结点位移为未知量，利用杆件和结构的刚度关系（即力和位移的关系）求出结点位移，各杆件内力也可同时给出。结构分析的计算机方法中，最常用的是矩阵位移法，也就是采用矩阵运算的位移法。矩阵方程的求解可化为与之等价的变分问题，与传统最速下降法相比，共轭梯度法可以加速一维极小搜索，提高问题求解效率。

本发明实施例提出的基于数字孪生的风机安装力学分析方法如图1所示，包括如下步骤：

S1：获取杆件结构的基本信息。

S2：根据所述杆件结构的基本信息生成杆件单元集合和荷载集合。

S3：根据所述杆件单元集合和所述荷载集合集成压缩稀疏行矩阵格式的整体刚度矩阵与外荷载向量。

S4：根据所述整体刚度矩阵和所述外荷载向量，采用共轭梯度法获得整体位移向量。

S5：根据所述整体位移向量得到变形后的杆件单元集合。

S6：根据所述变形后的杆件单元集合获得杆件结构应力分布结果。

其具体实施操作如下：

S1、获取和确定待分析的杆件结构的基本信息。

其中待分析的杆件结构基本信息包括：结构中各杆件的物理性能参数，主要包括杆件长度L、截面面积A、抗拉刚度EA、抗弯刚度EI、抗扭刚度、抗弯截面系数W；结构中各杆件的摆放位置、连接方式与约束方式，摆放位置主要指杆件两端点的空间坐标，连接方式主要指杆件与杆件之间的刚接、铰接的，约束方式主要指杆件与地面、墙体等固定物体的刚接、铰接等；施加在结构节点上的集中荷载，主要包括荷载类型（力或弯矩）、荷载大小、荷载方向以及施加节点等信息。其中端点指单一杆件单元的端部，节点指两个及以上杆件单元端点连接形成的部位。

本实施例通过空间点、节点六自由度以及单元的概念来实现结构体系的数值化，具体包括如下步骤：

首先，结构体系以杆件单元作为基本结构组成，其主要属性包括2个杆端点、2个杆端自由度、抗拉刚度EA、抗弯刚度EI、抗扭刚度、截面面积A以及抗弯截面系数W。

杆端点以空间点的形式表示，包括在空间直角坐标系下的X、Y、Z坐标；杆端自由度以节点六自由度的形式表示，包括节点沿X、Y、Z方向的3个平动以及绕X、Y、Z轴的3个转动。

结构体系除了杆件单元组成本身，还需要考虑杆件与杆件之间的连接方式与约束关系，例如两根杆件在相交端点处是刚接还是铰接，杆端是否被固定在地面或墙体上等。

本发明实施例主要通过节点六自由度编码的方式来实现连接方式与约束关系的表示，每个节点六自由度都具有X、Y、Z以及RX、RY、RZ六个不同的编号，分别用于表示节点沿X、Y、Z方向的3个平动以及绕X、Y、Z轴的3个转动。

根据上述节点六自由度编码的方式定义如下两个规则：

1）当自由度编号为预设的特殊值时，认为该自由度被约束。

例如，某个节点六自由度编号为，代表平动的三个自由度编号被赋值为FixID（即特殊值），表明其沿3个轴方向的平动被约束；代表转动的三个自由度编号分别被赋值a、b、c（a、b、c均为正整数），表示该节点可以绕3个轴转动（可以理解为与地面铰接）。

2）当两个节点六自由度在某个自由度上具有相同编号时，表明两个节点在此方向上必须具有相同位移或转角。

例如，两个节点六自由度的6个编号值完全相同时，表明对应的杆端在此处刚接（所有位移、转角必须一样）；两个节点六自由度的前3个编号值完全相同、但后3个编号值完全不同时，表明对应的杆端在此处铰接（位移必须一样、转角可以不同）。

按照以上规则，除了表示常见的铰接、刚接方式，也可以表示特殊的滑移铰接、单向铰等连接或约束方式。从以上描述中可以看到，本发明实施例所使用的空间点、节点六自由度以及单元的概念能够将复杂结构进行合理数值化，为结构在计算机中的存储、表示以及变形应力计算打下基础。

S2、根据杆件结构的基本信息生成杆件单元集合和荷载集合。

具体是包括以下两个步骤

1）对步骤S1中的所有杆件构造对应的杆件单元，形成杆件单元集合。

其中，杆件单元的物理参数性能以及2个杆端点在流程（1）中给出，而2个杆端节点六自由度则根据前面给出的自由度设定规则进行确定。

例如，可以用整数作为自由度编号，其中用-1表示，其他自由度从0开始依次赋予到杆件单元中，即：若待编号的自由度被约束，则赋予编号-1；若待编号的自由度与前面已编号的自由度表示同一位移或转角，则直接赋予共用的自由度编号，否则从未使用的自然数中选择最小的一个作为该自由度编号。

完成所有杆件单元的自由度编号后，杆件单元集合即保存了结构体系的所有信息。

S3、根据杆件单元集合和荷载集合集成压缩稀疏行矩阵格式的整体刚度矩阵与外荷载向量，具体按照矩阵位移法集成整体刚度矩阵K与外荷载向量F，其步骤如下：

首先构造一个元素全为0的n阶矩阵K，其中n为编号值不为的自由度个数，并定义好自由度的顺序。

例如，按照从0开始依次赋予自由度的方法，则可直接按编号值从小到大确定自由度的顺序。

对杆件单元集合中的每一个杆件单元执行如下整体刚度矩阵的集成步骤：

B1：计算其局部单元刚度矩阵与坐标转换矩阵/>。

B2：求得其全局单元刚度矩阵。

B3：按照自由度编号对应的规则将的元素添加到整体刚度矩阵K之中，具体是按照全局单元刚度矩阵中各自由度对应的编号，将全局单元刚度矩阵的元素添加到元素全为0的n阶矩阵，得到整体刚度矩阵。其中整体刚度矩阵K的表达式如下：

;

其中，的意义为第j个自由度的单位位移引起的第i个自由度的杆端力。

例如，假设中第3个自由度编号为u，在整体刚度矩阵中自由度排序为第7；第5个自由度编号为v，在整体刚度矩阵中自由度排序为第9，则需将/>矩阵中的（3,3）、（3,5）、（5,3）、（5,5）位置处元素分别添加到矩阵中的（7,7）、（7,9）、（9,7）、（9,9）位置处。对所有杆件单元进行一次该操作，就集成得到了整体刚度矩阵K。

同理，对所有节点上施加的外荷载向量的集成步骤如下：

A1：将所有节点上施加的外荷载分解成沿各个自由度方向的荷载分量。

A2：构造一个元素全为0的n维向量F。

A3：将所有荷载分量按自由度对应的方式添加到n维向量F对应分量上，就集成得到了外荷载向量F。

其中矩阵位移法的基本思路如下：

以所有自由度的位移或转角为未知量，通过求解外加荷载与位移/转角形成的线性方程组来得到位移或转角，可写成如下矩阵形式：

；

其中F为所有外加荷载组成的向量，为待求的所有位移或转角，K为整体刚度矩阵。另外，对于静定或超静定结构，整体刚度矩阵均满足对称正定性质。

要应用矩阵位移法，需要先集成外荷载向量F与整体刚度矩阵K。

（1）对于外荷载向量F，主要考虑节点集中荷载。

首先将所有自由度按照某种规则（例如编号值大小）进行排序，然后将所有节点荷载按分量添加到外荷载向量上即可。

例如，某超静定结构共有8个自由度，编号值为0~7，按顺序形成的编号向量为，某一杆端节点六自由度编号为/>，在该节点上施加了XY平面内45°方向、大小1N的外加力，则沿自由度编号3和4方向的荷载分量均为/>，将其添加到整体荷载向量的相应编号位置即可。

（2）对于整体刚度矩阵K，还需要考虑局部坐标系与全局坐标系的转换问题。

此处的局部坐标系是指以杆件单元轴线方向为X向的直角坐标系，不同单元通常具有不同的局部坐标系；全局坐标系则是指用于描述整体结构的坐标系，也用于确定自由度的方向（例如节点六自由度的X平动是节点沿全局坐标系X方向的运动）。对于杆件单元，局部坐标系下的荷载与位移存在如下关系：

其中为12阶方阵，被称为局部单元刚度矩阵，该矩阵与单元的性质有关，可表示如下，式中L为杆件的长度。

=

；

另外，局部坐标系下的单元荷载、位移/>与全局坐标系下的单元荷载/>、位移存在如下转换关系：

；

；

其中为坐标转换矩阵，与杆件单元的轴向有关，也是一个12阶方阵。杆件轴向可用单位方向向量（从杆件起始端点指向结束端点、长度为1的向量）表示为

其中，分别为单位方向向量/>在X，Y，Z三个轴方向的坐标。在球坐标系下，可以进一步用经度Lat、纬度Lon描述该方向，满足

；

；

若为0，则认为Lon=0；否则

；

在此表示法下，坐标转换矩阵可记为

；

；

该坐标转换矩阵为正交阵，且满足。进一步有

；

在全局坐标系下，有，其中/>为全局单元刚度矩阵，故有

；

因此，已知杆件单元的坐标转换矩阵和局部单元刚度矩阵，利用上式即可得到全局单元刚度矩阵，此后将对应位置元素添加到整体刚度矩阵K中即可。

例如，某超静定结构共有8个自由度，编号值为0~7，按顺序形成的编号向量为，某一单元的两端节点六自由度编号为与/>，则将该单元的全局单元刚度矩阵中左上角2×2子矩阵添加到整体刚度矩阵中3、4编号处的子矩阵中即可。

按照上述方法集成得到外荷载向量F与整体刚度矩阵K后，求解形成的线性方程组即可得到位移向量。

S4、根据整体刚度矩阵和外荷载向量，采用共轭梯度法求解获得整体位移向量，具体将整体刚度矩阵K转化为CSR格式（或在前一步集成时就采用CSR格式），利用共轭梯度法求解线性方程组得到整体位移向量/>。其中，初始迭代向量可取为n维零向量，误差限值limit通常可按相对误差进行设定，例如取为外荷载向量模长的一亿分之一，即/>。

本发明实施例中的共轭梯度法包括如下步骤：

对于具有大量杆件组成的复杂结构体系，上述整体刚度矩阵通常具有大型稀疏的特点，即矩阵的阶数很大，并且大部分元素都是0。对于这类矩阵，虽然能够采用常规的高斯消去法进行求解，但是过程中无法利用矩阵的稀疏性，这将导致十分巨大的计算量。实际上稀疏矩阵通常采用迭代法进行求解，本发明采用适合对称正定矩阵的共轭梯度法来求解上述线性方程组。

对于待求解的线性方程组，其中A为对称正定矩阵，x为未知向量，b为与x维数相同的常数向量，在给定误差限值limit下共轭梯度法的计算流程如下：

C1：初始化向量x为任意向量，其中上标表示迭代次数，（0）即初始值。

C2：计算剩余向量。

C3：若为0向量或/>（即误差限值），则结束迭代，输出所述迭代向量即为满足误差限值的近似解；否则，继续执行以下步骤。其中k=0,1，...，表示迭代次数；此处误差限值的取值可根据实际所需的计算精度自行选择设置。

C4：计算前进步长，更新向量/>，计算新的剩余向量/>

计算共轭系数，更新/>

完成一次迭代，返回到判断步骤C3。

共轭梯度法具有所需存储量小、收敛速度较快、稳定性高等优点，是解大型线性方程组的高效算法之一。

由于共轭梯度法在迭代过程中不会更改系数矩阵A，因此能够保持其稀疏性。提高求解效率的另一个要点是避免矩阵中0元素的运算，这一思路可以通过利用稀疏矩阵结构进行实现。

本发明实施例采用CSR（压缩稀疏行）矩阵格式来进行系数矩阵的存储，CSR格式仅存储矩阵中的非零元素值及其所在位置的列号。具体的，CSR格式用3个顺序列表存储矩阵信息，分别为值列表V、列号列表C以及行起始索引列表RS。

其中，值列表V存储了矩阵中所有非0元素的值，满足行号小的元素在前，行号相同的元素则可任意放置；列号列表C与值列表V对应，存储了每个元素所在列的列号；行起始索引列表则按顺序存储了每一行元素在V中的最小索引。

下面给出了一个稀疏矩阵的表示格式（索引从0开始）：

;

例如，第3行只有1个元素0.5，处在第2列，列索引号为1，故V中0.5在C中对应的元素为1；第3行前面共有3个元素，故值列表V中首个第3行元素的索引号为3，RS中第3个元素即为3。因为CSR格式只对稀疏矩阵中的非零元素进行存储，所以能够在极大程度上减少矩阵占用的存储空间。

另一方面，采用CSR格式便于进行稀疏矩阵与向量的乘法，矩阵A乘以向量p将得到一个新向量，新向量的每个元素实际上是矩阵每一行与向量p的点积。即

；

例如，用CSR矩阵的索引号为dex的行点乘向量p时，按下式计算即可

；

其中，表示/>在V中的最小索引，由于k的取值范围为/>到/>，故/>即为/>中的非零元素。/>为/>对应的列号，/>即为p中与/>位置相应的元素。可以看到，该向量点乘中仅考虑/>中的非零元素，计算时间取决于稀疏矩阵中非零元素的个数而不是所有元素的个数，因此具有极高的效率。此外，共轭梯度法中的系数矩阵A只涉及到与向量p或x的乘法，而不涉及其他矩阵相乘等操作，因此能够很好适配CSR矩阵与向量相乘高效这一特点。

S5、计算杆件结构变形，根据整体位移向量得到变形后的杆件单元集合；

S51：按照单元位移应力函数对每个杆件单元按照自由度对应的方式，从整体位移向量中提取出单元杆端的全局位移向量。

S52：计算求出杆件单元的局部位移向量。

S53：杆件结构中杆件的变形通常可采用离散的方式进行表示，具体通过将杆件单元沿轴向切分为有限多个杆件截面。

S54：计算局部坐标系下每个杆件截面中心的位移并进行相应移动。

S55：将杆件截面重新连接为变形后的曲杆（即变形后的杆件单元）即可。

例如，将杆件沿轴向均匀分成10份，各截面的局部x坐标分别为0、01L、0.2L……L，对应的比例值即分别为0、0.1、0.2……1，根据局部位移向量以及下式可以很快求出局部坐标系下各个截面处的位移值，其表达式如下：

;

最后将截面重新组合为杆件后就得到了变形后的结构。

以上内容给出了求解整体位移向量的方法，得到的仅为结构上杆端在每个自由度上的位移或转角值，例如杆端沿全局x方向的位移量，而无法给出杆件单元内部的位移与应力等，因此需要使用单元位移应力函数对解得的位移向量/>做进一步处理。

首先，对于指定单元，按照自由度编号对应规则，从整体位移向量中提取出单元两端的12个位移分量（被约束的分量置为0）形成全局坐标系下的单元位移向量，然后利用坐标转换矩阵按下式计算局部坐标系下的单元位移表达式如下：

；

杆件单元的局部坐标系如图2所示，此时单元位移可表示为如下公式：

；

其中表示线位移，/>表示角位移；x/y/z表示该位移沿或绕X/Y/Z轴方向；下标a/b表示该位移发生于杆件单元的a/b端。如/>即表杆件a端沿b轴方向的线位移。设杆件单元上任意一点M的坐标为/>的，则M的局部位移仅与/>，L为杆件单元的长度。则点M的空间位移/>与杆端局部位移关系如下：/>

;

其中分别为点M沿X，Y，Z轴方向的空间位移分量。按上式即可得到杆件单元内部任意一点的位移，例如杆件单元在仅有一端侧向位移情况下的轴线变形曲线如图3a所示，仅有一端发生转角情况下的轴线变形曲线如图3b所示。

得到杆件单元在局部坐标系下的杆端位移之后，还可以利用局部单元刚度矩阵按下式计算局部坐标系下的杆端力表达式如下：

；

其中杆端力可写成如下形式：

；

其中代表N轴力，V代表剪力，M代表弯矩；x/y/z表示该力沿或绕X/Y/Z轴方向；下标a/b表示该力作用于杆件单元的a/b端。在仅有杆端节点荷载的情况下，杆端力将满足下述条件：

在此条件下，杆件单元上任意局部x坐标截面下，记，可知截面处的轴力与/>相同（受压为正），绕局部y方向和z方向的弯矩可按如下公式计算：

在局部x坐标截面下，杆件单元外表面上任意一点的拉应力可按下式计算，W为前面提到的单元抗弯截面系数，角度/>的含义如图4所示，拉应力/>的表达式如下：

；

S6、根据所述变形后的杆件单元集合计算获得杆件结构应力分布结果，与步骤S5一致，按照每个杆件单元对应的自由度编号，从整体位移向量中提取出单元杆端的全局位移向量，之后包括如下步骤：

S61：计算求出其每个杆件单元局部坐标系下的杆端力。

S62：同样采用离散方式进行应力表示，沿轴向切分出有限多个截面，然后将每个圆截面切分成有限多个扇形。

S63：按下式计算扇形每个边缘角点处的应力值：

计算完毕后即得到杆件表面一些离散点处的应力值，应力点分布如图5所示。

S64：通常应力值大小可通过颜色对杆件结构进行颜色标识，得到所述应力分布可视化结果；

例如应力为正的区域用红色标识、应力为负的区域用蓝色标识，且应力绝对值越大的地方颜色值越深。在上述离散应力点的基础上，按照应力值利用线性插值函数对结构表面进行上色即可展示结构的应力分布。

S65：展示所述应力分布可视化结果。

本发明实施例基于上述风机安装杆件结构力学分析的方法，实现了可用于高效求解复杂杆系结构变形与内力的杆件结构分析的基于数字孪生的风机安装力学分析系统。该系统由如下模块组成：

基于数字孪生的风机安装力学分析系统可划分为以下5个模块：结构生成模块、荷载管理模块、结构求解模块、变形应力分析模块、结果展示模块。各个模块之间的关系如图6所示，其功能分别如下：

（1）结构生成模块：

其用于生成数值化的结构体系，主要实现杆件单元集合及自由度编号的管理功能。其中，杆件单元的节点六自由度利用装置中的DOF6对象进行表示，而所有自由度通过DOFManager对象进行管理，在DOFManager对象中内置了记录已分配自由度数量的变量Count，可以通过New方法自动分配一个未被使用的自由度编号，通过New6方法返回一个新的DOF6对象，该对象中的6个自由度编号均未被使用，此外还可通过Fix或Fix6方法返回被约束的自由度。每根杆件单元用Element对象表示，包含了杆件单元的2个杆端点、2个杆端自由度、抗拉刚度EA、抗弯刚度EI、抗扭刚度、截面面积A以及抗弯截面系数W。杆件单元集合通过装置中的Structure对象进行管理，主要通过AddElement方法添加杆件单元。Structure对象中内置了一个DOFManager对象用于结构自由度的管理。

（2）荷载管理模块：

其用于管理作用在结构节点上的荷载集合。单个荷载通过DOFLoad对象存储，包含作用的自由度编号以及荷载值。荷载集合使用LoadCase对象进行管理，主要通过AddDOFLoad方法添加单个荷载。另外，LoadCase对象还支持按系数添加功能，以支持作用于结构的不同荷载工况组合。在Structure对象中内置了一个默认的LoadCase对象用于记录临时工况。

（3）结构求解模块：

其用于矩阵位移法中线性方程组的构造与求解。基于整体刚度矩阵和外荷载向量的集成方法与步骤，从结构生成模块中的Structure对象中获取所有杆件单元信息生成整体刚度矩阵，从荷载管理模块中的LoadCase对象中获取所有荷载信息生成外荷载向量。基于共轭梯度法和CSR矩阵，求解得到包含所有自由度对应位移或转角的整体位移向量。

（4）变形应力分析模块：

其用于计算结构的变形以及各杆件上的应力分布。对于每一杆件，根据其端部对应的节点编号，从整体位移向量中提取出相应编号的元素，即为杆件的端部位移，按前述方法与步骤计算出杆件内部给定位置处的位移值，得到变形后的杆件形状，进一步计算杆件的端部受力以及内部给定位置处的应力值。杆件端部的位移及受力，杆件内部指定位置处的位移及应力值均作为分析结果传递给结果展示模块。

（5）结果展示模块：

其用于分析结果的可视化。利用三维可视化渲染引擎，将变形后的结构以三维立体模型的方式绘制到屏幕上，直观展现各杆件的变形情况；通过对结构杆件表面按应力值进行着色，直观给出结构不同位置处的应力大小。

本实施例通过C#编程语言实现，一些具体实现细节如下：

（1）关于杆件截面属性的描述如下：

本实施例利用Section对象实现常见杆件截面形状的属性自动计算，例如调用静态的Ring方法可生成环形杆件截面，只需传入材料的弹性模量E、剪切模量G以及外半径R以及内半径r即可，相关截面属性将按以下公式计算：

截面面积A：

截面惯性矩I：

抗弯截面系数W：

截面极惯性矩：

；

（2）基于SharpDX的三维可视化：

本实施例利用SharpDX（DirectX的封装库）进行三维模型的绘制与显示，其中杆件三维体通过表面模型进行绘制，即将杆件表面划分成一系列空间小三角形，再根据应力值对三角形进行逐个渲染，因此在变形应力分析模块中只需给出三角形角点处的离散位移值与应力值即可。

（3）基于预处理加快求解速度：

当整体刚度矩阵病态时，共轭梯度法的收敛较慢。为此，可以先通过预处理方法降低系数矩阵的条件数，改善其收敛性。具体的，选择一个与整体刚度矩阵K同阶的非奇异方阵S，设向量，即u满足/>，从而可将原方程/>改写为如下式子：

；

通过选择合适的S，可以让的条件数比起K有所改善，且/>仍保持对称正定，从而更快地解方程得到u，之后代入上式即得到整体位移向量/>。/>

实施例

如图7所示的待分析杆件结构，图中A，B，C，D指示各杆件端部，数字编号则为图中各节点位移相应的自由度编号。杆件的属性均为，各杆件端点之间、端点与地面之间均刚接，且YZ平面外的转角、位移均被约束，欲求在B端添加水平向右大小为1N的力时的结构变形应力。

在步骤S1确定所需的待分析结构的基本信息后，在结构生成模块生成杆件单元集合，构建待求解结构；在荷载管理模块确定荷载工况（即荷载信息）。其中步骤S1的具体代码如下：

1.double EA = 1, EI = 1; //杆件单元属性值。

2.var st = new Structure3D(); //定义待分析结构实例。

3.var gen = st.DOFManager; //用于管理结构自由度。

4.var p1 = new Point3(0, 0, 0); //A点坐标。

5.var p2 = new Point3(0, 0, 1); //B点坐标。

6.var p3 = new Point3(0, 1, 1); //C点坐标。

7.var p4 = new Point3(0, 1, 0); //D点坐标。

8.var dof1 = gen.Fix6(); //A点6个自由度被完全约束。

9.var dof2 = new DOF6(y: gen.New(), z: gen.New(), rx: gen.New()); //B点具有3个自由度。

10.var dof3 = new DOF6(y: gen.New(), z: gen.New(), rx: gen.New());//C点具有3个自由度。

11.var dof4 = gen.Fix6(); //D点具有6个自由度。

12.st.AddElement(p1, p2, dof1, dof2, EA, EI); //向结构中添加杆件AB。

13.st.AddElement(p2, p3, dof2, dof3, EA, EI); //向结构中添加杆件BC。

14.st.AddElement(p3, p4, dof3, dof4, EA, EI); //向结构中添加杆件CD。

15.st.AddDOFLoad(dof2.Y, 1); //向结构中添加节点荷载。

S2：根据从结构生成模块获取的杆件单元信息和从荷载管理模块获取的荷载信息，之后执行步骤S3，在结构求解模块中集成的整体刚度矩阵K和外荷载向量F的表达式如下：

S4：应用共轭梯度法，取初始迭代向量为零向量，取误差限值为，迭代得到整体位移向量的表达式如下：/>

；

S5：在变形应力分析模块，根据整体位移向量，进一步求出杆件结构变形以及杆件结构上的应力分布结果。

S6：将变形及应力分布结果传递给结果展示模块，所实现的第一角度可视化结果示意图如图8a所示，第二角度可视化结果示意图如图8b所示，第三角度可视化结果示意图如图8c所示，其中灰色无颜色变化的是未承受荷载时的杆件，与之对照的是承受荷载后的杆件，两者间的偏离程度表示杆件承受荷载后的变形大小，杆件颜色深浅表示相应位置应力的大小。

实施例

使用本发明实施例提出的基于数字孪生的风机安装力学分析方法及系统对项目中设计的风机结构进行分析。该杆件结构由490根杆件构成，各杆件端点之间、端点与地面之间均刚接，共有1026个自由度。

S1：确定所需的待分析结构的基本信息后，执行步骤S2，在结构生成模块生成杆件单元集合，构建待求解结构；在荷载管理模块确定荷载工况（即荷载信息）。

S3：根据从结构生成模块获取的杆件单元信息和从荷载管理模块获取的荷载信息，在结构求解模块中集成的整体刚度矩阵K和外荷载向量F。

S4：应用共轭梯度法，取初始迭代向量为零向量，取误差限值为，迭代得到整体位移向量/>。

S5：在变形应力分析模块，根据整体位移向量，进一步求出杆件结构变形以及杆件结构上的应力分布结果。

S6：将变形及应力分布传递给结果展示模块，所实现的第一角度整体可视化结果示意图如图9a所示，第二角度整体可视化结果示意图如图9b所示，第二角度局部可视化结果示意图如图9c所示，其中杆件单元的弯曲、倾斜等即表示杆件承受荷载后的变形，杆件颜色深浅表示相应位置应力的大小。

对于上述两个实施例，实施例1结构简单，同其他方法相比，难以体现出明显的效率优势。但在求解实施例2的复杂结构时，采用本发明实施所提出的基于数字孪生的风机安装力学分析方法，即使用共轭梯度法并结合CSR矩阵格式，求解时间不到1秒。而如果采用高斯消去法等方法，同时使用传统的矩阵格式，则在获得相同求解精度的前提下，该杆件结构求解需耗时20秒左右。且可见，本发明实施例中的基于数字孪生的风机安装力学分析方法可在兼顾计算精度的同时有效提升计算效率。

与传统的高斯消去法等方法在使用传统的矩阵格式相比，本发明实施例中的基于数字孪生的风机安装力学分析方法，采用了共轭梯度法进行方程迭代求解，可以有效降低迭代次数，而并不降低求解精度；使用CSR（压缩稀疏行）矩阵格式来进行系数矩阵的存储，仅存储矩阵中的非零元素值及其所在位置，加速了矩阵的存储和求解效率，但并不损失有效信息。因此本发明可在兼顾计算精度的同时可有效提升计算效率，实现复杂大型结构的快速建模与计算。同时，本发明对结构体系的表示方法及结构求解方法，可推广应用至风机安装结构之外的任意杆系结构，具有良好的实用性和可扩展性。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干等同替代或明显变型，而且性能或用途相同，都应当视为属于本发明的保护范围。

Claims (7)

1.一种基于数字孪生的风机安装力学分析方法，其特征在于，包括如下步骤：

S1：获取杆件结构的基本信息；

S2：根据所述杆件结构的基本信息生成杆件单元集合和荷载集合；

S3：根据所述杆件单元集合和所述荷载集合集成压缩稀疏行矩阵格式的整体刚度矩阵与外荷载向量；

S4：根据所述整体刚度矩阵和所述外荷载向量，采用共轭梯度法获得整体位移向量；

S5：根据所述整体位移向量得到变形后的杆件单元集合；

S6：根据所述变形后的杆件单元集合获得杆件结构应力分布可视化结果；

其中，所述杆件结构的基本信息包括：杆件结构中各个杆件单元的物理性能参数，杆件结构中各个杆件单元的摆放位置，杆件单元与杆件单元之间的连接方式以及杆件单元的约束方式，施加在杆件结构节点上的集中荷载信息；

步骤S2中，所述杆件单元与杆件单元之间的连接方式和所述杆件单元的约束方式采用节点六自由度编码的方式表示，杆端自由度以节点六自由度的形式表示，用整数作为自由度编号，所述节点六自由度具有X、Y、Z以及RX、RY、RZ六个不同的编号，分别用于表示节点沿X、Y、Z方向的3个平动以及绕X、Y、Z轴的3个转动，其中被约束的自由度的自由度编号为，用-1表示，其他自由度从0开始依次赋予到杆件单元中，即：若待编号的自由度被约束，则赋予编号-1；若待编号的自由度与前面已编号的自由度表示同一位移或转角，则直接赋予共用的自由度编号，否则从未使用的自然数中选择最小的一个作为该自由度编号；

步骤S3中所述外荷载向量的集成步骤如下：

A1：将所有节点上施加的外荷载分解成沿各个自由度方向的荷载分量；

A2：构造n维零向量；

A3：将所有所述荷载分量按自由度对应的方式添加到所述n维零向量对应分量上，得到外荷载向量F；

步骤S3中，具体按照矩阵位移法集成整体刚度矩阵K与外荷载向量F，其步骤如下：

首先构造一个元素全为0的n阶矩阵K，其中n为编号值不为的自由度个数，并定义好自由度的顺序；

按照从0开始依次赋予自由度的方法，直接按编号值从小到大确定自由度的顺序；

对杆件单元集合中的每一个杆件单元执行如下整体刚度矩阵的集成步骤如下：

B1：计算其局部单元刚度矩阵与坐标转换矩阵/>；

B2：求得其全局单元刚度矩阵；

B3：按照自由度编号对应的规则将的元素添加到整体刚度矩阵K之中，具体是按照全局单元刚度矩阵中各自由度对应的编号，将全局单元刚度矩阵的元素添加到元素全为0的n阶矩阵，得到整体刚度矩阵；

步骤S4具体包括：

根据压缩稀疏行矩阵格式的整体刚度矩阵和共轭梯度法求解整体位移线性方程组，得到整体位移向量，求解所述整体位移向量的表达式如下：

；

其中F为所有外加荷载组成的外荷载向量，为待求的整体位移向量的所有位移或转角，K为整体刚度矩阵；

所述共轭梯度法包括以下步骤：

C1：初始化未知向量为迭代向量；

C2：计算剩余向量和前进方向；

C3：若所述前进方向为零向量或所述剩余向量的模小于误差限值，则结束迭代，输出所述迭代向量；反之执行步骤C4；

C4：计算前进步长，更新所述迭代向量，计算更新所述迭代向量后的剩余向量；计算共轭系数，更新所述前进方向，之后执行步骤C3。

2.如权利要求1所述的基于数字孪生的风机安装力学分析方法，其特征在于，

刚性连接为两个节点六自由度的6个编号值完全相同，表明对应的杆端在此处刚性连接；

铰链连接为两个节点六自由度的代表平动的前3个编号值完全相同，代表转动的后3个编号值完全不同，表明对应的杆端在此处铰链连接；

所述节点六自由度满足以下规则：

当所述节点六自由度的某个自由度编号为预设的特殊值时，则表示该自由度编号的自由度被约束；

当两个节点六自由度在某个自由度上具有相同编号时，表明两个节点在此方向上必须具有相同位移或转角。

3.如权利要求2所述的基于数字孪生的风机安装力学分析方法，其特征在于，步骤S2中所述生成杆件单元集合的步骤如下：采用结构体系的数值化方法，对杆件结构的所有杆件构造对应的杆件单元，形成杆件单元集合，完成所有所述杆件单元的自由度编号后，所述杆件单元集合保存了杆件结构的所有信息。

4.如权利要求1所述的基于数字孪生的风机安装力学分析方法，其特征在于，所述压缩稀疏行矩阵格式包括值列表V，列号列表C以及行起始索引列表RS；

其中，值列表V存储了矩阵中所有非0元素的值，满足行号小的元素在前，行号相同的元素则可任意放置；列号列表C与值列表V对应，存储了每个元素所在列的列号；行起始索引列表则按顺序存储了每一行元素在V中的最小索引。

5.如权利要求2所述的基于数字孪生的风机安装力学分析方法，其特征在于，步骤S5包括以下步骤：

S51：对每个杆件单元，按照自由度对应的方式从整体位移向量中提取出所述杆件单元的杆端全局位移向量；

S52：计算所述杆件单元的局部位移向量；

S53：将所述杆件单元沿轴向切分为有限多个杆件截面；

S54：计算局部坐标系下每个杆件截面中心的位移值并进行相应移动；

S55：将所述杆件截面重新连接，得到变形后的杆件单元。

6.如权利要求5所述的基于数字孪生的风机安装力学分析方法，其特征在于，步骤S6包括以下步骤：

S61：计算每个杆件单元局部坐标系下的杆端力；

S62：将所示杆件截面切分成有限多个扇形；

S63：计算所述扇形每个边缘角点处的应力值；

S64：根据所述应力值对所述杆件结构进行颜色标识，得到所述应力分布可视化结果；

S65：展示所述应力分布可视化结果。

7.一种采用如权利要求1-6任一项所述的基于数字孪生的风机安装力学分析方法的基于数字孪生的风机安装力学分析系统，其特征在于，包括：

结构生成模块：用于获取杆件结构的基本信息，生成杆件单元集合；

荷载管理模块：用于根据所述基本信息生成荷载集合；

结构求解模块：用于根据所述杆件单元集合和所述荷载集合集成整体刚度矩阵与外荷载向量，求解整体位移向量；

变形应力分析模块：用于根据所述整体位移向量计算结构的变形，获得变形后的杆件单元集合，计算结构应力分布，杆件端部的位移及受力，杆件内部指定位置处的位移及应力值均作为分析结果传递给结果展示模块；

结果展示模块：用于将所述分析结果可视化生成应力分布可视化结果，并展示所述应力分布可视化结果。