绝对节点坐标法中压电层合板的大变形研究

郭永彬，孙靖尧，黎 亮，章定国; GUO Yong-bin; SUN Jing-yao; LI Liang; ZHANG Ding-guo

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

绝对节点坐标法中压电层合板的大变形研究 PDF

- ORCID：
郭永彬
✉
- ORCID：
孙靖尧
- ORCID：
黎亮
✉
- ORCID：
章定国

南京理工大学理学院，江苏南京 210094

中图分类号： O322； O313.7

最近更新：2023-01-04

DOI：10.16385/j.cnki.issn.1004-4523.2022.06.007

摘要

压电层合板可以通过对作动器施加电压变形为各种形态，在智能可折叠结构领域具有潜在的应用前景。理解这种结构的大变形作动机理是软体智能结构设计的基础。利用等效单层模型，基于绝对节点坐标法（ANCF），建立了一种柔性压电层合薄板单元。引入了压电材料的本构方程来推导弹性力和压电力，建立了压电层合薄板单元的动力学方程，并对比了在大变形范围内三种曲率表达形式的收敛效果。通过与ABAQUS有限元软件求得的结果进行比较，验证了压电层合薄板单元的正确性，并给出了一些大变形算例。结果表明，在压电材料全覆盖、部分覆盖、多段覆盖等情形下，该单元均能得到较为稳定的结果，表明此单元可以与单层线弹性单元进行耦合。此外还研究了压电层合板受集中力状态下的动态响应。该研究有助于理解受智能压电材料驱动的柔性或软体结构的复杂耦合非线性力学行为。

关键词

非线性振动; 绝对节点坐标法; 层合板; 压电材料; 智能结构; 大变形

1 概述

压电材料是一种同时具备正逆压电效应的智能材料，可被用作传感器和作动器，并被广泛用于振动控制、压电驱动、健康监测、能量收集等领域。由于近年来对轻型结构的需求不断增多，在上述的应用中，特别是压电驱动方面，包含大变形的压电层合结构具有越来越重要的研究价值。当结构的表面被压电复合材料覆盖并且施加外部电压到作动器时，可折叠的智能结构（例如折纸结构和软体机器人）可以变形到所需的形状。

目前，已有大量论文进行过包含压电材料的结构的非线性分析^［

1⁃6］，这些研究工作主要关注传统的压电陶瓷。实际上，压电陶瓷通常具有较高的结构刚度，这使其在粘结过程中（尤其是在粘贴弯曲表面上时）易于损坏。而压电纤维基复合材料不仅具有柔韧性，还具有很高的驱动力，因此受到了研究人员的关注。即使经历大的变形，压电复合材料仍可以紧密地贴合基层的表面。郭翔鹰等^{［参考文献 7

百度学术}7］研究了压电纤维复合材料（MFC）层合壳的非线性动力学问题。研究人员还证实，使用基于压电纤维的复合材料（如PVDF和MFC）在结构振动抑制方面的性能要优于传统的压电陶瓷^{［参考文献 7⁃8}7⁃8］。上述工作采用的理论均是基于小变形假设，因此并不适用于研究结构的大变形问题。

绝对节点坐标法（ANCF法）由Shabana等建立^［

9⁃10］，与假设模态法以及传统有限元法相比，ANCF法更适合柔性系统大变形建模。在过去的20年中，大量学者尝试了不同的新单元来描述更为复杂的机械模型^{［参考文献 11

百度学术}11］或解决早期研究中出现的诸如泊松锁定和剪切锁定等问题^{［参考文献 12

百度学术}12］，并提高计算的精度^{［参考文献 13

百度学术}13］。ANCF法具有常质量矩阵，零离心力和科氏力等特点，并且可以精确描述刚体的运动^{［参考文献 14

百度学术}14］。唯一形式复杂的部分是刚度阵，它是由连续介质力学或弹性线方法推导出的^{［参考文献 15

百度学术}15］。随着ANCF的广泛应用，该方法也可用于创建更精确的模型来分析复合结构的力学特性。张炜华等^{［参考文献 16

百度学术}16］对大变形复合材料薄板的多体系统动力学建模方法进行了研究。游斌弟等^{［参考文献 17

百度学术}17］进行了空间望远镜层合材料镜片展开过程非线性动力学行为分析。目前对于层合结构的位移场，主要有两种描述方法：等效单层（ESL）模型和分层模型^{［参考文献 18

百度学术}18］。在ESL模型中，层合结构的变形将被视为与单层结构相同。不同层的运动和变形都使用相同的自由度，该自由度一般在多层结构的中间轴上。这种建模方法的优点是，沿厚度方向的位移场是平滑的，并且无需层与层之间的边界条件，因此计算效率大为提高^{［参考文献 19

百度学术}19］。在分层模型中，位移场在每个层中都是平滑的，但在层的界面处不是平滑的。ESL模型和分层模型的示意图如图1所示。使用分层模型得到的结果一般更准确，但由于自由度过多，会影响计算效率。

图1 层合结构模型

Fig.1 The model of laminated structures

基于ANCF法的ESL模型也可用于小变形或大变形的压电层合结构。Gilardi等^［

20］基于ANCF法建立了压电作动器覆盖的柔性梁模型，研究了其静态和动态特性，并优化了高速旋转下的振动控制效果。Jiang等^{［参考文献 21

百度学术}21］基于高阶位移场，二次电势场和线性温度场建立了压电层合板结构的数学模型，并提供了数值模拟计算。Liu等^{［参考文献 22

百度学术}22］提出了一种新的计算方法用于由层合板组成的大型刚柔多体系统的动态分析。文献［23⁃24］基于ANCF法建立了多层压电薄板单元模型，并分析了由压电作动器驱动的三层薄板结构的振动特性。将计算的数值结果与商业有限元软件COMSOL的仿真结果进行对比，结论较为一致。但是文献［23⁃24］在推导过程的弹性力项中忽略了拉伸项应变，这使得该单元只能用于低张力板；同时，其建立的层合板离散模型使用的单元数量较少，有可能会在横向变形更大的情况下达不到计算精度要求。易灿明等^{［参考文献 25

百度学术}25］基于绝对节点坐标法，通过变形协调条件，建立了以基层为薄板单元，压电层为梁单元的梁板耦合模型。

本文采用ANCF法，在文献［

23⁃24］的基础上建立了一种考虑板的拉伸项应变的压电层合薄板单元，推导了该单元的动力学方程，并进行了几组静力学和动力学仿真。结论指出，压电层合薄板单元也可以很好地和线弹性单元耦合，即可以表达压电材料覆盖在层合薄板上任意位置的情况。

2 动力学模型

2.1 缩减板单元

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

图2 层合薄板单元模型

Fig.2 The model of thin laminated plate

r = {[\begin{matrix} r_{1} & r_{2} & r_{3} \end{matrix}]}^{T} = S q

（1）

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

S = [\begin{matrix} S_{1} I_{3 \times 3} & S_{2} I_{3 \times 3} & \dots & S_{12} I_{3 \times 3} \end{matrix}]

（2）

其中：

\begin{array}{l} S_{1} = - (ξ - 1) (η - 1) (2 η^{2} - η + 2 ξ^{2} - ξ - 1), \\ S_{2} = - l ξ {(ξ - 1)}^{2} (η - 1), \\ S_{3} = - w η {(η - 1)}^{2} (ξ - 1), \\ S_{4} = ξ (2 η^{2} - η - 3 ξ + 2 ξ^{2}) (η - 1), \\ S_{5} = - l ξ^{2} (ξ - 1) (η - 1), \\ S_{6} = w ξ η {(η - 1)}^{2}, \\ S_{7} = - ξ η (1 - 3 ξ - 3 η + 2 η^{2} + 2 ξ^{2}), \\ S_{8} = l ξ^{2} η (ξ - 1), \\ S_{9} = w ξ η^{2} (η - 1), \\ S_{10} = η (ξ - 1) (2 ξ^{2} - ξ - 3 η + 2 η^{2}), \\ S_{11} = l ξ η^{2} (ξ - 1), \\ S_{12} = - w η^{2} (ξ - 1) (η - 1) \end{array}

（3）

式中 l和w分别代表未变形状态下板单元的长度和宽度；I_3×3是一个大小为3×3的单位矩阵。ξ=x/l，η=y/w。x和y分别代表板单元未变形时中面上任意一点P在板单元坐标系下沿长度方向和宽度方向的坐标。

对于四节点缩减板单元，每个绝对坐标应包含9个自由度，分别为：

q_{n o d e} = {[r_{1} r_{2} r_{3} \frac{\partial r_{1}}{\partial x} \frac{\partial r_{2}}{\partial x} \frac{\partial r_{3}}{\partial x} \frac{\partial r_{1}}{\partial y} \frac{\partial r_{2}}{\partial y} \frac{\partial r_{3}}{\partial y}]}^{T}

（4）

这样每个单元应包含36个自由度。

2.2 等效单层模型

在对层合单元的积分过程中，需要分成三个部分来逐层积分。假设板单元长度为l，宽度为w，基层厚度为h_p，压电作动器厚度为h_a，压电传感器厚度为h_s，则：

\int_{V} d V = \int_{- \frac{h_{p}}{2}}^{\frac{h_{p}}{2}} \int_{0}^{l} \int_{0}^{w} d V_{p} + \int_{\frac{h_{p}}{2}}^{\frac{h_{p}}{2} + h_{a}} \int_{0}^{l} \int_{0}^{w} d V_{a} + \int_{- \frac{h_{p}}{2} - h_{s}}^{- \frac{h_{p}}{2}} \int_{0}^{l} \int_{0}^{w} d V_{s}

（5）

式中　 $V_{p}$ ， $V_{a}$ 和 $V_{s}$ 分别为板单元中基层、压电作动器和压电传感器的单元体积。

如果是单层板单元，则只需要公式（5）中的第一项，即：

\int_{V} d V = \int_{- \frac{h_{p}}{2}}^{\frac{h_{p}}{2}} \int_{0}^{l} \int_{0}^{w} d V_{p}

（6）

2.3 压电材料的本构方程

在三维空间中，如果压电材料仅在z轴方向产生极化效应，那么压电材料的本构方程将表示为：

\{\begin{matrix} σ = c^{E} ε - e E \\ D = e^{T} ε + ς E \end{matrix}

（7）

式中 σ表示应力矢量；ε表示应变矢量；c^E表示弹性矩阵；E表示电场强度；D表示电位移；e表示压电常数矩阵； $ς$ 表示介电常数矩阵。

各分量写成矩阵或向量形式为：

σ = {[\begin{matrix} σ_{x x} & σ_{y y} & σ_{x y} \end{matrix}]}^{T}

（8）

ε = {[\begin{matrix} ε_{x x} & ε_{y y} & ε_{x y} \end{matrix}]}^{T}

（9）

c^{E} = \frac{c_{E}}{1 - υ^{2}} [\begin{matrix} 1 & υ & 0 \\ υ & 1 & 0 \\ 0 & 0 & (1 - υ) / 2 \end{matrix}]

（10）

式中 υ和 $c_{E}$ 分别代表压电材料的泊松比和弹性模量。

e = [\begin{matrix} 0 & 0 & e_{31} \\ 0 & 0 & e_{31} \\ 0 & e_{15} & 0 \end{matrix}]

（11）

E = {[\begin{matrix} 0 & 0 & E \end{matrix}]}^{T}

（12）

D = {[\begin{matrix} D_{x x} & D_{y y} & D_{x y} \end{matrix}]}^{T}

（13）

ς = [\begin{matrix} ς_{11} & 0 & 0 \\ 0 & ς_{11} & 0 \\ 0 & 0 & ς_{33} \end{matrix}]

（14）

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

[\begin{matrix} σ_{x x} \\ σ_{y y} \\ σ_{x y} \end{matrix}] = \frac{c_{E}}{1 - υ^{2}} [\begin{matrix} 1 & υ & 0 \\ υ & 1 & 0 \\ 0 & 0 & (1 - υ) / 2 \end{matrix}] [\begin{matrix} ε_{x x} \\ ε_{y y} \\ ε_{x y} \end{matrix}]

（15）

另外，需要注意的是，当基层弯曲时，覆盖在顶层的压电传感器中的电荷将在内部移动，从而导致压电传感器的上端和下端的电荷分布不均匀，传感器将收集此电压信号作为反馈。同时，如果对作动器施加电压则会导致基层弯曲。压电层中的电压和电场强度之间的关系为：

E_{s} = \frac{ϕ_{s}}{h_{s}}, E_{a} = \frac{ϕ_{a}}{h_{a}}

（16）

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

2.4 广义质量阵

将公式（1）对时间求导，可以得到：

\dot{r} = S \dot{q}

（17）

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

T_{e} = \frac{1}{2} ρ \int_{V} {\dot{r}}^{T} \dot{r} d V = \frac{1}{2} {\dot{q}}^{T} M_{e} \dot{q}

（18）

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

M_{e} = ρ \int_{V} S^{T} S d V

（19）

式中　 $ρ$ 为板单元的等效密度。

2.5 广义弹性力

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

{W_{E}}_{e} = \frac{1}{2} \int_{V} ε σ d V

（20）

将压电材料的本构方程代入公式（20）可得：

{W_{E}}_{e} = \frac{1}{2} \int_{V} ε c^{E} ε d V - \frac{1}{2} \int_{V} ε e_{31} E d V

（21）

将公式（21）对广义坐标q求导后可得：

\begin{array}{l} {Q_{E}}_{e} = \frac{\partial W_{E e}}{\partial q} = \int_{V} \frac{\partial ε}{\partial q} c^{E} ε d V - \\ \frac{1}{2} \int_{V} \frac{\partial ε}{\partial q} e_{31} E d V = Q_{E e 1} - Q_{E e 2} \end{array}

（22）

其中

Q_{E e 1} = \int_{V} \frac{\partial ε}{\partial q} c^{E} ε d V

（23）

Q_{E e 2} = \frac{1}{2} \int_{V} \frac{\partial ε}{\partial q} e_{31} E d V

（24）

根据连续介质力学，应变表示为：

ε = ε_{m} - z κ

（25）

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

ε_{m} = \frac{1}{2} [\begin{matrix} ε_{11} \\ ε_{22} \\ ε_{12} \end{matrix}] = \frac{1}{2} [\begin{matrix} r_{x}^{T} r_{x} - 1 \\ r_{y}^{T} r_{y} - 1 \\ 2 r_{x}^{T} r_{y} \end{matrix}] = \frac{1}{2} [\begin{matrix} q^{T} S_{x}^{T} S_{x} q - 1 \\ q^{T} S_{y}^{T} S_{y} q - 1 \\ 2 q^{T} S_{x}^{T} S_{y} q \end{matrix}]

（26）

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

ε'_{m} = \frac{\partial ε_{m}}{\partial q} = [\begin{matrix} q^{T} S_{x}^{T} S_{x} \\ q^{T} S_{y}^{T} S_{y} \\ q^{T} S_{x}^{T} S_{y} + q^{T} S_{y}^{T} S_{x} \end{matrix}]

（27）

ε''_{m} = \frac{\partial ε'_{m}}{\partial q} = [\begin{matrix} S_{x}^{T} S_{x} \\ S_{y}^{T} S_{y} \\ S_{x}^{T} S_{y} + S_{y}^{T} S_{x} \end{matrix}]

（28）

弯曲方向的应变由曲率的定义式导出，其中精确的曲率公式为：

κ = [\begin{matrix} κ_{11} \\ κ_{22} \\ κ_{12} \end{matrix}] = [\begin{matrix} (r_{x x} n) / {\bar{n}}^{3} \\ (r_{y y} n) / {\bar{n}}^{3} \\ 2 (r_{x y} n) / {\bar{n}}^{3} \end{matrix}]

（29）

为了简化计算，曲率可以简化为以下形式，本文中称之为一次简化公式：

κ = [\begin{matrix} κ_{11} \\ κ_{22} \\ κ_{12} \end{matrix}] = [\begin{matrix} (r_{x x} n) / \bar{n} \\ (r_{y y} n) / \bar{n} \\ 2 (r_{x y} n) / \bar{n} \end{matrix}]

（30）

或者简化为曲率和绝对节点坐标之间的线性关系，本文中称之为线性简化公式：

κ = [\begin{matrix} κ_{11} \\ κ_{22} \\ κ_{12} \end{matrix}] = [\begin{matrix} S_{x x} q \\ S_{y y} q \\ 2 S_{x y} q \end{matrix}]

（31）

接下来以精确的曲率公式为例进行推导，其他的形式同理。在公式（29）中，有：

n = r_{x} \times r_{y} = [\begin{matrix} q^{T} S_{x 2}^{T} S_{y 3} q - q^{T} S_{x 3}^{T} S_{y 2} q \\ q^{T} S_{x 3}^{T} S_{y 1} q - q^{T} S_{x 1}^{T} S_{y 3} q \\ q^{T} S_{x 1}^{T} S_{y 2} q - q^{T} S_{x 2}^{T} S_{y 1} q \end{matrix}] = [\begin{matrix} q^{T} S_{Δ 1} q \\ q^{T} S_{Δ 2} q \\ q^{T} S_{Δ 3} q \end{matrix}] = [\begin{matrix} {\bar{n}}_{1} \\ {\bar{n}}_{2} \\ {\bar{n}}_{3} \end{matrix}]

（32）

以及

\bar{n} = \sqrt[]{n^{T} n} = \sqrt[]{{\bar{n}}_{1}^{2} + {\bar{n}}_{2}^{2} + {\bar{n}}_{3}^{2}}

（33）

其中

\begin{array}{l} S_{Δ 1} = \frac{1}{2} [(S_{x 2}^{T} S_{y 3} - S_{x 3}^{T} S_{y 2}) + \\ (S_{x 2}^{T} S_{y 3} - S_{x 3}^{T} S_{y 2})^{T}], \\ S_{Δ 2} = \frac{1}{2} [(S_{x 3}^{T} S_{y 1} - S_{x 1}^{T} S_{y 3}) + \\ (S_{x 3}^{T} S_{y 1} - S_{x 1}^{T} S_{y 3})^{T}], \\ S_{Δ 3} = \frac{1}{2} [(S_{x 1}^{T} S_{y 2} - S_{x 2}^{T} S_{y 1}) + \\ (S_{x 1}^{T} S_{y 2} - S_{x 2}^{T} S_{y 1})^{T}] \end{array}

（34）

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

n' = \frac{\partial n}{\partial q} = [\begin{matrix} 2 q^{T} S_{Δ 1} \\ 2 q^{T} S_{Δ 2} \\ 2 q^{T} S_{Δ 3} \end{matrix}]

（35）

n'' = \frac{\partial n'}{\partial q} = [\begin{matrix} 2 S_{Δ 1} \\ 2 S_{Δ 2} \\ 2 S_{Δ 3} \end{matrix}]

（36）

将 $\bar{n}$ 对绝对节点坐标求一次和二次偏导，可得：

\bar{n}' = \frac{\partial \bar{n}}{\partial q} = \frac{n^{T} n'}{\bar{n}}

（37）

\bar{n}'' = \frac{\partial \bar{n}'}{\partial q} = \frac{n'^{T} n' + n^{T} n''}{\bar{n}} - \frac{\bar{n}'^{T} \bar{n}'}{\bar{n}}

（38）

曲率κ对绝对节点坐标求一次偏导为：

\begin{array}{l} κ' = \frac{\partial κ}{\partial q} = [\begin{matrix} κ'_{11} \\ κ'_{22} \\ κ'_{12} \end{matrix}] = \\ [\begin{matrix} (S_{x x}^{T} n + r_{x x}^{T} n') / {\bar{n}}^{3} - 3 κ_{11} \bar{n}' / \bar{n} \\ (S_{y y}^{T} n + r_{y y}^{T} n') / {\bar{n}}^{3} - 3 κ_{22} \bar{n}' / \bar{n} \\ 2 (S_{x y}^{T} n + r_{x y}^{T} n') / {\bar{n}}^{3} - 3 κ_{12} \bar{n}' / \bar{n} \end{matrix}] \end{array}

（39）

将公式（26）~（39）代入（23）和（24）可以得到：

\begin{array}{l} Q_{E e 1} = h_{a} \int_{V} ε'_{m} c^{E} ε_{m} d V + h_{1} \int_{V} κ' c^{E} κ d V + \\ h_{b} \int_{V} ε'_{m} c^{E} ε_{m} d V + h_{2} \int_{V} κ' c^{E} κ d V + \\ h_{s} \int_{V} ε'_{m} c^{E} ε_{m} d V + h_{3} \int_{V} κ' c^{E} κ d V \end{array}

（40）

\begin{array}{l} Q_{E e 2} = \frac{1}{2} h_{a} \int_{V} ε'_{m} e_{31} E d V - \frac{1}{2} h_{4} \int_{V} κ' e_{31} E d V + \\ \frac{1}{2} h_{s} \int_{V} ε'_{m} e_{31} E d V - \frac{1}{2} h_{5} \int_{V} κ' e_{31} E d V \end{array}

（41）

其中，公式（41）还可以进一步化简为：

\begin{array}{l} Q_{E e 2} = \frac{1}{2} \int_{V} (ε'_{11} + ε'_{22}) e_{31} ϕ_{a} d V + \\ \frac{1}{2} \int_{V} (ε'_{11} + ε'_{22}) e_{31} ϕ_{s} d V - \\ \frac{1}{2} h_{4,5} \int_{V} (κ'_{11} + κ'_{22}) e_{31} ϕ_{a} d V - \\ \frac{1}{2} h_{4,5} \int_{V} (κ'_{11} + κ'_{22}) e_{31} ϕ_{s} d V \end{array}

（42）

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

\begin{array}{l} h_{1} = h_{a} {(\frac{h_{p}}{2} + \frac{h_{a}}{2})}^{2} + \frac{h_{a}^{3}}{12}, \\ h_{2} = \frac{h_{p}^{3}}{12}, \\ h_{3} = h_{s} {(\frac{h_{p}}{2} + \frac{h_{s}}{2})}^{2} + \frac{h_{s}^{3}}{12} \end{array}

（43）

\begin{array}{l} h_{4} = \frac{1}{2} (h_{p} + h_{a}), \\ h_{5} = \frac{1}{2} (h_{p} + h_{s}) \end{array}

（44）

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

2.6 广义压电力

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

W_{D e} = \frac{1}{2} \int_{V} E^{T} D d V

（45）

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

W_{D e} = \frac{1}{2} \int_{V} E^{T} e^{T} ε d V + \frac{1}{2} \int_{V} E^{T} ς E d V

（46）

将电势能分别对传感器层电压ϕ_s和作动器层电压ϕ_a求导，可以得到广义压电力：

Q_{D a e} = \frac{e_{31}}{2 h_{a}} \int_{V} [(ε_{11} + ε_{22}) - z (κ_{11} + κ_{22})] d V + \int_{V} \frac{1}{h_{a}} ς_{33} \frac{ϕ_{a}}{h_{a}} d V

（47）

Q_{D s e} = \frac{e_{31}}{2 h_{s}} \int_{V} [(ε_{11} + ε_{22}) - z (κ_{11} + κ_{22})] d V + \int_{V} \frac{1}{h_{s}} ς_{33} \frac{ϕ_{s}}{h_{s}} d V

（48）

如果不考虑外界在作动器上输入的电压，则有Q_Dse=0以及Q_Dae=0，这样即可对公式（47）和（48）求解，得出由结构变形引起的压电材料上下端的电势差。

2.7 动力学方程

在求得单元内的上述各项之后，需要对每个单元进行组装。

M = \sum_{e = 1}^{n} B_{e}^{T} M_{e} B_{e}

（49）

Q_{E} = \sum_{e = 1}^{n} B_{e}^{T} {Q_{E}}_{e}

（50）

Q_{F} = \sum_{e = 1}^{n} B_{e}^{T} {Q_{F}}_{e}

（51）

Q_{D s} = \sum_{e = 1}^{n} B_{e}^{T} {Q_{D s}}_{e}

（52）

Q_{D a} = \sum_{e = 1}^{n} B_{e}^{T} {Q_{D a}}_{e}

（53）

式中 B_e表示薄板单元节点坐标向量对应总体绝对节点坐标列阵中的布尔定位矩阵；M表示整体质量矩阵； $Q_{E}$ 表示板的整体广义弹性力向量；Q_F表示整体广义外力向量；Q_Ds表示传感器层的整体广义压电力；Q_Da表示作动器层的整体广义压电力。

组装后则得到整体结构的动力学方程如下：

\{\begin{array}{l} M \ddot{q} + Q_{E} = Q_{F} \\ Q_{D s} = 0 \\ Q_{D a} = 0 \end{array}

（54）

如果求解静力学问题，则动力学方程中将不包含 $M \ddot{q}$ 这一项。

3 仿真算例

3.1 不同曲率计算形式的收敛性分析

在进行压电层合结构的相关仿真之前，有必要先明确在大变形范围内三种不同的曲率表达形式的适用性。如图3所示层合悬臂板，长度l=0.5 m，宽度w=0.15 m，厚度分别为中层0.4 mm，上下层各0.3 mm。所有层均为相同的线弹性材料，弹性模量E=207 GPa，泊松比为0.3。板的两个端部分别受到两个大小为F=15 N的集中力外载荷。分别取单元数为1~48个的情况，求解出集中力施加位置处薄板的横向位移，研究三种不同的曲率表达形式下各自的收敛性，其结果如图4所示。同时与商业有限元软件ABAQUS的结论进行对比，ABAQUS仿真中网格尺寸划分大小为0.005 m，仿真结果显示板的末端横向位移为-0.2868 m，结果如图5所示。

图3 三层线弹性板模型

Fig.3 Model of three-layer linear elastic laminated plate

图4 三种曲率计算公式在大变形状态下各自的收敛性

Fig.4 Convergence using three different forms of curvature in large deformation

图5 ABAQUS中的薄板位形图

Fig.5 The configuration of the plate in ABAQUS

从图4中可以看出，使用精确的曲率公式和分母简化为一次项的公式得出了一致的结果，而且都可以很好的收敛，并且该结论与ABAQUS得到的结论几乎相同，相对误差为1.3%。在单元数量为48个的情况下，使用精确的曲率公式仿真用时为67.54 s，使用分母简化为一次项的公式仿真用时为65.23 s，在相同的计算结果下计算效率有小幅度的提升。使用线性简化形式的曲率公式，虽然结果可以收敛，但产生的横向变形相比于真实值明显偏大。由此可见，线性简化公式是难以在大变形的条件下使用的。

3.2 薄板在集中力作用下产生大变形

算例2验证了压电层合薄板模型在不同大小的外力作用下的正确性。如图6所示的三层悬臂板，长度l=0.6 m，宽度w=0.4 m，上下两层完整铺设厚度为0.5 mm的压电材料，中间层为1 mm厚的线弹性材料。线弹性材料的弹性模量为73 GPa，泊松比为0.3；压电材料的弹性模量为107 GPa，泊松比为0.3，压电常数为-6.5 C/m²，介电常数为3.01×10^-8 F/m。离散单元为沿长度方向6个单元，沿宽度方向4个单元，共24个单元。本文所有算例中的单元编号规则均与图6中所示相同。在悬臂板末端一角施加一个沿z方向的集中力，其大小与该点横向位移的关系如图7所示。

图6 全覆盖的压电层合板模型

Fig.6 Laminated plate with fully covered piezoelectric patches

图7 不同集中力下薄板末端的横向位移

Fig.7 The tip displacement in different concentrated forces

从图7中可以看出，使用线性简化的曲率公式得到的位移也是线性变化的。在不超过30 N的范围内，使用线性简化的曲率可以和精确曲率的方法得出近乎一致的结果，且相对误差不会超过1%。但随着集中力的增大，使用线性简化的曲率计算出的结果会逐渐偏离理论值。

3.3 薄板在作动器施加电压情况下产生变形

算例3验证了压电层合薄板模型在外电压作用下的正确性，同时探索了部分覆盖压电材料的情况下该单元与普通线弹性材料单元耦合的可行性。基层板的尺寸和算例2中相同，仅第6，7个单元是压电层合单元，其余均为单层线弹性材料。材料参数与算例2相同。在两个覆盖压电材料的单元的作动器上施加逐渐增大的电压，计算板上节点13和35的横向位移（如图8所示），其对比结果如表1和表2所示。

图8 部分覆盖的压电层合板模型

Fig.8 Laminated plate with partly covered piezoelectric patches

表1 不同电压下节点13的横向位移

Tab.1 Displacement of the 13th node in different voltages

电压/V	ABAQUS/m	本文/m	误差/%
100	7.86×10^-6	8.24×10^-6	4.83
200	1.55×10^-5	1.65×10^-5	6.45
300	2.31×10^-5	2.48×10^-5	7.36
400	3.05×10^-5	3.31×10^-5	8.52

表2 不同电压下节点35的横向位移

Tab.2 Displacement of the 35th node in different voltages

电压/V	ABAQUS/m	本文/m	误差/%
100	-6.91×10^-5	-7.14×10^-5	3.18
200	-1.42×10^-4	-1.44×10^-4	1.40
300	-2.17×10^-4	-2.19×10^-4	0.92
400	-2.96×10^-4	-2.95×10^-4	-0.34

接下来讨论压电层合单元在上下两侧均覆盖作动器的情况，基层和压电材料厚度均为1 mm。层合薄板结构在沿长度方向的0.1~0.3 m，沿宽度方向的0~0.4 m处覆盖压电材料。基层线弹性材料的弹性模量E=0.49 GPa，泊松比为0.33；压电材料的弹性模量E=39.6 GPa，泊松比为0.33，压电常数为-6.948 C/m²，介电常数为3.01×10^-8 F/m。考虑到压电作动的情况下随着变形量数量级的增大，需要更多的单元数量才能接近收敛值，此处离散单元数量为沿长度方向24个，沿宽度方向12个，共288个单元。在作动器施加0~1000 V逐渐增大的电压，求解出薄板自由端角落位置（即图9中A位置）的横向位移如图10所示。可以看到，在施加电压的情况下，变形趋势也与施加集中力的情况相同，呈现一种非线性趋势。

图9 双侧部分覆盖的压电层合板模型

Fig.9 Laminated plate with segmented piezoelectric patches on two sides

图10 不同电压下薄板末端的横向位移

Fig.10 The tip displacement in different voltages

最后计算一组分段贴片的情况。层合薄板结构在沿长度方向的0.1~0.2 m，以及0.4~0.5 m处分别覆盖一组压电材料，一侧为传感器层，一侧为作动器层，施加电压为500 V，模型如图11所示。基层线弹性材料的弹性模量E=0.49 GPa，泊松比为0.33；压电材料的弹性模量E=3.96 GPa，泊松比为0.33，压电常数为-6.948 C/m²，介电常数为3.01×10^-8 F/m。离散单元同样取288个。变形后基层线弹性材料薄板部分的三维位形图如图12所示。

图11 分段覆盖的压电层合板模型

Fig.11 Laminated plate with multiple piezoelectric patches

图12 施加500 V电压下层合板的位形图

Fig.12 The configuration of laminated plate under 500 V

从图12中可以看出，层合薄板在受到电压驱动的状态下，变形方式并非简单地围绕y轴方向进行卷曲，而是朝向类似球形的趋势产生变形，导致矩形的薄板在电压作用下远端的角落处变形比中间处更大，而在薄板中心位置有一处明显的凹陷。

3.4 压电层合板动力学仿真

对算例4进行动力学仿真。压电层合薄板的压电材料为全覆盖，线弹性材料和压电材料参数与算例1相同。离散单元数量为沿x轴12个，沿y轴8个，共96个单元。仿真时间为1 s，求解出仿真时间内集中力处的横向位移随时间的变化如图13所示。第1单元，第44单元和第96单元的电压随时间的变化如图14所示。

图13 压电层合板在集中力作用下的位移随时间变化情况

Fig.13 The displacement of piezoelectric laminated plates varies with time under the concentrated force

图14 不同单元电压随时间变化情况

Fig.14 The voltages vary with time in different elements

通过图13和14可以看出，在越靠近约束端的位置，集中力作用下产生的电压越大，而在施加集中力的自由端电压值反而较小。电压的增减趋势与横向位移的变化具有大致相同的周期。

4 结论

本文基于绝对节点坐标法建立了一个压电层合薄板单元，并通过几组算例验证了单元的正确性，以及使用不同形式的曲率计算公式下的收敛性。在大变形范围内，线性的曲率计算方法将不再适用，其计算结果相比于真实值偏大。当电压施加到压电板单元中的作动器时，层合薄板将弯曲，其值与ABAQUS有限元软件提供结果基本一致，这表明本模型描述压电耦合效应时具有很好的实用性，且在压电材料全覆盖或是部分、分段覆盖的情况下皆可适用。

参考文献

Abdelkefi A， Nayfeh A H， Hajj M R. Effects of nonlinear piezoelectric coupling on energy harvesters under direct excitation［J］. Nonlinear Dynamics， 2012， 67（2）： 1221-1232. [百度学术]

Leadenham S， Erturk A. Unified nonlinear electroelastic dynamics of a bimorph piezoelectric cantilever for energy harvesting， sensing， and actuation［J］. Nonlinear Dynamics， 2015， 79（3）： 1727-1743. [百度学术]

Latalski J， Warminski J. Nonlinear vibrations of a rotating thin-walled composite piezo-beam with circumferentially uniform stiffness （CUS）［J］. Nonlinear Dynamics， 2019， 98（4）： 2509-2529. [百度学术]

Tajeddini V， Muliana A. Nonlinear deformations of piezoelectric composite beams［J］. Composite Structures， 2015， 132： 1085-1093. [百度学术]

Wang B T， Rogers C A. Modeling of finite-length spatially distributed induced strain actuators for laminated beams and plates［J］. Intelligent Materials and Systems， 1991， 2（1）： 38-58. [百度学术]

李华东，齐克奇，石强，等. 应用于空间惯性传感器的压电驱动器的性能测试［J］. 中国光学， 2020， 13（2）： 344-353. [百度学术]

Li Huadong， Qi Keqi， Shi Qiang， et al. Performance test of piezoelectric actuators for space inertial sensors［J］. Chinese Optics， 2020， 13（2）：344-353. [百度学术]

郭翔鹰，刘大猛，张伟. 压电纤维复合材料层合壳的非线性动力学研究［J］. 动力学与控制学报， 2017， 15（5）： 430-438. [百度学术]

Guo Xiangying， Liu Dameng， Zhang Wei. Nonlinear dynamic analysis of piezoelectric macro-fiber composite laminated shells［J］. Journal of Dynamics and Control， 2017， 15（5）： 430-438. [百度学术]

Acevedo P， Vázquez M， Durán J， et al. The measurement of temperature gradients in a soft tissue phantom using PVDF arrays： a simulation case using the finite element method （FEM）［J］. Physics Procedia， 2015， 70： 1100-1103. [百度学术]

Berzeri M， Shabana A A. Development of simple models for the elastic forces in the absolute nodal coordinate formulation［J］. Journal of Sound and Vibration， 2000， 235（4）： 539-565. [百度学术]

Shabana A A， Christensen A P. Three-dimensional absolute nodal coordinate formulation： plate problem ［J］. International Journal for Numerical Methods in Engineering， 1997， 40（15）： 2775-2790. [百度学术]

Kübler L， Eberhard P， Geisler J， et al. Flexible multibody systems with large deformations and nonlinear structural damping using absolute nodal coordinates［J］. Nonlinear Dynamics， 2003， 34（1）： 31-52. [百度学术]

Nachbagauer K. State of the art of ANCF elements regarding geometric description， interpolation strategies， definition of elastic forces， validation and the locking phenomenon in comparison with proposed beam finite elements［J］. Archives of Computational Methods in Engineering， 2014， 21（3）： 293-319. [百度学术]

章孝顺，章定国，陈思佳，等. 基于绝对节点坐标法的大变形柔性梁几种动力学模型研究［J］. 物理学报， 2016， 65（9）： 156-165. [百度学术]

Zhang Xiaoshun， Zhang Dingguo， Chen Sijia， et al. Several dynamic models of a large deformation flexible beam based on the absolute nodal coordinate formulation［J］. Acta Physica Sinica， 2016， 65（9）：156-165. [百度学术]

Langerholc M， Slavič J， Boltežar M. A thick anisotropic plate element in the framework of an absolute nodal coordinate formulation［J］. Nonlinear Dynamics， 2013， 73（1-2）： 183-198. [百度学术]

Gerstmayr J， Irschik H. On the correct representation of bending and axial deformation in the absolute nodal coordinate formulation with an elastic line approach［J］. Journal of Sound and Vibration， 2008， 318（3）： 461-487. [百度学术]

张炜华，刘锦阳. 大变形复合材料薄板多体系统动力学建模［J］. 振动与冲击， 2016， 35（8）： 27-35. [百度学术]

Zhang Weihua， Liu Jinyang. Dynamic modeling of composite thin-plate multibody systems with large deformation［J］. Journal of Vibration and Shock， 2016， 35（8）： 27-35. [百度学术]

游斌弟，梁栋，温晓雷，等. 空间望远镜层合材料镜片展开过程非线性动力学行为分析［J］. 振动与冲击， 2020， 39（2）： 1-8. [百度学术]

You Bindi， Liang Dong， Wen Xiaolei， et al. Nonlinear dynamics behavior of a space telescope flexible reflective mirror in the deployment process considering the effects of geometric nonlinearity and laminated composite material ［J］. Journal of Vibration and Shock， 2020， 39（2）： 1-8. [百度学术]

Saravanos D A， Heyliger P R. Mechanics and computational models for laminated piezoelectric beams， plates， and shells［J］. ASME Applied Mechanics Reviews， 1999， 52（10）： 305-319. [百度学术]

Carrera E. Theories and finite elements for multilayered， anisotropic， composite plates and shells［J］. Archives of Computational Methods in Engineering， 2002， 9（2）： 87-140. [百度学术]

Gilardi G， Buckham B J， Park E J. Finite element modeling of a highly flexible rotating beam for active vibration suppression with piezoelectric actuators［C］. ASME International Mechanical Engineering Congress and Exposition. Seattle， USA， 2007： 2041-2050. [百度学术]

Jiang J P， Li D X. Finite element formulations for thermopiezoelastic laminated composite plates［J］. Smart Materials and Structures， 2008， 17（1）： 015027. [百度学术]

Liu C， Tian Q， Hu H Y. Dynamics of a large scale rigid-flexible multibody system composed of composite laminated plates［J］. Multibody System Dynamics， 2011， 26（3）： 283-305. [百度学术]

Nada A A， El-Assal A M. Absolute nodal coordinate formulation of large-deformation piezoelectric laminated plates［J］. Nonlinear Dynamics， 2012， 67（4）： 2441-2454. [百度学术]

El-Assal A M， Nada A A. A non-incremental finite element formulation of large deformation piezoceramic-laminated-plates［C］. 1st Joint International Conference on Multibody System Dynamics IMSDLappeenranta， Finland， 2010， [百度学术]

易灿明，余海东，王皓. 基于绝对节点坐标法的压电驱动复合结构动力学特性［J］. 上海交通大学学报， 2019， 53（6）： 665-672. [百度学术]

Yi Canming， Yu Haidong， Wang Hao. Dynamic behavior of coupled flexible plate structure with piezoelectric actuator via absolute nodal coordinate formulation ［J］. Journal of Shanghai Jiao Tong University， 2019， 53（6）： 665-672. [百度学术]