CN113051786A - 基于绝对节点坐标法的压电层合薄板大变形的仿真方法 - Google Patents

Abstract

本发明公开了一种基于绝对节点坐标法的压电层合薄板大变形的仿真方法，包括，建立了压电层合悬臂薄板模型，该模型中间层为线弹性材料，下层为压电传感器，上层为压电作动器；在输入驱动电压，外力大小，材料参数和模型参数后建立非线性动力学方程，通过数值算法求解该动力学方程。本发明能够使用绝对节点坐标法计算出压电层合薄板在施加电压驱动后产生的变形，并输出相应的位形图，可视化压电驱动的效果；同时可以计算出压电层合薄板在收到外力作用下产生的变形以及传感器层产生的电压，并输出变形和电压随时间变化的关系。

Description

基于绝对节点坐标法的压电层合薄板大变形的仿真方法

技术领域

本发明属于柔性多体系统动力学领域，具体为一种基于绝对节点坐标法的压电层合薄板大变形的仿真方法。

背景技术

MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。基于 MATLAB对柔性多体系统进行动力学建模并计算和分析系统的动力学响应被学者们广泛地应用。

绝对节点坐标法是一种适合柔性系统大变形建模仿真的方法。在过去的20 年中，大量学者尝试了不同的新单元来描述更为复杂的机械模型或解决早期研究中出现的诸如泊松闭锁和剪切闭锁等问题，并提高计算的精度。绝对节点坐标法具有常质量矩阵，零离心力和科氏力等特点，并且可以精确描述刚体的运动。唯一形式复杂的部分是刚度阵，它是由连续介质力学或弹性线方法推导出的。

绝对节点坐标法法也可用于小变形或大变形的压电层合结构。Nada在《Absolutenodal coordinate formulation oflarge-deformation piezoelectric laminatedplates》中基于绝对节点坐标法建立了多层压电薄板单元模型，并分析了由压电作动器驱动的三层薄板结构的振动特性。其计算的数值结果与商业有限元软件 COMSOL的仿真结果进行对比结论较为一致。但是Nada在推导过程的弹性力项中忽略了拉伸项应变，这使得该单元只能用于低张力板；同时，其建立的层合板离散模型使用的单元数量较少，有可能会在横向变形更大的情况下达不到计算精度要求。

发明内容

本发明的目的在于提出了一种基于绝对节点坐标法的压电层合薄板大变形的仿真方法。

实现本发明目的的技术方案为：一种基于绝对节点坐标法的压电层合薄板大变形的仿真方法，包括以下步骤：

步骤1、设定压电层合薄板的几何参数和材料参数以及作动器中施加的电压大小，建立压电层合薄板模型；

步骤2、基于绝对节点坐标法对压电层合薄板进行离散，得到单元的质量矩阵、广义弹性力、广义压电力，得到压电层合薄板的动力学方程；

步骤3、利用广义α法对动力学方程求解，求得所需的变形量以及电压大小；

步骤4、通过MATLAB输出相关物理图像，包括压电层合薄板的位形图，变形和电压随时间变化的关系图。

本发明与现有技术相比，其显著优点为：

(1)本发明考虑了拉伸效应以及离散单元数量对结果的影响，使得模型更加精确；

(2)本发明使用国际先进的动力学基本理论“绝对节点坐标法”进行动力学建模和仿真，仿真结果更为准确。

附图说明

图1为本发明的压电层合薄板单元离散模型图。

图2为本发明一典型实施案例中基层分两段覆盖压电材料的模型图。

图3为本发明一典型实施案例中基层完整覆盖压电材料的模型图。

图4为本发明一典型实施案例中施加500V驱动电压求得薄板的位形图。

图5(a)为本发明一典型实施案例中施加集中力求得的薄板变形随时间变化的关系图。

图5(b)为本发明一典型实施案例中施加集中力求得的第1单元的电压随时间变化的关系图。

图5(c)为本发明一典型实施案例中施加集中力求得的第44单元的电压随时间变化的关系图。

图5(d)为本发明一典型实施案例中施加集中力求得的第96单元的电压随时间变化的关系图。

图6为本发明的流程图。

具体实施方式

一种基于绝对节点坐标法的压电层合薄板大变形的仿真方法，包括以下步骤：

步骤1、设定压电层合薄板的几何参数和材料参数以及作动器中施加的电压大小，建立如图1所示的压电层合薄板模型。

步骤2、基于绝对节点坐标法对压电层合薄板进行离散，进而推导出单元的质量矩阵M、广义弹性力Q_E、广义压电力Q_D，广义外力Q_F，最后得到压电层合薄板的动力学方程。该步骤可以细分为五个具体步骤；

步骤2.1、基于绝对节点坐标法对压电层合薄板进行离散。

在绝对节点坐标法中，节点坐标将定义在全局坐标系中。对于层合薄板单元，全局坐标将存在于三维空间中。如图1所示，板单元上任意一点的位置矢量r在全局坐标中描述为：

r＝[r₁ r₂ r₃]^T＝S_q (1)

式中，r₁，r₂和r₃分别为全局位置矢量r的三个分量，S为单元的形函数，具体形式为：

S＝[S₁I_3×3 S₂I_3×3 … S₁₂I_3×3] (2)

其中：

S₁＝-(ξ-1)(η-1)(2η²-η+2ξ²-ξ-1)

S₂＝-lξ(ξ-1)²(η-1)

S₃＝-wη(η-1)²(ξ-1)

S₄＝ξ(2η²-η-3ξ+2ξ²)(η-1)

S₅＝-lξ²(ξ-1)(η-1)

S₆＝wξη(η-1)²

S₇＝-ξη(1-3ξ-3η+2η²+2ξ²)

S₈＝lξ²η(ξ-1)

S₉＝wξη²(η-1)

S₁₀＝η(ξ-1)(2ξ²-ξ-3η+2η²)

S₁₁＝lξη²(ξ-1)

S₁₂＝-wη²(ξ-1)(η-1) (3)

其中，l代表未变形状态下板单元的长度，w代表未变形状态下板单元的宽度， I_3×3是一个大小为3×3的单位矩阵，式中ξ＝x/l，η＝y/w。x和y分别代表板单元中局部坐标沿长度方向和宽度方向的分量。

需要注意的是，本发明中所写的形函数顺序必须对应单元内的节点编号。此处假设在全局坐标中，未变形状态下，x值和y值均最小的点为单元内的1号节点，沿此节点按矩形薄板单元的逆时针顺序分别标记2号，3号，4号节点。按照此规则对单元内的节点编号才能得到与公式一致的形函数。

另外，公式中q表示单元的绝对节点坐标。对于四节点的矩形薄板单元，每个单元应该包含36个自由度，即：

其中的下标1，2，3和4分别对应板单元中的节点编号。每个节点对应的绝对节点坐标q_n应包含9个自由度，分别为：

其中，q_n矩阵内的前三个元素表示全局位置矢量r的三个分量，中间三个元素表示全局位置矢量r对局部坐标x求一阶偏导数，最后三个元素表示全局位置矢量 r对局部坐标y求一阶偏导数。

步骤2.2、求出单元的质量矩阵M。

将板单元上任意一点的位置矢量r在全局坐标中描述对时间求导，可以得到：

则层合板单元的动能可以表示为：

其中，单元质量阵可由动能导出：

M＝ρ∫_VS^TSdV (8)

另外说明，在层合板单元的动力学方程中，如果涉及积分运算，则应分成三个部分来逐层积分，分别对应其基层，作动器层和传感器层。假设矩形板单元长度为l，宽度为w，基层厚度为h_p，压电作动器厚度为h_a，压电传感器厚度为h_s，则：

公式中默认作动器粘贴在基层的上侧，而传感器粘贴在基层的下侧，如果对于一个未覆盖压电材料的单层的单元，则只需要保留公式的第一项，也即：

步骤2.3、求出单元的广义弹性力Q_E。

层合板单元的弹性势能为：

将压电材料的本构方程带入公式可得：

将公式对绝对节点坐标q求一次偏导后可得到单元的广义弹性力，表示为：

其中，K_E表示单元的刚度矩阵。将广义弹性力分为两部分Q_E1和Q_E2，其中：

根据连续介质力学中的介绍，单元内的应变表示为：

ε＝ε_m-zκ (16)

其中，拉伸方向的应变ε_m由拉格朗日应变张量导出，即：

分别将ε_m对绝对节点坐标求一次和二次偏导，可以得到公式和，分别表示为：

曲率公式具体表示为：

在公式中，有：

以及：

其中：

将n对绝对节点坐标求一次和二次偏导，可得：

将n对绝对节点坐标求一次和二次偏导，可得：

接下来即可求出曲率κ对绝对节点坐标求一次偏导：

将公式代入可以得到：

Q_E1＝h_b,a,s∫_Vε'_mc^Eε_mdV+h_1,2,3∫_Vκ'c^EκdV (29)

推导过程中一些与结构横截面几何形状有关的常数如下：

其中，下标b，a，s分别代表基层，作动器层和传感器层。

进一步地，在三维空间中，对于薄板单元，如果压电材料仅在z轴方向产生极化效应，那么压电材料的本构方程将表示为：

其中，σ表示应力矢量，ε表示应变矢量，c^E表示弹性矩阵，E表示电场强度，D

表示电位移，e表示压电常数矩阵，

表示介电常数矩阵，写为分量形式即：

其中υ代表泊松比。

对于基层的线弹性材料，可以将本构方程中的压电常数矩阵和介电常数矩阵均视为0矩阵，则可退化为线弹性材料的本构方程，即：

另外，压电层中的电压和电场强度之间的关系为：

其中，φ_s和φ_a分别代表传感器和作动器的电压；E_s和E_a分别表示传感器和作动器的电场强度。

步骤2.4、求出广义压电力Q_D。

在压电层合单元中，电势能为：

将本构方程代入公式可得：

将电势能对分别对传感器层电压φ_s和作动器层电压φ_a求导，可以得到广义压电力，其具体形式为：

步骤2.5、求出压电层合薄板结构的动力学方程。

在求得单元内的上述各项之后，对每个单元利用布尔矩阵B进行组装，其过程和有限元法相同。组装后则得到整体结构的动力学方程如下：

其中，表示M整体质量矩阵，K表示整体刚度矩阵，Q_F表示广义外力向量，Q_Ws表示传感器层的广义压电力，Q_Wa表示作动器层的广义压电力。动力学方程中的第一个式子用于求解所建立模型的运动过程，第二个和第三个式子用于求解由结构变形引起的压电材料上下端的电势差。

如果求解的问题是静力学问题，则动力学方程如下：

动力学方程中的第一个式子用于求解静力平衡位置，第二个和第三个式子用于求解平衡状态下由结构变形引起的压电材料上下端的电势差。

假设施加的广义外力为一个沿z轴的集中力，则Q_F表示为：

Q_F＝[0 0 F 0 0 0 0 0 0]^T (43)

其中F为外力大小。

步骤3、利用广义α法求解动力学方程，该算法已经是微分方程数值求解算法领域非常成熟的计算方法，故不再详细叙述。

步骤4、通过MATLAB输出相关物理图像，包括压电层合薄板的位形图，变形和电压随时间变化的关系图。

本发明基考虑了拉伸项应变，并使用了更多的单元数量，得到了更精确的结果。

实施例1

结合两个具体实施案例对本发明进行进一步的介绍，显然所描述的实施案例仅仅是本发明的某一具体实施案例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施案例，都属于本发明保护的范围。

步骤1、设定压电层合薄板的模型几何参数(如图2所示)和材料参数：薄板长度0.6m，宽度0.4m，厚度0.001m，沿长度方向的0.1～0.2m，以及0.4～ 0.5m处分别覆盖一组压电材料，一侧为传感器层，一侧为作动器层，厚度也为 0.001m；基层线弹性材料的弹性模量E＝4.9×10⁸GPa，泊松比为0.33；压电材料的弹性模量E＝3.96×10⁹GPa，泊松比为0.33，压电常数-6.948C/m²，介电常数3.01×10^-8F/m。作动器施加电压为500V，没有外力项，离散单元数量取288 个。

步骤2、基于绝对节点坐标法对压电层合薄板进行离散，得到单元的质量矩阵，广义弹性力，广义压电力，进而得到压电层合薄板的动力学方程，转入步骤 3。

步骤3、使用广义α法编写动力学求解算法，求得所需的变形量以及电压大小，转入步骤4。

步骤4、通过MATLAB输出相关物理图像，即包括压电层合薄板的位形图，如图4所示。

实施例2

步骤1、设定压电层合薄板的模型几何参数(如图2所示)和材料参数：薄板长度l＝0.6m，宽度w＝0.4m，上下两层完整铺设厚度为0.5mm的压电材料，中间层为厚度为1mm的线弹性材料。线弹性材料的弹性模量E＝7.3×10¹⁰GPa，泊松比为0.3；压电材料的弹性模量E＝1.07×10¹¹GPa，泊松比为0.3，压电常数 -6.5C/m²，介电常数3.01×10^-8F/m。离散单元为沿长度方向12个单元，沿宽度方向8个单元，共96个单元。在悬臂板末端一角施加一个沿z方向的集中力，集中力大小F＝10N，没有设置驱动电压。

步骤2、基于绝对节点坐标法对压电层合薄板进行离散，得到单元的质量矩阵，广义弹性力，广义压电力，进而得到压电层合薄板的动力学方程，转入步骤 3。

步骤3、使用广义α法编写动力学求解算法，求得所需的变形量以及电压大小，转入步骤4。

步骤4、通过MATLAB输出相关物理图像，即集中力处的横向变形随时间的关系，如图5(a)所示，第1单元，第44单元和第96单元的电压随时间的关系，如图5(b)(c)(d)所示。

Claims (8)

1.一种基于绝对节点坐标法的压电层合薄板大变形的仿真方法，其特征在于，包括以下步骤：

步骤1、设定压电层合薄板的几何参数和材料参数以及作动器中施加的电压大小，建立压电层合薄板模型；

步骤2、基于绝对节点坐标法对压电层合薄板进行离散，得到单元的质量矩阵、广义弹性力、广义压电力，得到压电层合薄板的动力学方程；

步骤3、利用广义α法对动力学方程求解，求得所需的变形量以及电压大小；

步骤4、通过MATLAB输出相关物理图像，包括压电层合薄板的位形图，变形和电压随时间变化的关系图。

2.根据权利要求1所述的基于绝对节点坐标法的压电层合薄板大变形的仿真方法，其特征在于，设定的压电层合薄板的几何参数包括基层、传感器层、作动器层的长度、宽度、高度。

3.根据权利要求1所述的基于绝对节点坐标法的压电层合薄板大变形的仿真方法，其特征在于，设定的压电层合薄板的材料参数包括：基层的弹性模量和泊松比，传感器层、作动器层的弹性模量，泊松比、压电常数、介电常数，以及压电作动器的驱动电压值或者是施加外力的大小。

4.根据权利要求1所述的基于绝对节点坐标法的压电层合薄板大变形的仿真方法，其特征在于，基于绝对节点坐标法对压电层合薄板进行离散，得到单元的质量矩阵、广义弹性力、广义压电力，得到压电层合薄板的动力学方程的具体方法为：

步骤2.1、基于绝对节点坐标法对压电层合薄板进行离散，具体为：

将层合薄板单元上任意一点的位置矢量r在全局坐标中描述为：

r＝[r₁ r₂ r₃]^T＝Sq

式中，r₁，r₂和r₃分别为全局位置矢量r的三个分量，S为单元的形函数，q表示单元的绝对节点坐标；

步骤2.2、求出单元的质量矩阵，具体为：

将层合薄板单元上任意一点的位置矢量r在全局坐标中的描述对时间求导，得到：

确定层合板单元的动能为：

式中，ρ为材料的质量密度，V为层合板单元的体积；

根据层合板单元的动能导出单元质量阵为：

M＝ρ∫_VS^TSdV

步骤2.3、求出单元的广义弹性力，具体为：

确定层合板单元的弹性势能W_E为：

式中，ε为层合板单元的应变项，σ为层合板单元的应力项；

将压电材料的本构方程带入弹性势能公式得到：

式中c^E为层合板材料的弹性矩阵，e₃₁为压电材料的压电常数，E为压电层电场强度；

上式对绝对节点坐标q求一次偏导后得到单元的广义弹性力，表示为：

其中，K_E表示单元的刚度矩阵；

步骤2.4、求出广义压电力，具体为：

确定压电层合单元中的电势能为：

式中E代表压电层的电场强度，D代表电位移向量。

将本构方程代入压电层合单元中的电势能得到：

将电势能对分别对传感器层电压φ_s和作动器层电压φ_a求导，得到广义压电力，其具体形式为：

式中。h_a和h_s分别代表作动器层和传感器层的厚度，ε₁₁和ε₂₂分别代表拉伸应变向量的第一个和第二个分量，κ₁₁和κ₂₂分别代表曲率向量的第一个和第二个分量，z代表层合板单元沿厚度方向的位置，

代表压电材料的介电常数；

步骤2.5、求出压电层合薄板结构的动力学方程；

对每个单元利用布尔矩阵B进行组装，得到整体结构的动力学方程为：

其中，M表示整体质量矩阵，K表示整体刚度矩阵，Q_F表示广义外力向量，Q_Ws表示传感器层的广义压电力，Q_Wa表示作动器层的广义压电力；

如果求解的问题是静力学问题，则动力学方程为：

5.根据权利要求4所述的基于绝对节点坐标法的压电层合薄板大变形的仿真方法，其特征在于，单元的形函数具体形式为：

S＝[S₁I_3×3 S₂I_3×3 … S₁₂I_3×3]

其中：

S₁＝-(ξ-1)(η-1)(2η²-η+2ξ²-ξ-1)

S₂＝-lξ(ξ-1)²(η-1)

S₃＝-wη(η-1)²(ξ-1)

S₄＝ξ(2η²-η-3ξ+2ξ²)(η-1)

S₅＝-lξ²(ξ-1)η-1)

S₆＝wξη(_η-1)²

S₇＝-ξη(1-3ξ-3η+2η²+2ξ²)

S₈＝lξ²η(ξ-1)

S₉＝wξη²(η-1)

S₁₀＝η(ξ-1)(2ξ²-ξ-3η+2η²)

S₁₁＝lξη²(ξ-1)

S₁₂＝-wη²(ξ-1)(η-1)

其中，l代表未变形状态下板单元的长度，w代表未变形状态下板单元的宽度，I_3×3是一个大小为3×3的单位矩阵，式中ξ＝x/l，η＝y/w，x和y分别代表板单元中局部坐标沿长度方向和宽度方向的分量。

6.根据权利要求4所述的基于绝对节点坐标法的压电层合薄板大变形的仿真方法，其特征在于，步骤2.2中的积分运算，分成三个部分来逐层积分，分别对应其基层、作动器层和传感器层，假设矩形板单元长度为l，宽度为w，基层厚度为h_p，压电作动器厚度为h_a，压电传感器厚度为h_s，则：

其中，作动器粘贴在基层的上侧，而传感器粘贴在基层的下侧，如果对于一个未覆盖压电材料的单层的单元，则：

7.根据权利要求4所述的基于绝对节点坐标法的压电层合薄板大变形的仿真方法，其特征在于，步骤2.3中所述广义弹性力包括两部分Q_E1和Q_E2，分别为：

Q_E1＝h_b，a，s∫_Vc′_mc^Eε_mdV+h_1，2，3∫_Vκ′c^EκdV

其中，ε_m和ε′_m分别代表拉伸应变向量及其对绝对节点坐标q的偏导数，κ和κ′分别代表曲率向量及其对绝对节点坐标q的偏导数，c^E代表材料的弹性矩阵，e₃₁代表压电材料的压电常数，E代表压电材料的电场强度，h代表厚度，下标b，a，s分别代表基层、作动器层和传感器层，下标1～5分别所对应的量分别为：

8.根据权利要求4所述的基于绝对节点坐标法的压电层合薄板大变形的仿真方法，其特征在于，压电材料的本构方程将表示为：

其中，σ表示应力矢量，ε表示应变矢量，c^E表示弹性矩阵，E表示电场强度，D表示电位移，e表示压电常数矩阵，

表示介电常数矩阵，写为分量形式即：

其中υ代表泊松比；

对于基层的线弹性材料，将本构方程中的压电常数矩阵和介电常数矩阵均视为0矩阵，则可退化为线弹性材料的本构方程，即：

压电层中的电压和电场强度之间的关系为：

其中，φ_s和φ_a分别代表传感器和作动器的电压；E_s和E_a分别表示传感器和作动器的电场强度。

Images