CN108763658A - 基于等几何方法的组合薄壁结构固有频率设计方法 - Google Patents

Abstract

一种基于等几何方法的组合薄壁结构固有频率设计方法，首先通过NURBS描述初始的薄壁组合结构各子结构的曲面形状，将子曲面归类为设计域和非设计域两类。设计域的形状通过控制点的移动来改变，将这些控制点的位置设定为设计变量。采用等几何基尔霍夫壳单元对结构进行网格划分，计算子结构的质量和刚度矩阵，进行模态分析。利用子结构的前若干阶特征值和特征向量，分别构建各子结构的缩减模型，通过固定界面模态综合法，形成组合结构的缩减模型，求解得到其固有频率。然后对结构进行设计灵敏度分析，获取固有频率对形状设计变量的解析灵敏度，最后通过基于梯度的优化算法求解设计问题。本发明简化了网格生成过程，同时提高了优化效率。

Description

基于等几何方法的组合薄壁结构固有频率设计方法

技术领域

本发明涉及薄壁结构优化设计领域，具体涉及一种基于等几何方法的组合薄壁结构固有频率设计方法，克服有限元薄壁组合结构优化设计中网格生成繁琐、优化效率低下的问题。

背景技术

薄壁组合结构在航空航天、汽车制造等行业中广泛使用，其固有频率与结构的形状密切相关。薄壁组合结构的固有频率反映了结构的动态特性，决定了其对动态载荷的响应。

航空航天和汽车工业中的薄壁组合结构往往需要通过形状调整，来改变其固有频率以达到相关目的。比如，通过改变汽车车盖的固有频率，使其振动噪声辐射最小。这类结构的形状设计一般采用基于有限单元的结构优化技术：首先对薄壁整体结构进行参数化，对参数赋值，驱动CAD模型生成，然后划分网格，通过有限元分析得到固有频率，再根据优化算法重复从赋值到分析这些环节，直至得到较优的固有频率。在这个流程中，有限元网格的自动化、高质量生成需要较长的时间，对于实际的结构分析，此部分所占用的时间大约是整个分析所耗费时间的80％，十分繁琐低效。设计灵敏度一般也采用差分来近似，即通过两次或者三次计算来获取一个设计变量在某点处的灵敏度，耗时较长，精度较低，拖慢了优化速度。另外实际中需要优化的往往是组合结构的某个部分，上述流程却对整体结构进行了网格划分，造成了分析模型较大，这也降低了分析效率，进而拖慢了优化设计的速度。

等几何分析方法是国外学者Hughes提出的一种融合CAD建模和CAE分析的方法(Hughes T J R,Cottrell J A,Bazilevs Y.Isogeometric analysis:CAD,finiteelements,NURBS,exact geometry and mesh refinement[J].Computer Methods inApplied Mechanics&Engineering,2005,194(39):4135-4195.)，通过采用CAD建模中常用的Non-uniform rational B-spline来描述模型上所依附的物理场，使得CAD模型与CAE模型能无缝衔接。在形状描述中NURBS的控制点位置用来控制模型的几何形状。在力学分析中，NURBS的控制点系数被用来表示节点位移，即整个位移场也是一个NURBS。该方法的网格划分是通过对CAD模型NURBS的节点插入和升阶得到的，该过程不改变几何形状，只是增加了控制点数目，且操作过程简单，网格划分几乎不占用时间。等几何基尔霍夫壳体基于等几何概念开发的薄壳单元，通过离散壳体的中面生成薄壳单元，每个节点仅有三个位移自由度来描述壳体的变形，所生成的刚度矩阵和质量矩阵的维度较少。采用这种单元对薄壁结构进行力学分析的效率高、精度好，对薄壁结构进行固有频率的优化设计能大幅减少网格划分的时间，得到更好的薄壁结构形式。

发明内容

本发明的目的在于针对上述现有技术中的问题，提供一种基于等几何方法的组合薄壁结构固有频率设计方法，简化传统有限元设计中的网格生成过程，同时提高优化效率。

为了实现上述目的，本发明基于等几何方法的组合薄壁结构固有频率设计方法包括：

步骤1，将初始的组合薄壁结构的各个子结构采用NURBS曲面进行形状描述，或者从CAD软件中直接导出曲面模型，给出模型的节点序列以及控制点位置；

步骤2，标示控制点位置能够变动的子结构为设计域，标示所要改变的控制点为设计变量，标示控制点位置不能够变动的子结构为非设计域；

步骤3，形状能够改变的子结构为设计域子结构，反之为非设计域子结构，将非设计域子结构采用等几何基尔霍夫壳单元进行网格划分，计算其刚度矩阵和质量矩阵；

步骤4，根据Craig-Bampton模态综合法，构建非设计域子结构的动力学缩减模型，给出缩减后的非设计域结构刚度矩阵和质量矩阵；

步骤5，将设计域子结构采用等几何基尔霍夫壳单元进行网格划分，计算其刚度矩阵和质量矩阵，解析计算单元刚度矩阵和质量矩阵对单元节点位置变化的导数；

步骤6，结合步骤4缩减后的非设计域结构刚度矩阵和质量矩阵，通过Craig-Bampton模态综合法构建组合结构整体的刚度矩阵和质量矩阵，计算整体结构的固有频率；

步骤7，采用解析及灵敏度映射方法计算整体结构对于设计变量的设计灵敏度；

步骤8，使用基于梯度的优化算法，针对结构的固有频率设计要求，进行形状优化设计。

所述的步骤7具体包括以下步骤：7.1、计算NURBS曲面进行形状描述后的初始曲面模型与网格划分后的模型之间的几何映射关系，构建初始模型控制点位置与细化模型控制点位置之间的关系，初始模型控制点位置与细化模型控制点位置之间为线性关系，能够通过NURBS节点插入或者升阶算法得到；7.2、根据步骤5单元刚度矩阵和质量矩阵对单元节点位置变化的导数，计算出整体结构对于细化模型控制点位置的设计灵敏度；7.3、根据初始模型控制点位置与细化模型控制点位置的线性关系，以及整体结构对于细化模型控制点的设计灵敏度，通过灵敏度映射方法获得整体结构对于设计变量的设计灵敏度。

所述的步骤8中若优化迭代收敛时，得到所需要的组合薄壁结构形状。

采用NURBS曲面进行形状描述时通过下式表征：

其中，ξ和η为两个参数，其定义域为[0,1]；R为NURBS基函数，p、q分别表示基函数在两个参数方向上的阶数；P_i,j为第(i，j)个NURBS控制点的坐标，包括x,y,z三个分坐标；

该曲面共有n×m个控制点；

N为B-spline基函数，其定义在节点矢量上，两个方向上节点矢量分别为：Ξ＝{ξ₁，ξ₂，…，ξ_n+p+1}和H＝{η₁，η₂，…，η_n+p+1}，p和q分别为基函数的阶数，w_i,j为对应于(i，j)控制点的权值；

结构应变能小于设定值，总重量小于设定值，第一阶固有频率最大化，问题描述如下：

find z＝[z₁，z₂，…，z₂5]

maximize λ₁(z)

s.t.m(z)≤m_d

W＝u^TK(z)u≤W_d

K(z)u＝F

其中，z表示设计变量，其控制结构的形状；M为结构的质量矩阵，K为结构的刚度矩阵，F为结构外载荷，u为结构的变形场，λ表示结构的固有频率(特征值)；m为结构的总质量，m_d表示许用的最大质量，W表示结构总的应变能，W_d为许可的最大应变能。

步骤3计算刚度矩阵和质量刚度矩阵时的系统方程为：

Mü＋Ku＝0. (4)

刚度矩阵为：

K_lm表示刚度矩阵的第l行m列元素，δΠ_int为结构的内力虚功，u_l表示节点l的位移；

质量矩阵为：

M_ij＝∫_sρN_iN_jd_s

(6)

其中，ρ为结构的密度，N_i表示基函数；

n和m分别表示壳中面等效集中力和转矩，u表示壳体中面位移，ε和κ分别表示壳体的中面应变和曲率应变；S表示壳体的中面参数化曲面，()_,α表示括号中变量对标号为a的参数求导数，()_,αβ表示括号中变量对α和β参数依次导数，α,β∈{1,2}，N₃表示曲面的单位外法线；C^αβγδ表示结构的弹性张量分量，E表示材料的弹性模量，ν表示泊松比。

步骤4构建缩减后的非设计域结构刚度矩阵和质量矩阵时：

首先对子结构进行模态分析，然后求解出来各个子结构的内节点固定界面模态和边界节点模态 其是子结构前k个固定界面特征向量的集合，k远远小于该子结构总的自由度数目，表示子结构第r个特征向量，其通过对该子结构所有边界节点全约束后进行模态分析得到；

其中，K^s表示第s个子结构的刚度矩阵，ii和ib分别指代其内节点和边界节点对应的刚度矩阵分块，I代表单位刚度矩阵。

步骤5通过下式解析计算单元刚度矩阵和质量刚度矩阵对单元节点位置变化的导数：

其中，M_ij为结构质量矩阵的第i行j列元素，P_k为壳中面S的第k个控制点坐标，K_lm为结构刚度矩阵的第l行m列元素，u_l为位移矢量的第l个元素，N为形函数；()_,α表示括号中变量对标号为a的参数求导数。

步骤6模态变换矩阵为模态综合后的系统方程为：

其中，T为结构的变换矩阵，Φ_C和Φ_N分别表示结构的边界节点模态矢量集合和固定界面内节点模态矢量集合，w为结构的固有圆角频率；因为固定界面内节点模态集为选的前k阶模态向量，k远远小于结构自由度数目，通过上述矩阵操作达到了方程维度缩减的目的；求解上述方程即得到简化系统的固有频率。

根据单元刚度矩阵和质量刚度矩阵对单元节点位置变化的导数，计算整体结构对于网格节点的设计灵敏度；

固有频率(特征值)的细化模型设计灵敏度为：

质量的细化模型设计灵敏度为：

应变能的细化模型设计灵敏度为：

W，＝-u^TK，u. (17)

上述式子中，sum(M)表示对结构质量矩阵M的所有行列元素进行加和，()_,α表示括号中的量对下标中的变量进行求导，P_k表示第k个控制点坐标，K表示结构的刚度矩阵，u表示结构的变形矢量。

步骤8使用MATLAB中的fmincon函数，针对结构的固有频率设计要求，进行形状优化设计。

与现有技术相比，本发明具有如下的有益效果：采用等几何基尔霍夫壳单元对薄壁组合结构进行力学分析，减少了网格划分的负担，从CAD模型生成分析所用网格模型所需的时间较少，步骤也简单。网格模型与CAD模型间也不存在离散误差，提高了形状描述的精度，也提高了分析的精度。在优化设计中，用此方法替换传统的有限元方法使得花在网格划分环节上的时间减少，分析效率提高。本发明采用模态综合法来获取薄壁组合结构的固有频率，省掉了对形状不发生变化的子结构进行重复分析的负担。采用模态综合法后，薄壁组合结构中只有形状发生变化的子结构的结构矩阵需要多次生成，形状不变的子结构模型通过模态综合法被简化，并且在优化中不发生变化，系统的固有频率通过综合可变的子结构的矩阵和不变结构的缩减模型来得到，实现了仅对组合结构一部分进行重分析，就得到系统固有频率的效果，大大提高了力学分析的效率。本发明采用基于灵敏度映射的解析灵敏度计算策略，灵敏度计算快速而精确，加速了设计的收敛，能够对存在较多设计变量的组合薄壁结构进行设计。通过首先计算网格层面的灵敏度，然后构建设计变量和网格控制点位置间的线性映射矩阵，最后映射得到设计参量的设计灵敏度。这种策略具有较大的通用型，编程方便，可处理较多类型的设计参量，能在更大范围内寻求结构的最优形状，优化过程也更快速。

附图说明

图1：本发明基于等几何法进行组合薄壁结构固有频率设计的流程图；

图2：本发明实施例的初始的组合薄壁曲面及其分析网格模型，其由两个子曲面构成，Ω1为设计域，Ω2为非设计域；

图3：本发明实施例的初始的设计域曲面的轮廓曲线的控制点z向坐标及其优化后的值，以及优化前后的结构固有频率、应变能和质量；

图4:本发明实施例的薄壁结构经过优化最终得到的最优薄壁组合结构的形状；

图5：本发明实施例薄壁组合结构第一阶固有频率设计迭代过程图。

具体实施方式

下面结合附图对本发明做进一步的详细说明。

参见图1，本发明基于等几何方法的组合薄壁结构固有频率设计方法，包括以下步骤：

步骤(1)将初始的组合薄壁结构的各个子结构采用NURBS曲面进行形状描述，该曲面模也可从Rhino软件或者其他CAD软件中导出。任意子曲面的几何描述为：

其中，R为NURBS基函数，P_i,j为NURBS控制点的坐标，其包括x,y,z三个分坐标，该曲面共有n×m个控制点，ξ和η为两个参数其定义域为[0,1]。

其中，N为B-spline基函数，其定义域节点矢量上，在两个方向上节点矢量分别为：Ξ＝{ξ₁，ξ₂，…，ξ_n+p+1}和H＝{η₁，η₂，…，η_n+p+1}。

p和q分别为基函数的阶数，w_i,j为对应与(i，j)控制点的权值。

参见图2，其中Ω1的轮廓曲线的控制点位置作为设计变量，要求结构在应变能小于一定值，总重量小于一定值，其第一阶固有频率最大化。设计问题描述如下：

find z＝[z₁，z₂，…，z₂5]

maximize λ₁(z)

s.t.m(z)≤m_d

W＝u^TK(z)u≤W_d

K(z)u＝F

步骤(2)标示控制点位置可以变动的子结构Ω₁为设计域，标示控制点位置不可以变动的子结构为非设计域Ω₂。选定某些子曲面控制点的位置作为设计变量，通过优化其位置来优化结构的固有频率。在实施例中选择Ω₁的轮廓曲线的控制点z向位置作为设计变量，见附图2，z的初始取值见附图3。

步骤(3)将非设计域子结构Ω₂采用等几何基尔霍夫壳单元进行网格划分，计算其刚度矩阵和质量刚度矩阵。其系统方程为：

Mü+Ku＝0. (4)

刚度矩阵为：

质量矩阵为：

M_ij＝∫_sρN_iN_jds

(6)

其中：

步骤(4)根据Craig-Bampton模态综合法，构建非设计域子结构Ω₂的动力学缩减模型，给出缩减后的非设计域结构的刚度矩阵和质量矩阵。首先对子结构进行模态分析，然后求解出来各个子结构的内节点固定界面模态和边界节点模态

其中，ii和ib分别指代内节点和边界节点的刚度矩阵块。

步骤(5)将设计域子结构采用等几何基尔霍夫壳单元进行网格划分，计算其刚度矩阵和质量刚度矩阵，解析计算单元刚度矩阵和质量矩阵对单元节点位置变化的导数。其质量/刚度矩阵计算方法同步骤(3)相同，其导数的解析计算如下：

步骤(6)结合步骤(4)给出的非设计域结构缩减的单元刚度矩阵和质量矩阵，通过Craig-Bampton模态综合法构建组合结构整体的单元刚度矩阵和质量矩阵。计算整体结构的固有频率。模态变换矩阵为模态综合后的系统方程为：

求解即可得到简化系统的固有频率。

步骤(7)计算整体结构对于设计变量的设计灵敏度。整体结构对于设计变量的设计灵敏度是解析计算的，且是通过灵敏度映射方法计算的。其包含以下步骤：

7.1计算初始NURBS模型与网格细化后的分析模型之间的几何映射关系，即构建初始模型控制点位置与细化模型控制点位置之间的关系，其可以表述为矩阵，通过NURBS节点插入或者升阶算法得到，将其标记为T_a。

7.2根据步骤(5)获得的细化模型的质量、刚度矩阵的导数，计算整体结构对于网格节点的设计灵敏度。

固有频率的细化模型设计灵敏度为：

质量的细化模型设计灵敏度为：

应变能的细化模型设计灵敏度为：

W，＝-u^TK，u. (17)

7.3根据步骤(1)的线性关系和步骤7.2的灵敏度，通过灵敏度加权(即灵敏度映射)获得初始模型的设计灵敏度。固有频率的设计灵敏度为：

质量和应变能的设计灵敏度也同样采用上述式子进行灵敏度映射得到。

步骤(8)基于步骤(7)，使用MATLAB中的fmincon函数，针对结构的固有频率设计要求，进行形状优化设计。设计变量在步骤(2)中已经标示。当优化进行迭代时，只有设计域子结构的形状会变化，需要重复步骤(6)-(8)。非设计域的缩减单元刚度矩阵和质量矩阵在迭代中没有变化，重复使用。当优化迭代收敛时候，即可得到设计后的组合薄壁结构形状,具体实施例的优化形状如图4所示。图2展示了结构优化前后的形状可变的子曲面Ω₁的轮廓控制点的z坐标，薄壁组合结构的第一阶固有频率得到了提升，结构的应变能和质量也没有超过设定值，实现了薄壁组合结构在保证结构刚度、质量的情况下提高其基频的效果。

Claims (10)

1.一种基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于，包括以下步骤：

步骤1，将初始的组合薄壁结构的各个子结构采用NURBS曲面进行形状描述，或者从CAD软件中直接导出曲面模型，给出模型的节点序列以及控制点位置；

步骤2，标示控制点位置能够变动的子结构为设计域，标示所要改变的控制点为设计变量，标示控制点位置不能够变动的子结构为非设计域；

步骤3，形状能够改变的子结构为设计域子结构，反之为非设计域子结构，将非设计域子结构采用等几何基尔霍夫壳单元进行网格划分，计算其刚度矩阵和质量矩阵；

步骤4，根据Craig-Bampton模态综合法，构建非设计域子结构的动力学缩减模型，给出缩减后的非设计域结构刚度矩阵和质量矩阵；

步骤5，将设计域子结构采用等几何基尔霍夫壳单元进行网格划分，计算其刚度矩阵和质量矩阵，解析计算单元刚度矩阵和质量矩阵对单元节点位置变化的导数；

步骤6，结合步骤4缩减后的非设计域结构刚度矩阵和质量矩阵，通过Craig-Bampton模态综合法构建组合结构整体的刚度矩阵和质量矩阵，计算整体结构的固有频率；

步骤7，采用解析及灵敏度映射方法计算整体结构对于设计变量的设计灵敏度；

步骤8，使用基于梯度的优化算法，针对结构的固有频率设计要求，进行形状优化设计。

2.根据权利要求1所述基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于，所述的步骤7具体包括以下步骤：7.1、计算NURBS曲面进行形状描述后的初始曲面模型与网格划分后的模型之间的几何映射关系，构建初始模型控制点位置与细化模型控制点位置之间的关系，初始模型控制点位置与细化模型控制点位置之间为线性关系，能够通过NURBS节点插入或者升阶算法得到；7.2、根据步骤5单元刚度矩阵和质量矩阵对单元节点位置变化的导数，计算出整体结构对于细化模型控制点位置的设计灵敏度；7.3、根据初始模型控制点位置与细化模型控制点位置的线性关系，以及整体结构对于细化模型控制点的设计灵敏度，通过灵敏度映射方法获得整体结构对于设计变量的设计灵敏度。

3.根据权利要求1所述基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于：所述的步骤8中若优化迭代收敛时，得到所需要的组合薄壁结构形状。

4.根据权利要求1所述基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于：采用NURBS曲面进行形状描述时通过下式表征：

其中，ξ和η为两个参数，其定义域为[0，1]；R为NURBS基函数，p、q分别表示基函数在两个参数方向上的阶数；P_i，j为第(i，j)个NURBS控制点的坐标，包括x，y，z三个分坐标；该曲面共有n×m个控制点；

N为B-spline基函数，其定义在节点矢量上，两个方向上节点矢量分别为：Ξ＝{ξ₁，ξ₂，…，ξ_n+p+1}和H＝{η₁，η₂，…，η_n+p+1}，p和q分别为基函数的阶数，w_i，j为对应于(i，j)控制点的权值；

结构应变能小于设定值，总重量小于设定值，第一阶固有频率最大化，问题描述如下：

find z＝[z₁，z₂…，z₂5]

maximize λ₁(z)

s.t.m(z)≤m_d

W＝U^TK(z)u≤W_d

K(z)u＝F

其中，z表示设计变量，其控制结构的形状；M为结构的质量矩阵，K为结构的刚度矩阵，F为结构外载荷，u为结构的变形场，λ表示结构的固有频率(特征值)；m为结构的总质量，m_d表示许用的最大质量，W表示结构总的应变能，W_d为许可的最大应变能。

5.根据权利要求1所述基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于：步骤3计算刚度矩阵和质量刚度矩阵时的系统方程为：

Mü＋Ku＝0. (4)

刚度矩阵为：

K_lm表示刚度矩阵的第l行m列元素，δΠ_int为结构的内力虚功，u_l表示节点l的位移；

质量矩阵为：

M_ij＝∫_sρN_iN_jds

(6)

其中，ρ为结构的密度，N_i表示基函数；

n和m分别表示壳中面等效集中力和转矩，u表示壳体中面位移，ε和κ分别表示壳体的中面应变和曲率应变；S表示壳体的中面参数化曲面，()_,α表示括号中变量对标号为a的参数求导数，()_,αβ表示括号中变量对α和β参数依次导数，α,β∈{1,2}，N₃表示曲面的单位外法线；C^αβγδ表示结构的弹性张量分量，E表示材料的弹性模量，ν表示泊松比。

6.根据权利要求1所述基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于，步骤4构建缩减后的非设计域结构刚度矩阵和质量矩阵时：首先对子结构进行模态分析，然后求解出来各个子结构的内节点固定界面模态和边界节点模态其是子结构前k个固定界面特征向量的集合，k远远小于该子结构总的自由度数目，表示子结构第r个特征向量，其通过对该子结构所有边界节点全约束后进行模态分析得到；

其中，K^s表示第s个子结构的刚度矩阵，ii和ib分别指代其内节点和边界节点对应的刚度矩阵分块，I代表单位刚度矩阵。

7.根据权利要求1所述基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于，步骤5通过下式解析计算单元刚度矩阵和质量刚度矩阵对单元节点位置变化的导数：

其中，M_ij为结构质量矩阵的第i行j列元素，P_k为壳中面S的第k个控制点坐标，K_lm为结构刚度矩阵的第l行m列元素，u_l为位移矢量的第l个元素，N为形函数；()_,α表示括号中变量对标号为a的参数求导数。

8.根据权利要求1所述基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于，步骤6模态变换矩阵为模态综合后的系统方程为：

其中，T为结构的变换矩阵，Φ_C和Φ_N分别表示结构的边界节点模态矢量集合和固定界面内节点模态矢量集合，w为结构的固有圆角频率；因为固定界面内节点模态集为选的前k阶模态向量，k远远小于结构自由度数目，通过上述矩阵操作达到了方程维度缩减的目的；求解上述方程即得到简化系统的固有频率。

9.根据权利要求1所述基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于，根据单元刚度矩阵和质量刚度矩阵对单元节点位置变化的导数，计算整体结构对于网格节点的设计灵敏度；

固有频率(特征值)的细化模型设计灵敏度为：

质量的细化模型设计灵敏度为：

应变能的细化模型设计灵敏度为：

W，＝-u^TK，u.

(17)

上述式子中，sum(M)表示对结构质量矩阵M的所有行列元素进行加和，()_,α表示括号中的量对下标中的变量进行求导，P_k表示第k个控制点坐标，K表示结构的刚度矩阵，u表示结构的变形矢量。

10.根据权利要求1所述基于等几何方法的组合薄壁结构固有频率设计方法，其特征在于，步骤8使用MATLAB中的fmincon函数，针对结构的固有频率设计要求，进行形状优化设计。