方形平板振动克拉尼图形可视化计算MATLAB程序（Chladni Patterns）

惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

0前言

克拉尼图形（Chladni Patterns）是在1787年，由克拉尼首先发现并命名的。他将一个金属薄板中央固定，然后把细沙撒在金属板上，用小提琴摩擦边缘，板子上的细沙便会形成各种不同的图案。

相关的实验非常多，很多科技馆或者大学实验课或者网上视频大家也都接触过。具体原理也不难，就是每种频率下方板的振动模态是固定的，有的地方振动幅度大，有的地方振动幅度小，沙子在不断上下颠的过程中肯定会逐渐往振动小的地方汇集，于是就形成了克拉尼图形。

所以要想数值求解克拉尼图形，最重要的还是要求解出相应的振动方程才行。

本博客将求解方法分为3种，分别对应网上或者论文中常见的求解方式。相应文献这次就不放到前言了，放到后面每一章的开头。

下面的代码涉及到简单的数值分析基础，有些基础知识可能不会详细涉及。

1 数值时域求解

最近解微分方程比较上头，所以这里也顺手写了一个，效率比较低。

参考文献：
[1] 弹性力学（下册）(徐芝纶版) （第十四章：用差分法及变分法解薄板的小挠度弯曲问题）
[2] Abdeljaber O , Rjoub Y A . Free and forced vibration of rectangular plates using the finite difference method[C]// Green Building, Materials and Civil Engineering. 2014.

1.1 方程建立

根据常见的克拉尼图形实验来看，这属于弹性力学里的薄板小变形问题。

薄板静变形方程为：
$${\nabla }^{4}u=\frac{{\partial }^{4}u}{\partial x^4}+2\frac{{\partial }^{4}u}{\partial x^2\partial y^2}+\frac{{\partial }^{4}u}{\partial y^4}=0$$
其中u为薄板变形位移。当然，本文只考虑垂直于平板表面的位移，其余位移忽略。
再考虑运动，还需要引入时间t和加速度对应的量，方程变为：
$${\nabla }^{4}u+\frac{\rho h}{D}\frac{{\partial }^{2}u}{\partial t^2}=0$$
其中$D=E{h}^3/12(1-\mu {)}^2$，E为平板材料的弹性模量，μ为材料的泊松比，ρ为材料的密度，这都是材料属性。h为薄板的厚度。

平板位置为在xy平面的[-L,L]区间上，四周为自由边界条件。
平板中心(0,0)点的位置为固定边界条件，中间固定点位移为已知的余弦函数振动。

为简化计算量，后面1.2节将平板沿x轴和y轴裁剪为1/4对称的小方形，中间设置为对称边界条件。

这种1/4模型简化有一定的问题，比如当平板做非对称振动时就无法模拟，而且Chladni图形的实验中也观察到这种图形。这个在后面第2章和第3章会涉及。

1.2 数值差分方程建立

首先${\nabla }^{4}u$，我们直接采用下面的差分格式：

这里假设dx和dy方向的网格尺寸相同。当然，上面的表达式也并非唯一，这种13点差分格式的优点为结构比较紧凑，用到的点比较少，其余格式也可以自行泰勒展开凑一凑。

时间项由于只有二阶，所以更简单一些：
$$\frac{{\partial }^{2}u}{\partial t^2}\to \frac{u_{t-2}-2u_{t-1}+u_t}{d{t}^2}$$
之后将上面所有单元的对应的差分方程列出后，将显示时间格式项加入方程中，然后联立这些方程组得到当前时刻的空间u。每个单元格点的方程u表示如下：

$(20+\frac{\rho h}{D}\frac{d{x}^4}{d{t}^2})u_{m,n,t}-8(u_{m-1,n,t}+u_{m+1,n,t}+u_{m,n-1,t}+u_{m,n+1,t})$
$+2(u_{m-1,n-1,t}+u_{m-1,n+1,t}+u_{m+1,n+1,t}+u_{m+1,n-1,t})$
$+(u_{m-2,n,t}+u_{}$