一种功能梯度板力学特性数值仿真方法

技术领域

本发明涉及功能梯度板力学特性数值仿真技术领域，尤其涉及一种功能梯度板力学特性数值仿真方法。

背景技术

目前，功能梯度材料是两种或多种材料复合且成分和结构呈连续梯度变化的一种新型复合材料，基于其性质、组成或结构在空间上逐渐变化，进而可实现特定的功能或性能，以及其具有耐高温、承载能力大、疲劳寿命长等力学特性，在航空航天、医疗器械、电子和光电子器件等许多领域都具有潜在的用途。

其中，现有技术中为了考虑功能梯度板材料性质沿厚度方向非均匀性与传统壳单元失效，而需要将板结构在整个三维计算域内进行离散，包括厚度方向。然而这类建模方法存在较大的局限性，例如，目前发展较为成熟的有限元法，为了适应板结构面内尺寸远大于厚度方向尺寸的实际情况，避免出现数值问题(如自锁、网格依赖等)，满足工程精度要求，需要将计算域离散为较细的网格，这很大程度上耗费了大量的计算自由度，极大地提高了计算成本。以及专利公开号为CN114492123A的发明专利公开了一种基于分层法对功能梯度材料结构进行拓扑优化方法，虽将功能梯度材料结构划分成若干层，每层的材料参数按函数形式变化，相邻层之间的材料参数有较好的连续性，每层按照常规的有限元方法进行网格划分，建立有限元模型求解功能梯度材料的力学问题，基于不同的需求对功能梯度材料结构进行拓扑优化设计。但实际的计算量还是偏大，计算效率不高。以进而，亟需提出一种较低计算成本的功能梯度板的力学特性数值仿真方法。

发明内容

本发明的主要目的在于提供一种功能梯度板力学特性数值仿真方法，旨在解决现有获取功能梯度板力学特性仿真方法的计算成本高的技术问题。

为实现上述目的，本发明提供一种功能梯度板力学特性数值仿真方法，所述方法包括以下步骤：

S1、获取目标功能梯度板对应的结构模型，并获取结构模型中的多个第一控制点的特征信息，各第一控制点的特征信息包括各第一控制点的坐标及其权重与节点向量；

S2、采用升阶或节点加密算法将结构模型依据各第一控制点进行单元细分，并将各细分单元进行插值近似，得到细化区域内的多个第二控制点的特征信息、各第一控制点与对应的第二特征点的索引信息；

S3、将步骤S2中的各第二控制点的特征信息基于三维弹性理论，并采用目标功能梯度板对应的弹性矩阵，创建目标功能梯度板的静力学等几何比例边界有限元系统方程；

S4、根据目标功能梯度板的位移与节点力将静力学等几何比例边界有限元系统方程转化为一阶常微分方程，并采用Padé级数进行矩阵指数运算，得到基于各第二控制点的位移和节点力的矩阵；

S5、根据基于各第二控制点的位移和节点力的矩阵与各第二控制点的坐标、目标功能梯度板的厚度得到目标功能梯度板内各第二控制点的场变量以实现根据目标功能梯度板内各第二控制点的场变量确定目标功能梯度板内部各位置的力学特性。

可选地，采用基函数R(η,ζ)对各细分单元的位移u场进行插值近似，得到位移u场在xyz方向的基函数表示：

其中i为第二控制点编号，η、ζ为第二控制点的参数坐标，Ne＝(p+1)×(q+1)为一个单元内第二控制点数量，p、q分别为基函数沿η、ζ方向的阶次。

可选地，所述步骤S3具体包括：

根据三维弹性理论以及采用胡克定律得到应力-应变的矩阵；

获取应变-位移的矩阵，根据应力-应变之间的矩阵和应变-位移的矩阵得到基于NURBS的应力-位移的矩阵；

根据应力-应变之间的矩阵沿着细分单元进行积分得到内节点力{q}矩阵；

根据基于NURBS的应力-位移的矩阵并采用虚功原理进行整理得到基于位移的方程组；

根据内节点力{q}矩阵和基于位移的方程组得到位移、内力关于z坐标的一阶微分方程，即为目标功能梯度板的静力学等几何比例边界有限元系统方程。

可选地，所述应力与应变的矩阵为：

{σ}＝{σ

其中τ表示切应力，{σ}表示应力与{ε}表示应变，[D(z)]表示弹性矩阵

式中

可选地，所述获取应变与位移的矩阵的步骤包括：

根据位移u场在xyz方向的基函数表示，可以得到单元内任一点坐标可表示为

并根据单元内任一点坐标获取应变-位移的矩阵；

所述应变-位移的矩阵为：

{ε}＝[B

[B

式中|J|＝x

可选地，所述基于NURBS的应力-位移的矩阵为：

可选地，所述步骤S3中的目标功能梯度板的静力学等几何比例边界有限元系统方程为：

可选地，在步骤S4中，所述根据目标功能梯度板的位移与节点力将静力学等几何比例边界有限元系统方程转化为一阶常微分方程的步骤，具体包括：

设定

可选地，在步骤S4中，所述并采用Padé级数进行矩阵指数运算，得到基于各第二控制点的位移和节点力的矩阵的步骤，具体包括：

根据z

[Q

可选地，在所述步骤S3之前，需要确定目标功能梯度板是静力分析或模态分析；

若目标功能梯度板是静力分析，则执行步骤S3；

若目标功能梯度板是模态分析，则将步骤S2中的各第二控制点的特征信息基于三维弹性理论，并采用目标功能梯度板对应的弹性矩阵和质量矩阵，创建目标功能梯度板的动力学等几何比例边界有限元系统方程。

有益效果：

(1)本发明采用单元细分并结合NURBS基函数对功能梯度板的位移场进行插值近似，进而实现将功能梯度板的三维问题转换成二维平面的计算，实现将功能梯度板仿真具体化呈现；

(2)本发明基于三维弹性理论、弹性矩阵、虚功理论和Padé级数将功能梯度板的静力分析转换成求解矩阵方程，进而降低计算量且有效提高计算效率和计算精度，最终可高效获取功能梯度板内部的力学特征；以及在模态分析下，引用质量矩阵，可实现对功能梯度板内部的力学特征的快速确定。

(3)此外，还可实现各种复杂形状的功能梯度板内部的力学特征的确定，并验证了本方法具有灵活适用性。

附图说明

图1为本发明功能梯度板力学特性数值仿真方法第一实施例的流程示意图；

图2为图1所示的示意简图；

图3为承受分布荷载正方形功能梯度板几何形状及尺寸简图；

图4为图3所示的正方形功能梯度板网格离散示意图；

图5为幂律指数为n＝1时不同网格剖分情况下正方形Al/Al

图6为幂律指数为n＝10时不同网格剖分情况下正方形Al/Al

图7为受均布荷载圆形Ti/ZrO

图8为圆形Ti/ZrO

图9为不同幂律指数情况下简支圆形Ti/ZrO

图10为不同幂律指数情况下固支圆形Ti/ZrO

图11为复杂正方形开孔板几何形状及尺寸示意图；

图12是为复杂正方形开孔板网格离散示意图；

图13是n＝5，h/L＝0.05时复杂正方形开孔板前6阶模态图。

具体实施方式

应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。

如图1-2所示，图1为本发明的一种功能梯度板力学特性数值仿真方法提供第一实施例的流程示意图，图2为图所示的示意简图，在图2中包括结构模型的建立及离散、系统方程的构造、方程的求解及后处理过程，各过程涉及目标功能梯度板实体1、控制网格2、控制点3、等几何比例边界有限单元4、母单元5、参数单元6，所述模型的建立及离散过程的操作对象为实体功能梯度板1。以及详细的方法具体包括以下步骤：

步骤S1、获取目标功能梯度板对应的结构模型，并获取结构模型中的多个第一控制点的特征信息，各第一控制点的特征信息包括第一控制点的坐标及其权重与节点向量。

具体地，如图2所示，采用Rhinoceros软件建立结构模型并保存为“fbx”格式，然后获取读取文件中原始数据，并根据原始数据将目标目标功能梯度板进行离散化并获取多个第一控制点及其特征信息，例如，设置xy平面在功能梯度板底面，z轴垂直于xy平面向上，板厚为h，进而各第一控制点的特征信息包括各第一控制点的坐标xy及其权重w与节点向量knotvector。

步骤S2、采用升阶或节点加密算法将结构模型依据各第一控制点进行单元细分，并将各细分单元进行插值近似，得到细分区域的多个第二控制点的特征信息。具体地，根据各第一控制点的分布情况将结构模型细分为多个更小单元区域，并采用NURBS基函数R(η,ζ)对各细分区域的位移u场进行插值近似，得到位移u场在xyz方向的基函数表示：

其中i为第二控制点编号，η、ζ为第二控制点的参数坐标，Ne＝(p+1)×(q+1)为一个单元内第二控制点数量，p、q分别为基函数沿η、ζ方向的阶次。

进一步地，还可根据位移u场在xyz方向的基函数表示，可以得到单元内任一点坐标可表示为

S3、将步骤S2中的各第二控制点的特征信息基于三维弹性理论，并采用目标功能梯度板对应的弹性矩阵，创建目标功能梯度板的静力学等几何比例边界有限元系统方程。其中，所述步骤S3具体包括以下内容：

步骤S31，根据三维弹性理论以及采用胡克定律得到应力-应变的矩阵，得到的应力-应变的矩阵为：

{σ}＝{σ

其中τ表示切应力，{σ}表示应力与{ε}表示应变，[D(z)]表示弹性矩阵。并且，功能梯度板的杨氏模量E、泊松比v和密度ρ可以表示为两种不同材料沿z坐标方向按照一定规则组合的函数形式：

E(z)＝E

式中

步骤S32，获取应变-位移的矩阵，根据应力-应变的矩阵和应变-位移的矩阵得到基于NURBS的应力-位移的矩阵。

具体地，根据单元内任一点坐标获取应变与位移的矩阵；所述应变-位移的矩阵为：

{ε}＝[B

其中[B

以及式中|J|＝x

然后将式(3)应变-位移的矩阵带入式(1)应变-位移的矩阵，得到基于NURBS的应力-位移的矩阵，其矩阵为：

步骤S33，根据应力-应变之间的矩阵沿着细分单元进行积分得到内节点力{q}矩阵，即

步骤S34，根据基于NURBS的应力-位移的矩阵并采用虚功原理进行整理得到基于位移的方程组。

具体地，基于虚功原理，即{q(z

(6)，式中，

上述式子中{q}为内部节点力，{F}为等效节点荷载，下标“t”、“b”分别表示变量对应于板顶、板底，

步骤S35，根据内节点力{q}矩阵和基于位移的方程组得到位移、内力关于z坐标的一阶微分方程，即为目标功能梯度板的静力学等几何比例边界有限元系统方程。

具体地，根据式(5)、式(6)得到位移、内力关于z坐标的一阶微分方程，即功能梯度板等几何比例边界有限元系统方程为

其中，在步骤S3之前，还包括确定目标功能梯度板是静力分析或模态分析；

若目标功能梯度板是静力分析，则执行步骤S3；

若目标功能梯度板是模态分析，则将步骤S2中的各第二控制点的特征信息基于三维弹性理论，并采用目标功能梯度板对应的弹性矩阵和质量矩阵，创建目标功能梯度板的动力学等几何比例边界有限元系统方程。

S4、根据目标功能梯度板的位移与节点力将静力学等几何比例边界有限元系统方程转化为一阶常微分方程，并采用Padé级数进行矩阵指数运算，得到基于各第二控制点的位移和节点力的矩阵。

具体地，设定

其中，

并得到一阶常微分方程的解可表示为

其中{c}为积分常数。基于z

其中

因此，若是对目标功能梯度板进行静力分析，则在将功能梯度板静力问题转化为求解矩阵方程

而对于模态分析，通过加速度引起的动能公式得到

由此得质量矩阵

最终功能梯度板模态分析广义特征方程可表示为([K]-ω

进一步地，将求得的方程解进行后处理，具体包括：

S5、根据基于各第二控制点的位移和节点力的矩阵与各第二控制点的坐标、目标功能梯度板的厚度得到目标功能梯度板内各第二控制点的场变量，进而可实现根据目标功能梯度板内各第二控制点的场变量确定目标功能梯度板内部各位置的力学特性。具体地，将求得的第二控制点位移、节点力、第二控制点的坐标以及目标功能梯度板的厚度信息带入至方程解终得到功能梯度板内任意点场变量，并将结果数据导入origin或者Tecplot软件，进行后处理，即画出位移、应力曲线或云图，进而可实现直观的确定功能梯度板内任意位置的力学特性。

本实施例中，通过获取目标功能梯度板对应的结构模型，采用升阶或节点加密算法将结构模型依进行单元细分，并将各细分单元进行插值近似，并基于三维弹性理论，采用目标功能梯度板对应的弹性矩阵，创建目标功能梯度板的静力学等几何比例边界有限元系统方程；通过求解方程得到基于位移和节点力的矩阵，进而根据基于位移和节点力的矩阵与各控制点的坐标、目标功能梯度板的厚度得到目标功能梯度板内各控制点的场变量；最后根据目标功能梯度板内各控制点的场变量确定目标功能梯度板内部各位置的力学特性。进而降低目标梯度板内部的力学特性的计算成本，进而提高计算效率。

进一步地，为了验证本发明在功能梯度板力学特性数值仿真方面的适用性、精确性和高效性，将给出从不同角度考察本发明的计算性能。

如图3所示，给定一个边长为L的正方形板，其弹性模量、泊松比分别为E和v，厚度为h，具体参数取值根据考察对象的不同分别确定。进而得到正方形梯度板的离散图，如图4所示，为了后续考察本发明提出方法的网格鲁棒性，给出了本模型的三种可能的网格剖分形式，分别定义为“网格1”、“网格2”、“网格3”。

表1-不同厚度正方形功能梯度板中心挠度对比表

如表1所示，将承受均布力q

进一步地，通过分析幂律指数为n＝1、10时不同网格剖分情况下正方形Al/Al

其中，如图5所示，其中幂律指数为n＝1时Al/Al

以及如图6所示，分析了幂律指数为n＝10时Al/Al

进一步地，通过分析在受均布荷载下的功能梯度板的情况，以Ti/ZrO

进一步，通过采用本方法确定幂律指数为n＝0、2、4、10情况下简支Ti/ZrO

进一步地，分析了幂律指数为n＝0、2、4、10情况下固支Ti/ZrO

进一步地，获取不同幂律指数情况下简支圆形Ti/ZrO

表2是不同幂律指数情况下简支圆形Ti/ZrO

由上表2可知，计算了厚度比为r/h＝10、20情况下简支圆形功能梯度板中心无量纲挠度，并与基于简化4变量剪切和法向变形理论等几何分析(sSNDT IGA)、简化一阶剪切变形理论(sFSDT)以及解析解进行了对比，可以看出在不同厚度比、幂律指数情况下本发明计算结果精度均高于传统数值方法。

进一步地，为了考察本发明所述方法对复杂几何形状的适用性，给出了一个四周固支的复杂开孔板，如图11所示，所述板为Ti/ZrO

需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。

上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。