一种多位移约束下的敏度求解方法与流程

本发明属于飞机结构设计领域，涉及一种飞机结构规划法优化设计中的多位移约束下的敏度求解方法。

背景技术：

在飞机结构优化设计中，基于变量敏度信息的规划法是解决考虑位移、扭角等全局约束优化问题的有效手段。采用解析法和半解析法进行位移约束敏度求解时需要对每个位移约束计算其对应虚位移，由于飞机结构分析需要模拟多种飞行情况，导致分析模型工况及边界条件多种多样。而在对如机翼等结构进行考虑结构变形连续性的优化时，常常需要批量施加位移约束，多工况多约束使得位移敏度求解直接影响了规划法优化的效率和精度。现有优化工具在位移敏度求解时大都只包含一套实工况，即使有多套工况，也只取某一套工况的边界条件来计算所有虚位移。此外，在虚位移求解时需要求解静力平衡方程，这个过程中需要按照边界条件对结构刚度矩阵进行约简，当多个虚工况采用相同边界条件时，约简后的刚度矩阵也相同，如果不考虑约简后刚度矩阵的重用，而是对每一个工况进行一次单独的刚度矩阵约简，随着模型规模和约束规模的增加，计算机不可接受。上述缺点严重制约了基于敏度信息的规划法在解决大规模优化问题中的应用。

技术实现要素：

本发明的目的在于提供一种多位移约束下的敏度求解方法，解决不同工况下多位移约束的敏度求解的问题。

本发明的目的通过如下技术方案实现：

步骤一，构建飞机结构优化模型，获取所述优化模型的工况、变量、边界条件、位移约束及位移约束施加的位移约束点的信息；

步骤二，对所述优化模型进行静力分析，获取结构的位移及结构的刚度矩阵K；

步骤三，找出使各所述位移约束点产生最大位移的所述工况，同时，将找到的使各所述位移约束点产生最大位移的所述工况按所述边界条件进行分类，将具有相同所述边界条件的所述工况划分在同一组，按照选取的所述位移约束点的顺序排列，形成实位移列阵U；

步骤四，对所述刚度矩阵K进行约简，消除所述刚度矩阵K的奇异性，形成Ki矩阵；

步骤五，在所述Ki矩阵中，选取与步骤三中相同的所述位移约束点，在所述位移约束点上施加单位载荷，形成虚载荷列阵Pri；

步骤六，计算得到不同所述位移约束及不同所述边界条件下的结构虚位移列阵Upsei，并将多个所述结构虚位移列阵Upsei构成虚位移矩阵Upse；

步骤七，计算所述变量的单元微分刚度矩阵Kde，并根据所述单元微分刚度矩阵Kde装配结构微分刚度矩阵Kd；

步骤八，依次提取所述实位移列阵U的各列并转置，与所述结构微分刚度矩阵Kd及所述虚位移列阵Upse的对应列相乘，得到当前变量的各位移约束敏度。

优选地是，所述步骤四中形成所述Ki矩阵能够按照步骤三中不同类型的所述边界条件的个数，作为约简所述刚度矩阵K的循环次数，对所述刚度矩阵K进行约简，形成所述Ki矩阵。

优选地是，所述步骤六中结构虚位移列阵Upsei根据求解平衡方程Ki*Upsei＝Pri得到。

优选地是，所述步骤八中当前变量的各位移约束敏度的求解公式为：其中，为第i个变量中的第r个位移约束的敏度，为实位移矩阵U中第r列的转置，为结构微分刚度矩阵Kd，为虚位移矩列阵Upse的第r列。

本发明所提供的多位移约束下的敏度求解方法的有益效果在于，不需要事先对虚工况的边界做假设，完全尊重工程实际，静力平衡方程求解时充分考虑了数据的重用性，提升了计算机效率和准确性，使得基于敏度信息的规划法能应用于大规模结构优化问题中。

附图说明

图1为本发明多位移约束下的敏度求解方法中敏度计算的详细流程图；

具体实施方式

为使本发明实施的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行更加详细的描述。在附图中，自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。所描述的实施例是本发明一部分实施例，而不是全部的实施例。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面结合附图对本发明的多位移约束下的敏度求解方法做进一步详细说明。

本发明的一种起落架航向刚度模拟方法，包括如下步骤：

步骤一，构建飞机结构优化模型，获取所述优化模型的工况、变量、边界条件、位移约束信息(含位移约束及位移约束施加的位移约束点)。

工况为实际工作状况，即所处的各种工作环境。

其中，变量可以是金属的厚度、杆和梁的截面面积、以及复合材料的总厚度和分层厚度等；

结构在不同的工作环境下会有不同的边界条件，例如一根杆，约束一头的节点与约束中间的节点就会形成两种边界条件；

位移约束信息是指优化中位移约束的定义，包括位移约束的施加节点编号和约束自由度(X、Y、Z)，其中，在三维空间一个节点有6个自由度，但对某一个特定节点位移约束只考虑3个平动自由度(X、Y、Z)，因此，对一个节点，最多只能施加3个位移约束，但是一个优化问题可以在不同节点施加位移约束，本发明方法理论上对施加位移约束的总个数没有限制，即位移约束个数等于节点个数乘以自由度数。需要说明的是位移约束点与节点含义相同，只是表述不同而已，能够替换。

步骤二，对飞机结构优化模型进行静力分析，获取结构的位移及结构的刚度矩阵K，其中，根据实际定义的工况对飞机结构优化模型进行一次完整的静力分析，得到结构所有实际工况下的位移结果。

假设优化模型包含n个工况，Sub₁、Sub₂、…、Sub_n，则Us可以表述为：则第一列是Sub₁工况的实位移列阵U₁，优化模型各节点的位移分别为a11、a21、…、ak1，第二列是Sub₂工况的实位移列阵U ₂，以此类推到Sub_n工况，由n个位移列阵U组成了实位移矩阵Us，即所有实际工况下的位移结果。(注：步骤二还未涉及位移约束信息，只是完整的静力分析得到的实际工况下的位移结果。)

刚度矩阵K是方阵，其行数和列数都与Us的行数完全一致，均表示结构所有节点的自由度数，即位移约束个数。

步骤三，找出使各位移约束点产生最大位移的工况，同时，将找到的使各位移约束点产生最大位移的工况按边界条件进行分类，将具有相同边界条件的工况划分在同一组，按照选取的位移约束点的顺序排列，形成实位移列阵U。

假设优化模型有n个工况，k个位移约束，分别挑选出各位移约束点最大位移所在的工况C₁，C₂，......，C_k，其中多个工况号可能重复，因此挑选出的工况一定是初始所有工况的子集。将挑选出的各工况按边界条件共分为m组，S₁，S₂，……S_m，每组包含相同的边界条件，按照这m组中各工况对应的位移约束顺序提取每个约束对应的实位移并按列存储到实位移列阵U，则实位移列阵U的列数为位移约束个数，各列位移值与数组S中对应的位移约束排序必须保持一致。

对如何挑选出各位移约束点最大位移所在的工况进行详细说明。假设优化模型有n个工况，k个位移约束，第一个位移约束点为2号节点，从Us知道其对应最大位移出现在Sub1，则挑选出的位移列阵为Us的第一列，第二个位移约束点为5号节点，从Us知道其对应最大位移也出现在Sub1，则挑选出的位移列阵仍旧为Us的第一列，如果第二个位移约束点为3号节点，从Us知道其对应最大位移出现在Sub2，则挑选出的位移列阵就会是Us的第二列，所选第一列和第二列都为假设，以此类推，将挑选出的各工况实位移按列存储到实位移列阵U。

步骤四，对刚度矩阵K进行约简，形成Ki矩阵。按照步骤三中不同类型的边界条件的个数，作为约简刚度矩阵K的循环次数，对刚度矩阵K进行约简，形成Ki矩阵。

对如何进行约简刚度矩阵K进行说明。将刚度矩阵K按边界条件共分为m组，只需对k个约束循环m次，即可剔除刚度矩阵K中重复的行和列，消除刚度矩阵K的奇异性，形成Ki矩阵。如不分组循环，将对整个刚度矩阵K进行筛选，运算量大，按相同边界条件循环便于计算。

步骤五，在矩阵Ki中，选取与步骤三中相同的位移约束点，在位移约束点上施加单位载荷，形成虚载荷列阵Pri。

单个位移约束是施加在某个节点的某一自由度上的，虚载荷就是在该节点的约束自由度上施加单位载荷，单位载荷为1，其余节点载荷为0，其中节点即位移约束点。

对如何进行虚载荷列阵Pri进行说明。矩阵Ki按边界条件共分为m组，假设某一组中有3个位移约束，分别施加在1号节点、2号节点、3号节点，则对应的虚位移载荷列阵分别为：

Pr1＝[1，0，…，0]^T

Pr2＝[0，1，…，0]^T

Pr3＝[0，0，1，…，0]^T

其中，每一个组内的循环次数取决于本组内位移约束个数。虚载荷列阵Pri中的i即为节点编号。

需要说明的是矩阵Ki的施加节点与上文步骤三中实位移列阵U中节点一一对应。

步骤六，根据平衡方程组Ki*Upsei＝Pri，分别得到不同位移约束及不同边界条件下的结构虚位移列阵Upsei，并将多个结构虚位移列阵Upsei构成虚位移矩阵Upse。

逐个计算每个位移约束对应虚载荷下得到的虚位移Upsei，并按列储存到Upse，最终Upse的列数为位移约束个数，顺序与数组U和S中位移约束的排序一致(U和S的含义见说明书具体实施方式步骤三中第二段)。

步骤七，计算变量的单元微分刚度矩阵Kde，并根据单元微分刚度矩阵Kde装配结构微分刚度矩阵Kd。

此处所述的变量与上文步骤一中提到的变量表示同一个事物，每一个变量对每一个位移约束都有一个敏度值，本步骤将刚度矩阵K关于当前变量求微分，每个变量都需要重复这个过程。此步骤刚度矩阵K关于当前变量求微分是为了带入位移敏度解析式求解位移敏度。

下面介绍位移敏度解析式的由来：

位移约束是工程结构优化中广泛使用的一种全局刚度约束形式，下文将从静力平衡方程出发，系统的推导位移敏度的解析表达式。

静力平衡方程一般表达式为：

KU＝F………………………………(1)

上式(1)两端对第i个变量x_i求导可得：

$\frac{\∂K}{\∂x_i}U+K\frac{\∂U}{\∂x_i}=\frac{\∂F}{\∂x_i}...\left(2\right)$

由于大多数静力模型中载荷并不随变量的变化发生变化，即认为载荷与变量不相关，故式(2)可简化为：

$\frac{\∂U}{\∂x_i}=-K^{-1}\frac{\∂K}{\∂x_i}U...\left(3\right)$

上式(3)即为位移敏度的解析式，但是实际工程问题中模型规模通常较大，结构刚度矩阵呈病态或准病态，难以直接求逆。本文引入虚载荷法来解决这一问题。过程如下：

设位移约束的自由度编号为r，则可生成一个虚载荷列阵，表示为式(4)：

将式(4)代入式(1)，变换后可以得到：

$U_{pse}^r=K^{-1}{P}_r...\left(5\right)$

式(3)做转置，可得：

$\frac{\∂U_r^T}{\∂x_i}=-U_r^T\frac{\∂K}{\∂x_1}{K}^{-1}...\left(6\right)$

式(6)两端同时乘虚载荷列阵得：

$\frac{\∂U_r}{\∂x_i}=-U_r^T\frac{\∂K}{\∂x_i}{K}^{-1}{P}_r...\left(7\right)$

将式(5)代入式(7)，第r个位移对第i个变量的位移导数可表述为：

$\frac{\∂U_r}{\∂x_i}=-U_r^T\frac{\∂K}{\∂x_i}{U}_{pse}^r...\left(8\right)$

上式(8)即为引入虚载荷法以后的位移敏度全解析式。

其中，观察式(3)和(8)，可发现两式都具有为结构微分刚度矩阵Kd，Kd为现有常规手段。

步骤八，依次提取实位移列阵U的各列并转置，与结构微分刚度矩阵Kd及虚位移列阵Upse的对应列相乘，得到当前变量的各位移约束敏度。

对于不同变量结构微分刚度矩阵Kd是不同的，依次取位移矩阵U和矩阵Upse的对应列以及当前变量的Kd，按照上文求得的式(8)计算当前变量对应的所有位移约束的敏度。

$\frac{\∂U_r}{\∂x_i}=-U_r^T\frac{\∂K}{\∂x_i}{U}_{pse}^r$

该公式一次计算一个变量关于一个位移约束的敏度，，为第i个变量中的第r个位移约束的敏度，为实位移矩阵U中第r列的转置，为结构微分刚度矩阵Kd，为虚位移矩列阵Upse的第r列。

将一个变量对应的各位移约束敏度依次计算完并存储，再对下一个变量进行计算，最后将所有敏度结果格式化输出。

下面举一个例子进行详细说明，假设飞机结构优化模型包含6个工况，6个变量，4个位移约束信息。

第一步中定义数组存储6个变量，4个位移约束信息。

第二步中得到实位移列阵Us，共6列，每列存储一个工况下的结构位移，假设第一列是第一种工况的实位移列阵U₁，节点位移分别为4、2、3、5，第二列是第二种工况的实位移列阵U ₂，节点位移分别为5、7、6、3，后三列以此类推，由五个位移列阵U组成了该实位移矩阵Us。

由于只有4个位移约束信息，则刚度矩阵K为4×4的方阵。

第三步该优化问题有4个位移约束，第1个位移约束点为3号节点，从Us知道其对应最大位移(即第3行中的数值8)出现在第4列(即Sub4)，则挑选出的位移列阵为Us的第4列，第2个位移约束点为2号节点，从Us知道其对应最大位移(即第2行中的数值6)也出现在第3列(即Sub3)，则挑选出的位移列阵为Us的第3列，第3个位移约束点为1号节点，从Us知道其对应最大位移(即第1行中的数值8)仍旧出现在第3列(即Sub3)，则挑选出的位移列阵为Us的第3列，第4个位移约束点为4号节点，从Us知道其对应最大位移(即第4行中的数值9)出现在第5列(即Sub5)，则挑选出的位移列阵为Us的第5列。即实位移列阵U有4列，分别储存Us的第4、3、3、5列。

假设将挑选出的4个工况按边界条件共分为2组，相同边界条件分一组，假设第1组包含约束1和4，第2组包含约束2和3。

用表格表述为：

U需要根据工况对应的边界号排列，先排列相同组别的不同约束条件，再排列不同组别，其中，按照选取的位移约束点的顺序排列，所选例子即按照约束号1、4、2、3的顺序来排列，得到

第四步将刚度矩阵K按边界条件共分为2组，只需对4个约束循环2次，即可消除刚度矩阵K的奇异性，剔除刚度矩阵K中重复的行和列，形成Ki矩阵。

第五步矩阵Ki按照边界条件分为2组，第一组中有2个位移约束，分别施加在1号节点、4号节点，则对应的虚位移载荷列阵分别为：Pr1＝[1，0，0，0]^T和Pr4＝[0，0，0，1]^T；第二组中有2个位移约束，分别施加在2号节点、3号节点，则对应的虚位移载荷列阵分别为：Pr2＝[0，1，0，0]^T和Pr3＝[0，0，1，0]^T，需要注意的是矩阵Ki的施加节点与上文第三步中实位移列阵U中节点一一对应。

第六步根据平衡方程Ki*Upsei＝Pri，可以计算出每一个虚载荷下的结构虚位移列阵Upsei，并组装成Upse；则Upse1为Upse的第一列，Upse4为Upse的第二列，Upse2为Upse的第三列，Upse3为Upse的第四列。

第七步，计算变量的单元微分刚度矩阵Kde，并根据单元微分刚度矩阵Kde装配结构微分刚度矩阵Kd，其中，变量为6个。

第八步依次提取实位移列阵U的各列并转置，与结构微分刚度矩阵Kd及虚位移列阵Upse的对应列相乘，得到当前变量的各位移约束敏度。

实位移列阵U、虚位移列阵Upse、结构微分刚度矩阵Kd都为4×4矩阵，可知一个变量有对应的4个位移约束敏度，将一个变量对应的4个位移约束敏度依次计算完并存储，再对剩余5个变量进行计算，最后将所有敏度结果格式化输出。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。