方形平板振动克拉尼图形可视化计算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_{m+2,n,t}+u_{m,n-2,t}+u_{m,n+2,t})$
$=\frac{\rho h}{D}\frac{d{x}^4}{d{t}^2}(2u_{m,n,t-1}-u_{m,n,t-2})$

最终，整理为：
$$K\ast U=F$$
其中U为列向量，由前面网格位置$u_{m,n}$组成，是未知数。K为一个系数方阵，每一行为对应的$u_{m,n}$网格相应的系数。F为每个网格的受力，也是一个列向量。
最终根据K和F，求解线性方程组，就可以求解出网格位置U。

当然，只有方程还是不行的，微分方程另一个老大难就是边界条件（甚至很多问题边界条件比方程本身还难）。

因为前面${\nabla }^{4}u$用的13点差分格式，用到了上下左右向外方向上的2个距离的格点，所以为了能够顺利计算，就需要在平板的边缘外再补2圈虚拟节点。虚拟节点只与平板边缘的边界条件有关，所以时间项也都忽略为0。

首先是对称边界条件，这个很简单，就是沿着对称轴对称就行。下面以左右对称边界条件的格点方程为例：

然后是外侧自由无约束边界条件，第一层虚拟节点可以用内力M=0列方程得到：

对于平板四个角，只用Mx和My不能完全得到对应的方程，还需要引入Mxy=0

平板的角上第一层虚节点还差两个格点没有给出，再根据对称假设，假设沿对角线对称且Mx=My=0，得到：

第二层虚节点，我们用分布反力V=0来列出方程：

第二层虚节点在右上角还有3个节点没有方程涉及到，因为不涉及计算，所以全部强制赋值0即可。

因此，我们便得到了所有网格单元的方程，比如1/4平板网格数为6×6，加上两层虚节点，网格数变为10×10。算上前面列的方程，每个网格点都有一个独立的方程，因此我们便得到了100个线性方程。之后求解这个100个方程KU=F的方程组，就得到每个网格点的位置U。

这一段可能稍微啰嗦一些，但是如果想要动手编程计算的话，还是很难跳过的。

1.3 计算结果

之后是MATLAB计算仿真的结果，为完整展示振动图形，将1/4板做完整对称处理。虚拟节点没有展示。

这里只计算0.03s几个振动周期的平均结果，时间步长为1e-6s，空间网格点为16（对称完后为31个点）。

下图为650hz对应的结果：

黑色越深，代表振幅越小。

下图为792hz对应的结果：

下图为850hz对应的结果：

可以看到振动频率越高，图形越复杂。

这意味着频率越高，需要的空间网格量越大，而且频率越高，需要更小的时间步长。因此这种方法在高频率振动时，不是可行的展示办法。但在低频的时候，还是比较精准的。

代码如下，这里自己添加了一个阻尼项，因为如果不添加，平板振幅会越振越大。

clear
clc
close all

%平板信息
L=0.150;%边长300mm
h=0.001;%厚度10mm
mu=0.33;%泊松比0.33
rho=2.7e3;%板子的密度
E=70e9;%弹性模量
D=E*h^3/12/(1-mu^2);

%振动信息
Amp=0.01;
Freq=650;%频率越高，需要的网格越密。参考Freq=792,N=16。Freq=400,N=8。

%构建网格
dt=1e-6;%时间步长%1000hz的话，就是
N=16;%偶数%网格总数量为N*N，节点数目N+1*N+1
dx=L/N;%网格大小
x=dx*(0:N);

[X,Y]=meshgrid(x,x);%生成方形网格

%列方程组
[Sq,Bq]=Equation_Sq0(N,mu,-Amp);%初始值为-0.1的无外力分布
%已知Sq*U=Bq;%求U的分布，U即是当前平板上每一个点的位移
U0=Sq\Bq;
U1=U0;
U2=U1;

U_Out1=zeros(N+1,N+1);
U_Out2=zeros(N+5,N+5);
U_Save=zeros((N+1)^2,200);t_Save=1;
%计算动态方程
t_start=0;t_end=0.03;%起始时间和终止计算时间
jishu=1;%用于计数的一个变量
N_jishu=(t_end-t_start)/dt;
for t_k=t_start:dt:t_end
    %0点处的运动位置
    u0=Amp*cos(2*pi*Freq*t_k+pi);%以正弦方式运动
    L_Sq=N+5;%实际计算时网格的尺寸（包含外侧扩展的两层）
    if jishu==1%第一步计算，把矩阵Sq计算出来
        [Sq,Bq]=Equation_Sq0(N,mu,u0);%初始值为-0.1的无外力分布
        %补充运动项，把U1，U2代入，计算U3
        for k=1:L_Sq^2
            [r_k,c_k]=ind2sub([L_Sq,L_Sq],k);
            if (r_k>=3 && r_k<=L_Sq-2) && (c_k>=3 && c_k<=L_Sq-2)
                Sq(k,k)=20+rho*h*dx^4/D/dt^2;
                Sq(k,[k+1,k-1,k+L_Sq,k-L_Sq])=-8;
                Sq(k,[k+L_Sq+1,k+L_Sq-1,k-L_Sq+1,k-L_Sq-1])=2;
                Sq(k,[k+2,k-2,k-2*L_Sq,k+2*L_Sq])=1;
                Fd=-100*sign(U2(k)-U1(k))*(U2(k)-U1(k))^2/dt^2;%增加阻尼项
                Bq(k)=dx^4/D*(rho*h/dt^2*(2*U2(k)-U1(k))+Fd);
            end
        end
        %再重新定义一下中心约束点
        Indx_Center=sub2ind([L_Sq,L_Sq],3,3);
        Sq(Indx_Center,:)=0;Sq(Indx_Center,Indx_Center)=1;
        Bq(Indx_Center)=u0; %初始值为-0.1
        %Sq=sparse(Sq);转换为稀疏矩阵的形式，会让计算稍微快一些
    else
        %其余情况Sq矩阵不会变化，所以不用重复计算
        for k=1:L_Sq^2
            [r_k,c_k]=ind2sub([L_Sq,L_Sq],k);%这一块还可以优化，让速度更快一点
            if (r_k>=3 && r_k<=L_Sq-2) && (c_k>=3 && c_k<=L_Sq-2)
                Fd=-100*sign(U2(k)-U1(k))*(U2(k)-U1(k))^2/dt^2;%增加阻尼项
                Bq(k)=dx^4/D*(rho*h/dt^2*(2*U2(k)-U1(k))+Fd);
            end
        end
        Bq(1+2+2*L_Sq)=u0; %固定位置随u0变化
    end
    U3=Sq\Bq;
    %定义前两个时间步长下的位置信息
    U1=U2;
    U2=U3;
    %储存，用作输出用
    U_Out2(:)=U3(:);
    U_Out=U_Out2(3:end-2,3:end-2);


%     if mod(jishu,100)==1
%     figure(1)
%     clf
%     %pcolor(X,Y,U_Out)
%     mesh(X,Y,U_Out)
%     caxis([-0.2,0.2])
%     zlim([-0.2,0.2])
%     colorbar
%     pause(0.1)
%     disp(t_k)
%     end
    jishu=jishu+1;%时间步加一

    %记最后200个数据储存
    if jishu+50*200>=N_jishu
        if mod(jishu,50)==1
            U_Save(:,t_Save)=U_Out(:);
            t_Save=t_Save+1;
        end
    end
end

%取一些特征点，观察计算情况
figure()
hold on
for k=1:size(U_Save,1)
    plot(U_Save(k,:))
end
hold off

%填充对称图形
U_Out_A=U_Out;
U_Out_A(:)=max(U_Save,[],2)-min(U_Save,[],2);
U_Out_A2=[fliplr(U_Out_A(:,2:end)),U_Out_A];
U_Out_A3=[flipud(U_Out_A2(2:end,:));U_Out_A2];


%绘制
figure()
x3=dx*(-N:N);
[X3,Y3]=meshgrid(x3,x3);
mesh(X3,Y3,U_Out_A3)

%插值绘制（网格太稀疏了，绘图效果不好看，还是要插值一下）
xp=dx/2*(-2*N:2*N);
[Xp,Yp]=meshgrid(xp,xp);
F=griddedInterpolant(X3',Y3',U_Out_A3','spline');
U_Out_p=F(Xp,Yp);

figure()
sp=pcolor(Xp,Yp,U_Out_p);
mcp=[[linspace(0,0.5,16)',linspace(0,0.5,16)',linspace(0,0.5,16)'];
    [linspace(0.5,1,32)',linspace(0.5,1,32)',linspace(0.5,1,32)'];
    [1,1,1];[1,1,1]];
colormap(mcp)
sp.FaceColor = 'interp';


function [Sq,Bq]=Equation_Sq0(N,mu,u0)
%N网格数目
%mu泊松比
%u0初始平板位置
%外拓展两圈后平板网格的索引
L_Sq=N+5;
%角落边界点，都设置为0
Point_Corner0=[L_Sq,L_Sq;L_Sq-1,L_Sq;L_Sq,L_Sq-1];
%自由角垂直外边界，共2个
Point_CornerC=[L_Sq-1,L_Sq-2;L_Sq-2,L_Sq-1];
%第一层边界点（非对称）
Point_Out1=[(L_Sq-1)*ones(L_Sq-5,1),(3:L_Sq-3)';...
    (3:L_Sq-3)',(L_Sq-1)*ones(L_Sq-5,1)];
%对角线外边界，共1个
Point_Corner=[L_Sq-1,L_Sq-1];
%第二层边界点
Point_Out2=[L_Sq*ones(L_Sq-4,1),(3:L_Sq-2)';...
    (3:L_Sq-2)',L_Sq*ones(L_Sq-4,1)];
%第一层对称边界
Point_Mirror1=[2*ones(L_Sq-2,1),(3:L_Sq)';...
    (3:L_Sq)',2*ones(L_Sq-2,1)];
%第二层对称边界
Point_Mirror2=[1*ones(L_Sq-2,1),(3:L_Sq)';...
    (3:L_Sq)',1*ones(L_Sq-2,1)];
%左上角对称边界
Point_MirrorC=[1,1;1,2;2,1;2,2];

Sq=zeros(L_Sq^2);
Bq=zeros(L_Sq^2,1);

for k=1:L_Sq^2

    [r_k,c_k]=ind2sub([L_Sq,L_Sq],k);
    %四周边界点
    %四周角落边界点，都设置为0
    if IsRowInRowList(Point_Corner0,[r_k,c_k])
        Sq(k,k)=1;
        Bq(k)=0;
    end
    %自由角垂直外边界
    if IsRowInRowList(Point_CornerC,[r_k,c_k])
        if r_k==2
            Sq(k,k:k+2)=[1,-2,1];
        elseif r_k==L_Sq-1
            Sq(k,k-2:k)=[1,-2,1];
        elseif c_k==2
            Sq(k,[k,k+L_Sq,k+2*L_Sq])=[1,-2,1];
        elseif c_k==L_Sq-1
            Sq(k,[k-2*L_Sq,k-L_Sq,k])=[1,-2,1];
        end
        Bq(k)=0;
        %计算第一层边界点
    elseif IsRowInRowList(Point_Out1,[r_k,c_k])
        if r_k==2 %My=0
            Sq(k,[k+1-L_Sq,k+1,k+1+L_Sq])=[-mu,2+2*mu,-mu];
            Sq(k,k)=-1;Sq(k,k+2)=-1;
        elseif r_k==L_Sq-1 %My=0
            Sq(k,[k-1-L_Sq,k-1,k-1+L_Sq])=[-mu,2+2*mu,-mu];
            Sq(k,k)=-1;Sq(k,k-2)=-1;
        elseif c_k==2 %Mx=0
            Sq(k,[k,k+L_Sq,k+2*L_Sq])=[-1,2+2*mu,-1];
            Sq(k,k+L_Sq-1)=-mu;Sq(k,k+L_Sq+1)=-mu;
        elseif c_k==L_Sq-1 %Mx=0
            Sq(k,[k,k-L_Sq,k-2*L_Sq])=[-1,2+2*mu,-1];
            Sq(k,k-L_Sq-1)=-mu;Sq(k,k-L_Sq+1)=-mu;
        end
        Bq(k)=0;
        %自由角对角线外边界，每个角1个
    elseif IsRowInRowList(Point_Corner,[r_k,c_k])
        if r_k==2 && c_k==2
            Sq(k,[k,k+2*L_Sq+2])=[1,1];
            Sq(k,[k+2,k+2*L_Sq])=[-1,-1];
        elseif r_k==L_Sq-1 && c_k==2
            Sq(k,[k,k+2*L_Sq-2])=[1,1];
            Sq(k,[k-2,k+2*L_Sq])=[-1,-1];
        elseif r_k==2 && c_k==L_Sq-1
            Sq(k,[k,k-2*L_Sq+2])=[1,1];
            Sq(k,[k+2,k-2*L_Sq])=[-1,-1];
        elseif r_k==L_Sq-1 && c_k==L_Sq-1
            Sq(k,[k,k-2*L_Sq-2])=[1,1];
            Sq(k,[k-2,k-2*L_Sq])=[-1,-1];
        end
        Bq(k)=0;
        %4计算第二层边界点
    elseif IsRowInRowList(Point_Out2,[r_k,c_k])
        if r_k==1
            Sq(k,k)=1;
            Sq(k,[k+1-L_Sq,k+1,k+1+L_Sq])=[2-mu,2*mu-6,2-mu];
            Sq(k,[k+3-L_Sq,k+3,k+3+L_Sq])=[mu-2,-2*mu+6,mu-2];
            Sq(k,k+4)=-1;
        elseif r_k==L_Sq
            Sq(k,k)=1;
            Sq(k,[k-1-L_Sq,k-1,k-1+L_Sq])=[2-mu,2*mu-6,2-mu];
            Sq(k,[k-3-L_Sq,k-3,k-3+L_Sq])=[mu-2,-2*mu+6,mu-2];
            Sq(k,k-4)=-1;
        elseif c_k==1
            Sq(k,k)=1;
            Sq(k,[k+L_Sq-1,k+L_Sq,k+L_Sq+1])=[2-mu,2*mu-6,2-mu];
            Sq(k,[k+3*L_Sq-1,k+3*L_Sq,k+3*L_Sq+1])=[mu-2,-2*mu+6,mu-2];
            Sq(k,k+4*L_Sq)=-1;
        elseif c_k==L_Sq
            Sq(k,k)=1;
            Sq(k,[k-L_Sq-1,k-L_Sq,k-L_Sq+1])=[2-mu,2*mu-6,2-mu];
            Sq(k,[k-3*L_Sq-1,k-3*L_Sq,k-3*L_Sq+1])=[mu-2,-2*mu+6,mu-2];
            Sq(k,k-4*L_Sq)=-1;
        end
        Bq(k)=0;
        %计算除边界点外的正常平板上的点
    elseif (r_k>=3 && r_k<=L_Sq-2) && (c_k>=3 && c_k<=L_Sq-2)
        Sq(k,k)=20;
        Sq(k,[k+1,k-1,k+L_Sq,k-L_Sq])=-8;
        Sq(k,[k+L_Sq+1,k+L_Sq-1,k-L_Sq+1,k-L_Sq-1])=2;
        Sq(k,[k+2,k-2,k-2*L_Sq,k+2*L_Sq])=1;
        Bq(k)=0;%dx^4/D*(rho*h/dt^2*(2*0-0));
        %计算对称边界的外插点
    elseif IsRowInRowList(Point_Mirror1,[r_k,c_k])
        if r_k==2
            Sq(k,k)=1;Sq(k,k+2)=-1;
        elseif c_k==2
            Sq(k,k)=1;Sq(k,k+2*L_Sq)=-1;
        end
        Bq(k)=0;
    elseif IsRowInRowList(Point_Mirror2,[r_k,c_k])
        if r_k==1
            Sq(k,k)=1;Sq(k,k+4)=-1;
        elseif c_k==1
            Sq(k,k)=1;Sq(k,k+4*L_Sq)=-1;
        end
        Bq(k)=0;
    elseif IsRowInRowList(Point_MirrorC,[r_k,c_k])
        if r_k==1 && c_k==1
            Sq(k,k)=1;Sq(k,k+4+4*L_Sq)=-1;
        elseif r_k==1 && c_k==2
            Sq(k,k)=1;Sq(k,k+4+2*L_Sq)=-1;
        elseif r_k==2 && c_k==1
            Sq(k,k)=1;Sq(k,k+2+4*L_Sq)=-1;
        elseif r_k==2 && c_k==2
            Sq(k,k)=1;Sq(k,k+2+2*L_Sq)=-1;
        end
        Bq(k)=0;
    end
end
% Sq([32,41,42],:)=[];Bq([32,41,42])=[];
% rank(Sq)%检查是否满秩

%补充两个边界约束（中心点已知，）
%初始已知中心点坐标
Sq(1+2+2*L_Sq,:)=0;Sq(1+2+2*L_Sq,1+2+2*L_Sq)=1;
Bq(1+2+2*L_Sq)=u0; %初始值为-1
end

function TF=IsRowInRowList(List,Point)
TF1=(List(:,1)==Point(1));
TF= any(List(TF1,2)==Point(2));
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294

方程建立有几点小的心得体会：
1是要寻找一种让矩阵Sq保持不变的方程形式，这种矩阵在力学求解器或者其它方程求解器中很常见，比如刚度矩阵、气动矩阵等等。
2这里为了方便展示和检查，在代码中把Sq=sparse(Sq);这个注释掉了。实际上Sq矩阵是一个稀疏矩阵，在matlab中转换成稀疏矩阵格式，可以显著的减少内存加速计算。
3方程的建立初期还是比较痛苦的，需要检查是否满秩，约束是否刚好够。比如全模型的方程很容易列出少约束的情况，这时就需要找到哪些方程是重复无用的方程（比如中间固定点旁边的几个点，就可以用对称或反对称边界条件替换原物理方程）。

2 简单的波动解

参考文献：
[1]关于Chladni图案
https://zhuanlan.zhihu.com/p/108448193
[2]Creating Digital Chladni Patterns
https://thelig.ht/chladni/

这里我们做一个大胆的假设，平板上每一个点都在做周期运动sin(wt+k)，而且相位差不多，只是振幅大小不同。

那么，我们不需要知道上一章所列的平板振动方程，也可以大概列出平板方程的通解。

设平板上每一点的振幅为u(x,y,t)，可以表示为：
$$u(x,y,t)=w(x,y)\ast \sin (wt+k)$$
w(x,y)是平板上每个点的振幅。当然后面那个w(频率)和k(相位)我们不关心，因为后面用不到。

如果我们知道每个点对应的振幅，那么振幅等于0的位置就是Chladni图案。所以下一步我们就计算振幅w(x,y)的通解。

现在建立平板模型，设平板的范围为[-0.5,0.5]，中心点为(0,0)。平板的边缘自由振动。

假定w(x,y)是由若干个余弦波叠加而成的。而平板要想形成共振，需要承载较为完整的驻波。所以X方向大致满足条件的解可以有如下的表示方式：
$$X(x)=\cos (n\pi (x-0.5))$$
其中n为方形板上x方向波峰和波谷总和的数量。这里假定平板的边缘在余弦函数的峰或谷上。Y方向同理。

振动时，具体w(x,y)的振幅由当地的X(x)和Y(y)相乘决定，即：
$$\cos (n\pi (x-0.5))\ast \cos (m\pi (y-0.5))$$
n和m不一定相同。因此考虑x和y的对称性，m和n对换的情况也会同时叠加出现。最终振幅w(x,y)的方程为
$$\begin{gathered} w(x,y)=\cos (n\pi (x-0.5))\ast \cos (m\pi (y-0.5)) \\ \pm \cos (m\pi (x-0.5))\ast \cos (n\pi (y-0.5)) \end{gathered}$$
上面式子中正号代表对称图案情况，负号代表反对称图案情况。

最终不同m和n取值下，生成的Chladni图案如下：

根据某些文献中给出的Chladni’s Law，可以根据生成的图案来大概反推频率：
$$f\sim (m+2n{)}^2$$
也就是频率f和(m+2n)的平方成正比。

上图的matlab代码如下：

clear
clc
close all
%绘制所有可能波节的输出结果
N=9;
figure(1)
for m=1:N
    for n=1:N
        %ax=subplot(N,N,N*(m-1)+n);
        dp=1/N;
        ax=subplot('Position',[(m-1)*dp 1-(n-0)*dp 0.9*dp 0.9*dp]);
        
        fun1=@(x,y) cos(n*pi*(x-0.5)).*cos(m*pi*(y-0.5));
        fun2=@(x,y) cos(m*pi*(x-0.5)).*cos(n*pi*(y-0.5));
        fun=@(x,y) fun1(x,y)+sign(m-n)*fun2(x,y);
        fimplicit(fun,[-0.5,0.5])
        ax.XTick=[];
        ax.YTick=[];
    end
end
set(gcf,'position',[100,100,800,600])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

3 理论求解

参考文献：
[1]方奕忠, 王钢, 沈韩,等. 方形薄板二维驻波的研究[J]. 物理实验, 2014
[2]Neville H. Fletcher , Thomas D. Rossing .The Physics of Musical Instruments [M] ISBN-13:978-0-387-94151-6
[3]严琪琪, 陈彦, 胡湘. 自由边界条件下方形平板受迫振动模式的探究[J]. 大学物理, 2018

当然，上面的方程还是太过于简化，下面给出了一种能够较好预测出结果的一种理论解方法。下面方程参考参考文献[1]方形薄板二维驻波的研究。（其实[3]自由边界条件下方形平板受迫振动模式的探究 那篇文章的效果看上去更好，但是我没太看懂 (╯°Д°)╯︵┻━┻，先用简单的顶上）。

首先还是需要求解下面这个方程：
$${\nabla }^{4}u+\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}=0$$
依然假设方程的解为：
$$u(x,y,t)=w(x,y)\ast \sin (\omega t+k)$$
代入回方程，可以得到关于振幅w(x,y)的方程：
$${\nabla }^{4}w-{\gamma }^{2}w=0$$
在x=0和x=1的边界条件为：
$$\begin{gathered} \frac{{\partial }^{2}w}{\partial x^2}+\mu \frac{{\partial }^{2}w}{\partial y^2}=0 \\ \frac{\partial }{\partial x}(\frac{{\partial }^{2}w}{\partial x^2}+(2-\mu )\frac{{\partial }^{2}w}{\partial y^2})=0 \end{gathered}$$
在y=0和y=1的边界条件同理，为：
$$\begin{gathered} \frac{{\partial }^{2}w}{\partial y^2}+\mu \frac{{\partial }^{2}w}{\partial x^2}=0 \\ \frac{\partial }{\partial y}(\frac{{\partial }^{2}w}{\partial y^2}+(2-\mu )\frac{{\partial }^{2}w}{\partial x^2})=0 \end{gathered}$$

参考文献[1]的求解过程中也用到了一些简化，解的形式和上一章节比较像。由于本人没有具体实验，所以如何利用频率反推波节这一块，公式精度具体有多少，需要怎样修正，都是未知数。本文也不是论文，就不再做实验验证了。
理论上看该方法在较高频率预测效果应该会比较好。

clear
clc
close all

%绘制理论求解图形
%平板信息
L=0.5;
h=0.001;%厚度1mm
mu=0.33;
rho=2.7e3;%板子的密度
E=70e9;%弹性模量
D=E*h^3/12/(1-mu^2);
c=sqrt(E/rho)*h/sqrt(3)/sqrt(1-mu^2);%计算一个常数
%振动频率
f=1560;
%寻找可能的波节m和n
n_t=round(sqrt(f/pi/c*L^2)-0.5);
n_t1=n_t+(-8:8);
[m_t2,n_t2]=meshgrid(n_t1,n_t1);

MN2Freq=pi/2/L^2*c*((n_t2+0.5).^2+(m_t2+0.5).^2);%已知mn，反推频率的一个公式
[~,Ind_nt2]=sort(abs(MN2Freq(:)-f));
%得到接近当前频率的几组m和n。
[x,y]=meshgrid(0:0.01:L,0:0.01:L);
for k=1:16
    %计算方程的解
    mm=m_t2(Ind_nt2(k));
    nn=n_t2(Ind_nt2(k));
	%得到可能的mm值和nn值，然后计算出bn和bm。
    bn=nn+0.5;
    Ps1_n1=sqrt(2)/2*(exp(-pi*bn*x/L)-(-1)^nn*exp(pi*bn*(x-L)/L));
    Ps1_n2=-sin(pi*bn/L*x-pi/4);
    Ps1_n=Ps1_n1+Ps1_n2;

    bm=mm+0.5;
    Ps1_m1=sqrt(2)/2*(exp(-pi*bm*y/L)-(-1)^mm*exp(pi*bm*(y-L)/L));
    Ps1_m2=-sin(pi*bm/L*y-pi/4);
    Ps1_m=Ps1_m1+Ps1_m2;
	
    Ps1_mn=Ps1_m.*Ps1_n;
    
	%参照上一章，做出x和y交换后的另一组解
    Ps2_n1=sqrt(2)/2*(exp(-pi*bn*y/L)-(-1)^nn*exp(pi*bn*(y-L)/L));
    Ps2_n2=-sin(pi*bn/L*y-pi/4);
    Ps2_n=Ps2_n1+Ps2_n2;

    bm=mm+0.5;
    Ps2_m1=sqrt(2)/2*(exp(-pi*bm*x/L)-(-1)^mm*exp(pi*bm*(x-L)/L));
    Ps2_m2=-sin(pi*bm/L*x-pi/4);
    Ps2_m=Ps2_m1+Ps2_m2;
    Ps2_mn=Ps2_m.*Ps2_n;
    
	%绘图
    figure(1)
    subplot(4,4,k)
    contour(Ps1_mn+Ps2_mn,1)
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

当频率等于1560hz时，比较可能的几组解如下：

当频率等于2340hz时，比较可能的几组解如下：