里兹法求解薄板振动方程

我的研究方向是对微机电压电水听器进行建模，查阅相关文献时往往会遇到这个偏微分方程：

\[ D\nabla^{4}w+\rho_{h}w_{tt}=0 \]

先不管参数 \(D,\rho_{h}\)，让我们关注方程本身，把这个式子当作单纯的“计算对象”。我没能在 Wolfram 的帮助文档里找到“一键式”的数值求解方法，而解析解只能处理固定边圆板、简支边圆板和简支边方板等特殊情况，如果是蜂巢形的水听器单元呢？

面对偏微分方程难解的问题，里兹法的思路大概是：先通过能量分析，把方程问题变成优化问题，再用基函数逼近真实解。在基函数设置时，就可以考虑薄板形状和边界条件的问题，从而得到一个比较通用的程序，利用 Wolfram 语言的内置函数，还可以把这个程序写得更简单。

最小作用量，最小物理知识

于我而言，“最小作用量”这套方法是在最小化物理知识的基础上，解决所面临的建模问题。

首先，还是要计算出薄板的动能和势能表达式：

\[ \begin{align} T&=\frac{\omega^{2}}{2}\rho_{h}a^{2}\iint_{\bar A}\bar w^{2}(x,y)\mathrm{d}x\mathrm{d}y\\ U&=\frac{1}{2}Da^{-2}\iint_{\bar A}(\nabla^{2}\bar w)^{2}\mathrm{d}x\mathrm{d}y \end{align} \]

其中 \(T\) 就是众所周知的 \(m v^2/2\) 这个形式，\(U\) 则需要参考基尔霍夫薄板理论的应力、应变表达式，并代入线弹性材料应变能的公式里。

补充两点：

积分区域经过尺度变换，所以可以将尺度系数 \(a\)（不妨理解为“半径”）提到积分之外。
为了简化过程，选取了 \(U\) 的最简形式，这仅对固定边界条件成立，如果是更一般的情况，应该添加一项：

\[ -2(1-\nu)\left[\frac{\partial^{2} w}{\partial x^{2}} \ \frac{\partial^{2} w}{\partial y^{2}}-\left(\frac{\partial^{2} \ w}{\partial x \partial y}\right)^{2}\right] \]

然后，用二者之差代表最小作用量：

\[ \Pi=T-U \]

“最小”并不意味着“最小”，只是取得极值，或者说对于自变函数 \(w(x,y)\) 来说，变分 \(\delta\Pi=0\)。所以，虽然里兹采用的能量泛函为 \(U-T\)，但这一点区别无关紧要。

另外，最小作用量的被积函数部分，就是拉格朗日量（或拉格朗日密度）。将其代入欧拉-拉格朗日方程后，同样可以得到前面给出的偏微分方程。对我来说，这要比从微元体出发的推导容易太多。

1 2	Needs["VariationalMethods`"] EulerEquations[1/2 \[Rho]h D[w[x, y, t], t]^2 - 1/2 dd (Laplacian[w[x, y, t], {x, y}])^2, w[x, y, t], {x, y, t}] // FullSimplify

用什么基函数来逼近真实解？

很多文章选取的基函数都很巧妙，既要满足边界条件，又要起到逼近真实解的作用。但我们打算得到比较通用的程序，就不能把基函数形式限制在一种特殊场景中。

所以，不妨将基函数写作两项相乘的形式¹：

边界方程

灰色部分为薄板区域，红色部分为薄板边界：

1	Graphics[{LightGray, #, Red, RegionBoundary@# }] & /@ {Disk[], Rectangle[{-1, -1}, {1, 1}], RegularPolygon[6]} // GraphicsRow

对于圆形薄板来说，其边界方程为：

\[ (x^2+y^2-1)^{i} \]

对于方形薄板来说，其边界方程为：

\[ (x-1)^{i}(x+1)^{j}(y-1)^{k}(y+1)^{l} \]

注意式子里的上标代表该边界的边界条件，\(0\) 对应于自由边界条件，\(1\) 对应于简支边界条件，\(2\) 对应于固定边界条件。

如果是正六边形板，这个边界方程就很难直接写出。可以先写一个函数，能输入各顶点坐标，输出边界方程：

(* 计算两点之间的线方程 *)
lineEquation[{p1_, p2_}, {x_, y_}] := 
  (x - p2[[1]]) (p1[[2]] - p2[[2]]) - (y - p2[[2]]) (p1[[1]] - p2[[1]])

(* 将点列表转换为成对的点 *)
pointsPair[points_List] := 
  Partition[Riffle[#, RotateLeft[#]], 2] & @ points

(* 计算多边形的边界 *)
polygonBorder[points_List, {x_, y_}] := 
  Times @@ (lineEquation[#, {x, y}] & /@ pointsPair @ points)

正六边形的六个顶点即为：

1	CirclePoints[6]

可计算得到边界方程：

1	polygonBorder[%, {x, y}]

逼近函数

使用勒让德多项式逼近函数：

多项式两两相乘，得到用于函数逼近的二元基函数：

以固定边界圆板为算例，其基函数形式为：

1	basicFunctions = Table[LegendreP[i - 1, x] LegendreP[j - 1, y], {i, 4}, {j, 4}] (x^2 + y^2 - 1)^2 // Flatten

终于变成矩阵问题

为方便计算，先把近似解写成基函数加权求和的形式：

\[ \hat w(x,y)=\sum\limits_{k=1}^{K}c_ {k}\psi_{k}(x,y) \]

代入动能 \(T\)、势能 \(U\) 的表达式中，积分部分即为系数矩阵的“二次型”：\(\boldsymbol{c}^\mathrm{T}\boldsymbol{Mc}\) 和 \(\boldsymbol{c}^\mathrm{T}\boldsymbol{Kc}\)，其中 \(\boldsymbol{M},\boldsymbol{K}\) 分别称为质量矩阵和刚度矩阵，矩阵元素的表达式为²：

\[ \begin{align} k_{ij}&=\iint_{\bar A}\nabla^{2}\psi_{i}\cdot\nabla^{2}\psi_{j} \ \mathrm{d}x\mathrm{d}y\\ m_{ij}&=\iint_{\bar A}\psi_{i}\cdot\psi_{j}\mathrm{d}x\mathrm{d}y \end{align} \]

继续前面的例子：Wolfram 语言允许我们在区域上进行积分，可计算出刚度矩阵和质量矩阵：

1
2

kmatrix = Table[Integrate[Laplacian[basicFunctions[[i]], {x, y}] Laplacian[basicFunctions[[j]], {x, y}], Element[{x, y}, Disk[]]], {i, 1, 16}, {j, 1, 16}];
mmatrix = Table[Integrate[basicFunctions[[i]] basicFunctions[[j]], Element[{x, y}, Disk[]]], {i, 1, 16}, {j, 1, 16}];

甚至可以发现，对于圆板来说，两个矩阵的所有元素竟然都有精确解：

1	kmatrix // MatrixForm

1	mmatrix // MatrixForm

用矩阵热图进行可视化，不难发现这两个矩阵都是对称矩阵，利用这一特性可以减少运算量。

1	MatrixPlot[#[[1]], PlotLabel -> #[[2]]] & /@ {{kmatrix, "刚度矩阵"}, {mmatrix, "质量矩阵"}} // GraphicsRow

通过选取合适的系数，取得作用量的极值，即作用量对向量 \(\boldsymbol{c}\) 求导为零：

\[ \frac{\partial\Pi}{\partial c}=0 \]

化简得到：

\[ (\boldsymbol{K}-\lambda^{4}\boldsymbol{M})\boldsymbol{c}=0 \]

1	{vals, coefs} = Eigensystem[{N@kmatrix, N@mmatrix}, -4];

把优化问题转化为质量矩阵和刚度矩阵的广义特征值求解问题，将特征向量 \(\boldsymbol{c}\) 与基函数向量 \(\boldsymbol{\psi}\) 点乘，即得近似解表达式。可将前四个模态的密度图绘制如下：

其中的密度图为模态振型，上面的数字为特征值 \(\lambda\)，与以下微分特征方程中的特征值 \(k\) 与半径 \(a\) 的乘积近似相等：

\[ \nabla^{4}w(r)=k^{4}w(r) \]

程序简化并总结为函数

参考上述分析过程，有以下几点优化空间：

第一，质量矩阵、刚度矩阵均为对称矩阵，可削减计算量：

1
2
3

Table[N@Integrate[Laplacian[basicFunctions[[i]], {x, y}] Laplacian[basicFunctions[[j]], {x, y}], Element[{x, y}, Disk[]]], {i, 1, 16}, {j, 1, 16}]; // RepeatedTiming

(*{1.52804, Null}*)

1
2
3

Table[N@Integrate[Laplacian[basicFunctions[[i]], {x, y}] Laplacian[basicFunctions[[j]], {x, y}], Element[{x, y}, Disk[]]], {i, 1, 16}, {j, 1, i}]; // RepeatedTiming

(*{0.813599, Null}*)

第二，计算积分时可并行计算，以减少程序耗时：

1
2
3

$ProcessorCount

(*4*)

1
2
3

ParallelTable[N@Integrate[Laplacian[basicFunctions[[i]], {x, y}] Laplacian[basicFunctions[[j]], {x, y}], Element[{x, y}, Disk[]]], {i, 1, 16}, {j, 1, i}]; // RepeatedTiming

(*{0.355318, Null}*)

第三，为编写更具一般性的程序，应使用数值积分。对于矩阵中存在的大量零元素，需要设定合适的准确度目标：

1
2
3

ParallelTable[NIntegrate[Laplacian[basicFunctions[[i]], {x, y}] Laplacian[basicFunctions[[j]], {x, y}], Element[{x, y}, Disk[]], AccuracyGoal -> 5], {i, 1, 16}, {j, 1, i}]; // RepeatedTiming

(*{1.47974, Null}*)

于是，将优化后的函数整合为一个程序包：

Clear[polygonBorder, lineEquation, pointsPair]

Get[FileNameJoin[{NotebookDirectory[], "RitzMethod.wl"}]]

?RitzSolver

对于上面的固定边界圆板，可以这样调用函数：

对于简支边界方板，需要多设置两个参数，一个是边界方程的指数，另一个是泊松比：

实际上，参考此算例的解析解表达式，可知泊松比并不影响结果：

1	SSSquare[{m_, n_}, {x_, y_}] := {Pi Sqrt[(m/2)^2 + (n/2)^2], Sin[Pi m (x + 1)/2] Sin[Pi n (y + 1)/2]}

绘图，并与里兹解作比较：

1	DensityPlot[#[[2]], {x, -1, 1}, {y, -1, 1}, PlotLabel -> N@#[[1]]] & /@Flatten[#, 1] &@Table[SSSquare[{i, j}, {x, y}], {i, 2}, {j, 2}] // GraphicsRow

最后是固定边界正六边形板，也就是我在文章最开始讲的“蜂巢形的水听器单元”：

第一个解看上去有些奇怪，可以绘制出立体图形，只是凹陷而已：

本文处理的方程是齐次的，所以不管解的幅值是多少，都能满足该方程，但为了方便后续处理，一般会对幅值作归一化：

为了降低理解难度，文中略去了许多细节。比如一般边界条件的刚度矩阵元素表达式，你可以在程序文件中找到。至于更详细的推导过程，可以参考下面列出的两篇文献，或者自行搜索“基尔霍夫薄板理论”的相关内容。

24/12/04

Zhao, J. (2023). Variable stiffness discrete Ritz method for free vibration analysis of plates in arbitrary geometries. Journal of Sound and Vibration, 553, 117662.↩︎
Leissa, A. W. (2005). The historical bases of the Rayleigh and Ritz methods. Journal of Sound and Vibration, 287(4–5), 961–978.↩︎

Joffoo's blog