粘性射流线性稳定性色散方程求解

发表于 更新于 分类于

色散方程

射流线性稳定性中色散方程（Dispersion relation）是描述射流受到扰动时，时间扰动增长率和空间扰动增长率的关系。假设粘性射流受到轴向波数为 \(k\) 、周向阶数为 \(n\)、时间增长率为 \(\omega\) 的扰动，粘性相关无量纲数为 \(Oh\)，则色散方程的无量纲形式 \(D(Oh;k,n,\omega)=0\) 如下所示 \[ \begin{equation}\label{eq1} \begin{vmatrix} D_{11} & D_{12} & D_{13} \\ 2\mathrm{i}kI_n^\prime(k) & -\frac{l^2}{k^2}I_n^\prime(l)-I_{n+1}(l) & \frac{l^2}{k^2}I_n^\prime(l)+I_{n-1}(l) \\ 2\mathrm{i}n[kI_n^\prime(k)-I_n(k)] & lI_{n+2}(l) & lI_{n-2}(l) \end{vmatrix} = 0, \end{equation} \]

其中 \[ \begin{equation} \begin{aligned} D_{11} &= \omega\left[ 2Ohk^{2}I_{n}^{ \prime \prime }(k)+ \omega I_{n}(k) \right] -(1-n^{2}-k^{2})kI_{n}^{ \prime }(k), \\ D_{12} &= 2\mathrm{i}\omega OhlI_{n+1}^{ \prime }(l) -(1-n^{2}-k^{2})\mathrm{i}I_{n+1}(l), \\ D_{13} &= -2\mathrm{i}\omega OhlI_{n-1}^{ \prime }(l) + (1-n^{2}-k^{2})\mathrm{i}I_{n-1}(l), \end{aligned} \end{equation} \]

以及 \[ l=\sqrt{k^2+\frac{\omega}{Oh}}. \] 以上方程出现第一类修正 Bessel 函数 \(I_n (x)\)，它是一个无穷级数，求解方程 \(\eqref{eq1}\) 时涉及求解 \(I_n (x)\) 的反函数。

符号求解

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
clear; close all; clc;
% 定义初始变量
Oh0 = 1;
n0 = 0;
k0 = 1e-2:1e-2:1;

% 定义符号量
syms Oh k omega n l
assume(in(real(omega),'positive'));

% 行列式
D11 = omega*(2*Oh*k^2*diff(besseli(n,k),k,2) + omega*besseli(n,k)) ...
    - (1 - n^2 - k^2)*k*diff(besseli(n,k),k);
D12 = 2i*omega*Oh*l*diff(besseli(n+1,l),l)...
    - (1 - n^2 - k^2)*1i*besseli(n+1,l);
D13 = -2i*omega*Oh*l*diff(besseli(n-1,l),l)...
    + (1 - n^2 - k^2)*1i*besseli(n-1,l);
D21 = 2i*k*diff(besseli(n,k),k);
D22 = -l^2/k^2*diff(besseli(n,l),l) - besseli(n+1,l);
D23 = l^2/k^2*diff(besseli(n,l),l) + besseli(n-1,l);
D31 = 2i*n*(k*diff(besseli(n,k),k) - besseli(n,k));
D32 = l*besseli(n+2,l);
D33 = l*besseli(n-2,l);
D = [D11, D12, D13;...
    D21, D22, D23;...
    D31, D32, D33];
eqn = det(D);
eqn = subs(eqn, l, sqrt(k^2+omega/Oh));
eqn = subs(eqn, Oh, Oh0);
eqn = subs(eqn, n, n0);

% 求解
omega0 = zeros(1,length(k0));
for i=1:length(k0)
    eqn0 = subs(eqn, k, k0(i));
    omega0(i) = vpasolve(eqn0, omega, 1);
end

% 保存数据
data = [k0', omega0'];
save(['Oh=', num2str(Oh0), '_n0=', num2str(n0), '.txt'], 'data', '-ascii');