设飞行器是细长体, 长l m, 轴向单位长度上的质量密度为ρ (x) , E (x) J (x) 是横坐标为x的截面上的抗弯刚度, 发动机推力记为P。设发动机推力与纵轴重合并与飞行器固联。飞行器发生横向振动或其它类型运动时发动机喷管推力中线与纵轴端点相切。因此, 发动机推力的横向分量将极大地影响振动参数。

空气的轴向阻力假定按某一规律C (x) 分布于飞行器全长上。换言之, C (x) 表示横坐标为x的点附近单位长度上的轴向空气阻力。在飞行中某一t时刻坐标为x的截面上的轴向内力记为N (x, t) , 剪力记为Q (x, t) , 弯曲力矩用M (x, t) 表示。今观察垂直于纵轴截下的微元上的力矩在t时刻的平衡方程式 (见图1) 。

用u (x, t) 表示横截面坐标为x处在t时刻的横向位移, 用v表示飞行器的速度。注意到轴向力是由推力P在轴向投影Pcos θ ( θ是弯曲轴线的尾端点上的切线与平衡态轴线ox的夹角) 、分布阻力C (x, t) Δx和惯性力-ρ (x, t) $\frac{\text{d}v}{\text{d}t}$合成的。由图1.1可写出

${\rm M}+\text{Δ}{\rm M}-{\rm M}-Q\text{Δ}x+{\rm N}\text{Δ}u=0。$

当Δx→0时

$\frac{\partial {\rm M}(x,t)}{\partial x}=Q(x,t)-{\rm N}(x,t)\frac{\partial u(x,t)}{\partial x}。(1.1)$

在轴线发生弯曲时, 推力中线将发生偏转。在小振动情况下, 此时推力在ox方向 (平衡态) 的投影Pcos θ≈P, 而在垂直方向的投影为${\rm P}\sin \theta \approx {\rm P}\tan \theta \approx {\rm P}\frac{\partial u}{\partial x}$。此时, 用G₀表示飞行器总质量

$G{}_0(t)={\displaystyle {\int }_0^l\rho }(x,t)\text{d}x。(1.2)$

依设, 飞行器受到的总轴向阻力为

$D{}_0(t)={\displaystyle {\int }_0^{l}C}(x,t)\text{d}x。(1.3)$

飞行器在t时刻的轴向加速度为

$\frac{\text{d}v}{\text{d}t}=\frac{{\rm P}-D{}_0}{G{}_0}g‚(1.4)$

作用到横坐标为x处的轴向力为

$\begin{array}{l}{\rm N}(x,t)={\rm P}(t)-\\ \frac{{\rm P}(t)-D{}_0(t)}{G{}_0(t)}G(x,t)-D(x,t),(1.5)\end{array}$

式中$\begin{array}{l}G(x,t)=g{\displaystyle {\int }_0^x\rho }(s,t)\text{d}s‚\\ D(x,t)={\displaystyle {\int }_0^{x}C}(s,t)\text{d}s\end{array}(1.6)$

如果把 (1.5) 式改写成

$\begin{array}{l}{\rm N}(x,t)={\rm P}(t)\left(1-\frac{G(x,t)}{G{}_0}\right)+\\ D{}_0(t)\left(\frac{G(x,t)}{G{}_0(t)}-\frac{D(x,t)}{D{}_0(t)}\right)‚(1.7)\end{array}$

容易看到N (x, t) 满足边界条件:

N (o, t) =P (t) , N (l, t) =0。

现将N (x, t) 的表达式 (1.7) 代入力矩平衡方程式 (1.1) :

$\begin{array}{l}\frac{\partial {\rm M}(x,t)}{\partial x}=Q(x,t)-{\rm P}(t)\left(1-\frac{G(x,t)}{G{}_0(t)}\right)\cdot \\ \frac{\partial u}{\partial x}-D{}_0(t)\left(\frac{G(x,t)}{G{}_0(t)}-\frac{D(x,t)}{D{}_0(t)}\right)\frac{\partial u}{\partial x}‚(1.8)\end{array}$

上式再对x求二次微商后得

$\begin{array}{l}\frac{\partial {}^2{\rm M}(x,t)}{\partial x{}^2}=\frac{\partial Q(x,t)}{\partial x}-{\rm P}(t)\frac{\partial }{\partial x}\left(1-\frac{G(x,t)}{G{}_0(t)}\right)\cdot \\ \frac{\partial u}{\partial x}-D{}_0(t)\frac{\partial }{\partial x}\left(\frac{G(x,t)}{G{}_0(t)}-\frac{D(x,t)}{D{}_0(t)}\right)\frac{\partial u}{\partial x}。(1.9)\end{array}$

显然, 剪力Q (x, t) 的变化可写成外分布力和侧向惯性力之和:

$\frac{\partial Q(x,t)}{\partial x}=q(x,t)-\rho (x,t)\frac{\partial {}^{2}u}{\partial t{}^2}。(1.10)$

又依梁的小弯曲理论知M (x, t) =EJ (x, t) ·$\frac{\partial {}^{2}u}{\partial x{}^2}$, 将M和Q的值代入 (1.9) 式后便得到弹性细长体在推力和阻力共同作用下的横向运动方程式 ^[1]:

$\begin{array}{l}\rho (x,t)\frac{\partial {}^{2}u}{\partial t{}^2}+\frac{\partial {}^2}{\partial x{}^2}E(x)J(x,t)\frac{\partial {}^{2}u}{\partial x{}^2}+\\ \frac{\partial }{\partial x}[{\rm P}(t)\left(1-\frac{G(x,t)}{G{}_0(t)}\right)+D{}_0(t)(\frac{G(x,t)}{G{}_0(t)}-\\ \frac{D(x,t)}{D{}_0(t)})]\frac{\partial u}{\partial x}=q(x,t),(1.11)\end{array}$

或者简写为

$\begin{array}{l}\rho (x,t)\frac{\partial {}^{2}u}{\partial t{}^2}+\frac{\partial {}^2}{\partial x{}^2}E(x,t)J(x,t)\frac{\partial {}^{2}u}{\partial x{}^2}+\\ \frac{\partial }{\partial x}{\rm N}(x,t)\frac{\partial u}{\partial x}=q(x,t),(1.12)\end{array}$

式内N (x, t) 由 (1.5) 或 (1.7) 式确定。

在具体问题的讨论中可以区分下列三种不同的情况。

${\rm N}(x,t)={\rm P}(t)\left(1-\frac{G(x,t)}{G{}_0(t)}\right)$。

${\rm N}(x,t)=D{}_0(t)\left(\frac{G(x,t)}{G{}_0(t)}-\frac{D(x,t)}{D{}_0(t)}\right)。$

$D{}_0(t)={\rm P}(t)‚{\rm N}(x,t)={\rm P}(t)-D(x,t)。$

总之, 由于飞行状态不同, 选取不同的轴向力函数N (x, t) , 基本方程式 (1.11) 可以描述飞行器在各种状态下的横向运动。

现在再来讨论基本方程式 (1.11) 的边界条件。首先, 无论是主动段或被动段, 在x=0和x=1的两个端点上力矩都等于零, 所以

$EJ(x,t)\frac{\partial {}^{2}u}{\partial x{}^2}|{}_{x=0}=EJ(x,t)\frac{\partial {}^{2}u}{\partial x{}^2}|{}_{x=l}=0。(1.13)$

至于两端的剪力情况与自由梁有根本不同。如果坐标原点置于发动机喷管后端面中心, 由于发动机随飞行器的横向振动而摆动将产生外加剪力${\rm P}(t)\frac{\partial u}{\partial x}|{}_{x=0}$, 故第二组边界条件可以写为:

$\begin{array}{l}\left[\frac{\partial }{\partial x}EJ(x,t)\frac{\partial {}^{2}u}{\partial x{}^2}+{\rm P}(t)\frac{\partial u}{\partial x}\right]{}_{x=0}=0\\ \frac{\partial }{\partial x}EJ(x,t)\frac{\partial {}^{2}u}{\partial x{}^2}|{}_{x=l}=0‚\end{array}(1.14)$

这样, (1.13) 和 (1.14) 联合起来构成所讨论问题的边界条件。注意到, 与自由梁不同的是 (1.14) 中的第1式, 而其他三个边界条件则与两端自由的梁是相同的。当阻力存在时, 锥型前端头部所产生的阻力是由激波后压力升高造成的, 故前端阻力实际上是分布在头部附近的锥面上, 而包括在阻力函数C (x, t) 中。在本节所作基本假定的条件下, 边界条件 (1.13) 和 (1.14) 可望相当准确地成立 ^[2]。

为了讨论方便, 常把 (1.12) 式改写成一阶方程组。设ρ (x, t) 总大于某一正常数, 定义算子

$\begin{array}{l}A(t)=\frac{1}{\rho (x,t)}[\frac{\partial {}^2}{\partial x{}^2}EJ(x,t)\frac{\partial {}^2}{\partial x{}^2}+\\ \frac{\partial }{\partial x}{\rm N}(x,t)\frac{\partial }{\partial x}]‚(1.15)\end{array}$

和引进符号$\frac{\partial u}{\partial t}=v$, 则 (1.12) 式可改写为:

$\begin{array}{l}\frac{\partial u}{\partial t}=v\\ \frac{\partial v}{\partial t}=-Au+q(x,t)\end{array}$

或者写成向量形式

$\frac{\partial }{\partial t}\left(\begin{array}{l}u\\ v\end{array}\right)=\left(\begin{array}{cc}0& 1\\ -A& 0\end{array}\right)\left(\begin{array}{l}u\\ v\end{array}\right)+\left(\begin{array}{l}0\\ 1\end{array}\right)q(x,t)。(1.16)$

综合第1节内的讨论, 基本方程式 (1.11) 和边界条件 (1.13) 、 (1.14) 一起构成飞行器的横向运动方程组, 当然也包括刚体横向运动部分。本节将讨论关于横向振动的振型和频率的理论。从 (1.11) 式中可看到, 飞行器的各种参数是随时间t而变化的。当主要研究振动问题, 而不是研究刚体运动时, 因为飞行器参数的变化比振动参数的变化慢, 我们有理由用固化系数的概念, 把诸系数固化。例如分别研究主动段、被动段、巡航段的横向振动, 可以认为在某一较长时间间隔中, 系数EJ (x, t) 、N (x, t) 等仅是横坐标x的函数。实践证明, 这样做至少对振动问题的研究是合理的。

现应用分离系数法来讨论本征值问题。设在 (1.11) 中令q=0, 存在有如y (x, t) =y (x) ·exp (λt) 形式的解, 其中y (x) 满足边界条件 (1.13) 和 (1.14) 。代入 (1.11) 、 (1.13) 及 (1.14) 并消去非零因子exp (λt) 后, 再采取固化系数的方法, 认为诸系数在某一时间间隔内不依赖于时间t, 就得到常微分方程边值本征值问题:

$\begin{array}{l}\lambda {}^{2}y(x)+\frac{1}{\rho (x)}\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}+\\ \frac{1}{\rho (x)}\frac{\text{d}}{\text{d}x}{\rm N}(x)\frac{\text{d}y}{\text{d}x}=0‚(2.1)\\ EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}|{}_{x=0}=EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}|{}_{x=l}=0‚\\ \left(\frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}+{\rm P}\frac{\text{d}y}{\text{d}x}\right)|{}_{x=0}=0‚\\ \frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}|{}_{x=l}=0。(2.2)\end{array}$

显然, 如果某一复数λ和y (x) 是上式的解, 则y (x) ·exp (λt) 必是前节得到的运动方程组 (1.11) , (1.13) , (1.14) 的一个特解。

下面应用希尔伯特空间的算子理论来讨论本征值问题。记 (2.1) 中的微分算符为:

$A=\frac{1}{\rho (x)}\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^2}{\text{d}x{}^2}+\frac{1}{\rho (x)}\frac{\text{d}}{\text{d}x}{\rm N}(x)\frac{\text{d}}{\text{d}x}。(2.3)$

当然不是随便什么函数都能为算符A所作用。从物理概念出发, 需要规定算符A的作用范围。为此引用下列定义:一切以ρ (x) 为权的定义于 (0, l) 上的平方可积的复值函数全体构成完备的内积空间称为希尔伯特空间, 并记为L²_ρ (0, l) 。在这个空间中两个函数y₁ (x) 和y₂ (x) 的内积定义为:

$\begin{array}{l}\langle y{}_1(x),y{}_2(x)\rangle ={\displaystyle {\int }_0^{l}y}{}_1(x)\overline{y{}_2(x)}\rho (x)\text{d}x‚(2.4)\\ \text{■}y(x)\text{■}=\langle y(x)‚y(x)\rangle {}^{\frac{1}{2}}=\\ \left({\displaystyle {\int }_0^l|}y(x)|{}^2\rho (x)\text{d}x\right){}^{\frac{1}{2}}。(2.5)\end{array}$

(2.5) 式称为函数y (x) 的范数。函数y (x) 属于L${}_{\rho }^2$ (0, l) 这一事实记为y (x) ∈L${}_{\rho }^2$ (0, l) 。

空间L${}_{\rho }^2$ (0, l) 中一切能为微分算符作用并满足边界条件的函数全体构成微分算子A的定义域, 记为:

$\begin{array}{l}D(A)=\{y(x)|y\in \text{L}{}_{\rho }^2(0‚l)；y^{\prime }、EJ(x)y^{″}、\\ (EJ(x)y^{″}{)}^{\prime }\in AC；\\ EJ(x)y^{″}|{}_{x=0,l}=0‚(EJ(x)y^{″}{)}^{\prime }|{}_{x=l}=0‚\\ [(EJ(x)y^{″}{)}^{\prime }+{\rm P}{y}^{\prime }]|{}_{x=0}=0\}‚(2.6)\end{array}$

式内y∈AC表示y (x) 是绝对连续的函数, 因而是一次可微的, 它的微商必属于L${}_{\rho }^2$ (0, l) 。这样定义的算子A和D (A) 是L²_ρ (0, l) 中的稠定线性算子。

首先我们证明, 按上述定义的线性微分算子A是自伴的, 这是它最重要的特征。设φ (x) ∈D (A) , ψ (x) ∈D (A^*) , D (A^*) 是A的伴随算子A^*的定义域。依伴随算子的定义 ^[3], 有恒等式成立:

$\langle A\phi ‚\psi \rangle =\langle \phi ‚A{}^{*}\psi \rangle ‚\phi \in D(A)‚\psi \in D(A{}^{*})。(2.7)$

为了证明A是自伴算子, 我们计算上式左端。用分部积分法可写出下列等式:

$\begin{array}{l}\langle A\phi ‚\psi \rangle ={\displaystyle {\int }_0^l\frac{1}{\rho (x)}}(\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^2\phi }{\text{d}x{}^2}+\\ \frac{\text{d}}{\text{d}x}{\rm N}(x)\frac{\text{d}\phi }{\text{d}x})\overline{\psi (x)}\rho (x)\text{d}x=\\ \left[(\frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}{}^2\phi }{\text{d}x{}^2})\overline{\psi (x)}\right]|{}_0^l-\end{array}$

$\begin{array}{l}\left(EJ(x)\frac{\text{d}{}^2\phi }{\text{d}x{}^2}\cdot \frac{\text{d}\overline{\psi (x)}}{\text{d}x}\right)|{}_0^l+\\ \left(\frac{\text{d}\phi }{\text{d}x}\cdot EJ(x)\frac{\text{d}\overline{\psi }}{\text{d}x}\right)|{}_0^l-\end{array}$

$\begin{array}{l}\left(\phi (x)\cdot \frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}{}^2\overline{\psi (x)}}{\text{d}x{}^2}\right)|{}_0^l+\\ \left({\rm N}(x)\frac{\text{d}\phi }{\text{d}x}\cdot \overline{\psi (x)}\right)\end{array}$

${}_0^l-\left({\rm N}(x)\cdot \phi (x)\frac{\text{d}\overline{\psi }}{\text{d}x}\right)$

$\begin{array}{l}{}_0^l+\\ {\displaystyle {\int }_0^l\phi }(x)\left(\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^2\overline{\psi }}{\text{d}x{}^2}+\frac{\text{d}}{\text{d}x}{\rm N}(x)\frac{\text{d}\overline{\psi }}{\text{d}x}\right)\text{d}x。(2.8)\end{array}$

注意到上式右端最后一项正是〈φ, A^*ψ〉。用F表示上式右端前六项的和, 注意函数N (x) 是绝对连续和可微的, 满足条件N (0) =P, N (l) =0。容易检查, 为了对任何φ (x) ∈D (A) 都有F=0, 必须且只须ψ (x) ∈D (A) 。另方面, 从 (2.8) 右端最后一项的表达式中又可看出伴随算子A^*的微分表达式也是由 (2.3) 定义的A, 即在D (A) 上A的表达式和A^*的表达式相同。于是D (A) =D (A^*) , A=A^*, A是自伴算子。

根据希尔伯特空间中的无界线性算子理论 ^[3,4]可作出下列结论:

1) 算子A的所有本征值λ${}_j^2$都是实数, 因而本征值问题 (2.1) 、 (2.2) 的所有本征值都位于实轴上。如果λ${}_j^2$<0, 那末λ_j=iω_j是纯虚数, φ_j (x) ·exp (iω_jt) 是系统可能产生的一种横向振动。如果在某种轴向力N (x) 作用下出现一个或几个本征值λ²>0, 那末至少出现一个λ>0的特解, 它对应于一种振型运动φ (x) ·exp (λt) , 此时, 当t→∞时exp (λt) →∞, 这个振型就是不稳定的。下面将看到, 只要发动机推力P足够大, 这种情况必然会发生。

2) 不同的本征值λ${}_i^2$和λ${}_j^2$, i≠j所对应的本征函数φ_i (x) 和φ_j (x) 按L${}_{\rho }^2$ (0, l) 中的内积是相互直交的, 即

$\begin{array}{l}\langle \phi {}_i(x)‚\phi {}_j(x)\rangle =\\ {\displaystyle {\int }_0^l\phi }{}_i(x)\overline{\phi {}_j(x)}\rho (x)\text{d}x=0‚i\ne j‚(2.9)\end{array}$

这就是说, 不同的振型在 (2.9) 的意下互相直交, 这一点与两端自由的梁的情况是相同的。

3) 不难证明, 算子A一定有无穷多个本征值, 如果按它的绝对值自小至大排列起来, 那末当n→∞时, λ${}_n^2$→∞, 这从物理概念上是极易理解的。对这一事实的严格证明只需应用算子的比较定理即可 ^[5], 这里不再赘述。

为了说明在推力 (阻力) 存在时飞行器横向运动的特点, 再回头讨论 (2.1) 式。设φ (x) 是A的定义域中的任一函数, 取φ (x) =1, 对 (2.1) 式两端分别取内积, 用分部积分法作简单运算后有:

$\begin{array}{l}\lambda {}^2=-\langle A\phi ‚\phi \rangle =-{\displaystyle {\int }_0^l\frac{1}{\rho (x)}}(\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^2\phi }{\text{d}x{}^2}+\\ \frac{\text{d}}{\text{d}x}{\rm N}(x)\frac{\text{d}\phi }{\text{d}x})\overline{\phi (x)}\rho (x)\text{d}x=\\ -{\displaystyle {\int }_0^{l}E}J(x)\left|\frac{\text{d}{}^2\phi }{\text{d}x{}^2}\right|{}^2\text{d}x+{\displaystyle {\int }_0^l{\rm N}}(x)\left|\frac{\text{d}\phi }{\text{d}x}\right|{}^2\text{d}x(2.10)\end{array}$

与没有推力或阻力作用的两端自由的梁比较, 上式右端多了一项与N (x) 有关的项, 它的存在使本征值向正的方向移动。只要函数N (x) 足够大 (只要增大推力就够了) , 就有可能出现λ²>0的情况, 这将对应不稳定的横向运动。这一事实早就有人指出过 ^[6,7]。这种情况当然是要努力避免的, 否则飞行器的横向振动就要失稳发散。这一点在下一节还要进一步讨论。

为了讨论方便, 有时把本征值问题写成一阶线性方程组的形式。参照第一节内的表达式 (1.16) , 置q=0有:

$\lambda y=z,\lambda z=-Ay；$

或者写成向量形式

$\lambda \left(\begin{array}{l}y\\ z\end{array}\right)=\left(\begin{array}{cc}0& 1\\ -A& 0\end{array}\right)\left(\begin{array}{l}y\\ z\end{array}\right)(2.11)$

于是二次本征值问题 (2.1) 就变为矩阵算子

$\mathit{A}=\left(\begin{array}{cc}0& 1\\ -A& 0\end{array}\right)(2.12)$

的一次本征值问题。如果再采用向量写法, 令φ={y, z}, 则 (2.11) 可以表达得更简单

$\lambda \phi =A\phi 。(2.13)$

一个简支梁两端受压时, 当压力大到一定程度就会发生失稳现象。对细长体飞行器来说, 当推力或阻力增大到一定程度也可能发生失稳现象。以后将看到, 对同一个细长体, 失稳推力 (阻力) 要比失稳压力 (两端受压) 大好几倍。下面将证明, 对均匀飞行器, 前者比后者约大2.6倍。

从 (2.10) 式可以看到, 当N (x) 很小时, 例如接近于零, 此时本征值将主要由弹性项决定。但是, 当N (x) 增大时, 某些λ²将由负值逐渐增大, 经过零点移向正值。因此, 使某些本征值λ²=0是飞行器开始失稳的标志。我们的目的是找到这种临界推力 (阻力) , 一旦推力或阻力大于这个临界值后, 就会出现λ²>0的本征值, 使飞行器的横向运动变得不稳定。所以, 应该找到算子A的所有的零本征值对应的非零本征元 (函数) , 也就是求出算子A的全部零子空间。然而, A的零子空间中包括刚体运动对应的本征函数, 它是φ (x) =c, c为常数, 与没有推力或阻力作用的自由飞行器不同, 这里φ (x) =a+bx这种线性函数不再是零本征值对应的本征元了。因为它不满足边界条件。

现在研究在什么条件下会出现除刚体运动以外的零本征值对应的非零本征元。根据上述理由, 在 (2.1) 式中置λ²=0, 还假定轴向力函数N (x) 中的阻力部分

$\frac{D(x)}{D{}_0}={\displaystyle {\int }_0^x\frac{C(\xi )}{D{}_0}}\text{d}\xi$

不依赖于t, 即相对阻力分布规律不随时间变化, 在函数N (x) 的表达式 (1.7) 中设各参数都不依赖于t, 再用η代替P, 用ξ代替总阻D₀, 就得到一个新的方程组:

$\begin{array}{l}\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}+\eta \frac{\text{d}}{\text{d}x}\left(1-\frac{G(x)}{G{}_0}\right)\frac{\text{d}y}{\text{d}x}+\\ \xi \frac{\text{d}}{\text{d}x}\left(\frac{G(x)}{G{}_0}-\frac{D(x)}{D{}_0}\right)\frac{\text{d}y}{\text{d}x}=0‚(3.1)\\ \left(EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}\right)|{}_{x=0}=\left(EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}\right)|{}_{x=l}=0‚\\ \left(\frac{\text{d}}{\text{d}x}EJ(x)+\eta \frac{\text{d}y}{\text{d}x}\right)|{}_{x=0}=0,(3.2)\\ (\frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2})|{}_{x=l}=0.\end{array}$

根据明显的物理意义, 任何非零正数η和ξ能使上式有y (x) ≠c的解存在都是一组失稳推力和阻力。可以推想, 这种组合是很多的, 甚至有无穷多组, 其中最小的一组是我们最感兴趣的, 它们将是保证飞行器横向运动稳定的基本限制。

如果阻力D₀和C (x) 是已知的 (外形和速度为已知) , 使 (3.1) 出现非常数解的最小推力η叫做固定阻力时的失稳推力。反之, 如果推力η=P是固定的, 使 (3.1) 式有非常数解的总阻力ξ就叫做失稳阻力。

推力远大于阻力时, 即高加速飞行时, 推力对横向振动频率的影响常常是最重要的问题, 下面将主要讨论这个情况。当阻力为零时, (3.1) 式变为

$\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}+\eta \frac{\text{d}}{\text{d}x}\left(1-\frac{G(x)}{G{}_0}\right)\frac{\text{d}y}{\text{d}x}=0。(3.3)$

现在证明, 在零阻力飞行中, 最小失稳推力由下列二次方程的最小非零正本征值η决定:

$\begin{array}{l}\frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}w}{\text{d}x}+\eta \left(1-\frac{G(x)}{G{}_0}\right)w=0‚\\ \left(EJ(x)\frac{\text{d}w}{\text{d}x}\right)|{}_{x=0}=\left(EJ(x)\frac{\text{d}w}{\text{d}x}\right)|{}_{x=l}=0。(3.4)\end{array}$

为证明这一事实, 在 (3.3) 中令$\frac{\text{d}y}{\text{d}x}=w$, 再注意到$\left(1-\frac{G(x)}{G{}_0}\right){}_{x=l}=0$, 则 (3.2) 和 (3.3) 可写为

$\begin{array}{l}\frac{\text{d}}{\text{d}x}\left[\frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}w}{\text{d}x}+\eta \left(1-\frac{G(x)}{G{}_0}\right)w\right]=0‚\\ \left(EJ(x)\frac{\text{d}w}{\text{d}x}\right)|{}_{\begin{array}{l}x=0\\ x=l\end{array}}=0‚\\ (\frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}w}{\text{d}x}+\eta (1-\frac{G(x)}{G{}_0})w)|{}_{\begin{array}{l}x=0\\ x=l\end{array}}=0。(3.5)\end{array}$

如果η和w (x) 是上式的一组非零解, 那末必有

$\frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}w}{\text{d}x}+\eta \left(1-\frac{G(x)}{G{}_0}\right)w=c。$

但是, (3.5) 中的第二组边界条件要求c=0。所以, 只要w (x) 满足 (3.4) , 则 (3.5) 就必被满足, 而不管η是什么值。最后, 由于w (x) 是可微的, 而函数G (x) 是绝对连续的 (可分段光滑) , 所以对 (3.4) 的第一式可以再微分一次, 从而得到 (3.5) 。这样即可看到 (3.4) 的本征值和本征元就是 (3.3) 的本征值和本征元。

最小失稳推力, 即 (3.4) 的最小正本征值η₁, 可以用优化方法用计算机求出。

令$A{}_1=\frac{\text{d}}{\text{d}x}$EJ (x) $\frac{\text{d}}{\text{d}x}‚B{}_1=1-\frac{G(x)}{G{}_0}$, 那末有

$\begin{array}{l}\eta {}_1=\underset{\begin{array}{l}w\in D(A{}_1)\\ w\perp {\rm Z}(A{}_1)\end{array}}{{\displaystyle \min }}\frac{-\langle Aw,w\rangle }{\langle Bw,w\rangle }=\\ \underset{\begin{array}{l}w\in D(A{}_1)\\ w\perp {\rm Z}(A{}_1)\end{array}}{{\displaystyle \min }}\frac{{\displaystyle {\int }_0^{l}E}J(x)\left|\frac{\text{d}w}{\text{d}x}\right|{}^2\text{d}x}{{\displaystyle {\int }_0^l\left(1-\frac{G(x)}{G{}_0}\right)}|w|{}^2\text{d}x}‚(3.6)\end{array}$

式中D (A₁) =ww、$\frac{\text{d}w}{\text{d}x}$、$EJ\frac{\text{d}w}{\text{d}x}\in AC$;

$EJ(x)\frac{\text{d}w}{\text{d}x}$_{x=0, l}=0, 而Z (A₁) 表示A₁的零子空间, 即Z (A₁) ={w|A_lw=0}。

注意本征值问题 (3.4) 所确定的算子A₁也是自伴的, 而且是非负的。B₁是一个正的普通函数, 因而也可看成为一个有界自伴算子。根据一个熟知的关于λB₁φ=-A₁φ这类本征值问题的特点, (3.6) 中第一个等式必然成立 (参看文献[11]第五章第34节中的定理) , 唯一需要检查的是下列等式

$\begin{array}{l}-\langle A{}_{1}w,w\rangle =-{\displaystyle {\int }_0^l\frac{\text{d}}{\text{d}x}}EJ(x)\frac{\text{d}w}{\text{d}x}\cdot \overline{w(x)}\text{d}x=\\ -EJ(x)\frac{\text{d}w}{\text{d}x}\overline{w}(x)|{}_0^l+{\displaystyle {\int }_0^{l}E}J(x)\frac{\text{d}w}{\text{d}x}\cdot \frac{\text{d}\overline{w}}{\text{d}x}\text{d}x=\\ {\displaystyle {\int }_0^{l}E}J(x)\left|\frac{\text{d}w}{\text{d}x}\right|{}^2\text{d}x‚\end{array}$

上式内在分部积分时用到了原方程 (3.4) 的边界条件。

用同样的方法可以讨论问题 (3.1) 中推力为零而仅有阻力的情况。这时只要在 (3.4) 式中用ξ代替η, 用${\rm N}(x)=\frac{G(x)}{D{}_0}-\frac{D(x)}{D{}_0}$代替$\left(1-\frac{G(x)}{D{}_0}\right)$即得到求最小失稳阻力公式。由于阻力和推力在 (3.4) 式中处于对称地位, 故不再重复对失稳阻力的讨论。

一般说来, 任何飞行器的截面不可能是均匀的, 即ρ (x) , EJ (x) 不可能是常数, 而是很不规则的函数, 用解析方法求解前节内讨论的边值本征值问题是毫无希望的。对均匀梁, 情况则完全不同, 可以用解析方法进行详细的讨论。通过对均匀梁的讨论可以从理论上得到启发, 这对了解问题的实质可能是很有帮助的。另外, 有了对均匀梁的精密数字结果, 可以利用自伴算子的比较定理去估计非均匀飞行器的失稳推力或失稳阻力的上下界。从这个意义来说, 弄清均匀梁的情况对具体的非均匀的情况可能提供有用的参考, 这就是本节的目的。

回头讨论本征值问题 (2.1) 、 (2.2) 。设质量密度ρ (x) 和弯曲刚度EJ (x) 都是常数, 仍记为ρ和EJ, 依定义 (1.6) ,

$G(x)=g\rho x,G{}_0=\rho gl‚\frac{G(x)}{G{}_0}=\frac{x}{l}。(4.1)$

引进新变量ξ=x/l, 将这些量代入 (2.1) 、 (2.2) 后得到均匀飞行器横向运动本征值问题:

$\lambda {}^{2}y+\frac{EJ}{\rho l{}^4}\frac{\text{d}{}^{4}y}{\text{d}\xi {}^4}+\frac{{\rm P}}{\rho l{}^2}\frac{\text{d}}{\text{d}\xi }(1-\xi )\frac{\text{d}y}{\text{d}\xi }=0‚(4.2)$

相应的边界条件是

$\begin{array}{l}\frac{\text{d}{}^{2}y}{\text{d}\xi {}^2}|{}_{\xi =0}=\frac{\text{d}{}^{2}y}{\text{d}\xi {}^2}|{}_{\xi =1}=0‚\\ \left(\frac{\text{d}{}^{3}y}{\text{d}\xi {}^3}+\frac{{\rm P}l{}^2}{EJ}\frac{\text{d}y}{\text{d}\xi }\right)|{}_{\xi =0}=\frac{\text{d}{}^{3}y}{\text{d}\xi {}^3}|{}_{\xi =1}=0。\end{array}(4.3)$

上式内依然假定阻力为零, 即仅考虑推力。依前节的证明, 失稳推力由下列边值本征值问题的最小本征值决定 (参看 (3.4) 式) :

$\begin{array}{l}\frac{\text{d}{}^{2}w}{\text{d}\xi {}^2}+\eta \left(\frac{l{}^2}{EJ}\right)(1-\xi )w=0‚\\ \frac{\text{d}w}{\text{d}\xi }|{}_{\xi =0}=\frac{\text{d}w}{\text{d}\xi }|{}_{\xi =1}=0。\end{array}(4.4)$

我们知道上式中第一式的通解是一种贝塞尔函数 (参看文献[8], 问题2.162) 。如果记$\mu =\frac{l{}^2}{EJ}\eta$, 则它的通解是

$\begin{array}{l}w(\xi )=c{}_1\sqrt{1-\xi }J{}_{\frac{1}{3}}\left(\frac{2}{3}\sqrt{\mu }(1-\xi ){}^{3/2}\right)+\\ c{}_2\sqrt{1-\xi }J{}_{-\frac{1}{3}}\left(\frac{2}{3}\sqrt{\mu }(1-\xi ){}^{3/2}\right)‚(4.5)\end{array}$

式中c₁和c₂为任意常数, 为了使它满足边界条件 (4.4) , 就要适当选择这些常数。依贝塞尔函数的定义, J_y (x) 有下列级数表达式:

$J{}_y(x)={\displaystyle \sum _{k=0}^{\infty }\frac{(-1){}^k}{k!}}\frac{(x/2){}^{v+2k}}{\text{Γ}(y+k+1)}‚(4.6)$

式内Γ (t) 是伽玛函数。引进新变量$s=\frac{2}{3}\sqrt{\mu }\cdot (1-\xi ){}^{3/2}‚(4.5)$式又可改写为

$\begin{array}{l}w(s)=c{}_{1}s{}^{\frac{1}{3}}J{}_{\frac{1}{3}}(s)+c{}_{2}s{}^{\frac{1}{3}}J{}_{-\frac{1}{3}}(s)=\\ c{}_{1}w{}_1(s)+c{}_{2}w{}_2(s)。(4.7)\end{array}$

将右端两个函数按 (4.6) 展成级数, 有

$\begin{array}{l}w{}_1(s)=s{}^{1/3}J{}_{\frac{1}{3}}(s)=\frac{s{}^{2/3}}{2{}^{1/3}Г\left(\frac{1}{3}+1\right)}-\\ \frac{s{}^{8/3}}{2{}^{7/3}\cdot Г\left(\frac{1}{3}+2\right)}+\frac{S{}^{14/3}}{2{}^{13/3}Г\left(\frac{1}{3}+1\right)}+\cdots ‚\\ w{}_2(s)=s{}^{1/3}J{}_{-\frac{1}{3}}(s)=\frac{2{}^{1/3}}{Г\left(1-\frac{1}{3}\right)}-\\ \frac{s{}^2}{2{}^{5/2}Г\left(2-\frac{1}{3}\right)}+\frac{s{}^4}{2{}^{11/3}Г\left(3-\frac{1}{3}\right)}-\cdots 。\end{array}$

注意到

$\begin{array}{l}\frac{\text{d}w{}_1}{\text{d}\xi }=\frac{\text{d}w{}_1}{\text{d}s}\cdot \frac{\text{d}s}{\text{d}\xi }=-\frac{2}{3}\frac{\sqrt{\mu }\left(\frac{3}{2}\frac{1}{\sqrt{\mu }}\right){}^{1/3}}{2{}^{1/3}\text{Γ}\left(\frac{1}{3}+1\right)}+{\rm O}(s{}^2)‚\\ \frac{\text{d}w{}_2}{\text{d}\xi }=\frac{\text{d}w{}_2}{\text{d}s}\cdot \frac{\text{d}s}{\text{d}\xi }=\frac{\sqrt{\mu }\left(\frac{3}{2}\frac{1}{\sqrt{\mu }}\right){}^{1/3}}{2{}^{2/3}\text{Γ}\left(2-\frac{1}{3}\right)}s{}^{4/3}+{\rm O}(s{}^{10/3})‚\end{array}$

可以推知w₁ (s) 不满足边界条件 (4.4) 中的任何一个, 所以在 (4.5) 中应取c₁=0。另一个函数w₂ (s) 在s=0处即ξ=1处满足边界条件。为了使它又满足另一个边界条件, 必须适当选择s的值, 使函数$\frac{\text{d}w(s)}{\text{d}\xi }$在ξ=0处等于零。记住s是μ的函数, 而μ又是η的函数, 选择s的值实际上就是选择η的值, 使w₂ (x, η) 对x的一次微商在x=0处等于零。这样我们既找到了 (4.4) 的本征函数w₂ (x) , 又找到了本征值η。

应用贝塞尔函数的微分恒等式 ^[9]

$\frac{\text{d}J{}_v(x)}{\text{d}x}=\frac{v}{x}J{}_v-J{}_{v+1}‚$

可算出

$\begin{array}{l}\frac{\text{d}}{\text{d}s}\left(s{}^{1/3}J{}_{-\frac{1}{3}}(s)\right)=\frac{1}{3}s{}^{-2/3}J{}_{-\frac{1}{3}}(s)+\\ s{}^{1/3}\left(\frac{-1}{3s}J{}_{-\frac{1}{3}}(s)-J{}_{\frac{2}{3}}(s)\right)=-s{}^{1/3}J{}_{\frac{2}{3}}(s)。\end{array}$

因此,

$\frac{\text{d}w{}_2(s)}{\text{d}\xi }=\frac{\text{d}w{}_2(s)}{\text{d}s}\cdot \frac{\text{d}s}{\text{d}\xi }=cs{}^{2/3}J{}_{\frac{2}{3}}(s)‚(4.8)$

式内c为常数, 并利用前面定义了的关系式$s=\frac{2}{3}\sqrt{\mu }(1-\xi ){}^{3/2}$, 故上式内$c=\sqrt{\mu }\left(\frac{3}{2}\frac{1}{\sqrt{\mu }}\right){}^{1/3}$。为了使 (4.8) 式满足ξ=0处的边界条件 (4.4) , 必须选择自变量

$s=\frac{2}{3}\sqrt{\mu }=\frac{2l}{3}\frac{1}{\sqrt{EJ}}\sqrt{\eta }‚(4.9)$

即选择η的值, 使它成为贝塞尔函数$J{}_{\frac{2}{3}}(s)$的零点。由于$J{}_{\frac{2}{3}}(s)$有无穷多个零点, 故本征值问题 (4.4) 有无穷多个不同的本征值和本征函数, 与$J{}_{\frac{2}{3}}(s)$的零点一一对应。

我们最感兴趣的是最小的一个本征值η₁, 它对应$J{}_{\frac{2}{3}}(s)$的最小的那个零点。从贝塞尔函数表中可查得${}^{[10]}‚\frac{2}{3}$阶贝塞尔函数的第一个零点s₁=3.377, 故$\mu {}_1=\left(\frac{3}{2}s{}_1\right){}^2=25.66$, 从而本征值问题 (4.4) 的第一个最小本征值为 1

$\eta {}_1=\mu {}_1\frac{EJ}{l{}^2}=25.66\frac{EJ}{l{}^2}。(4.10)$

这就是均匀飞行器的横向弯曲失稳推力。

为了与上述结果比较, 下面写出两端受压的自由梁在没有轴向过载时的失稳压力。容易证明, 此时失稳压力是下列本征值问题的最小值

$EJ\frac{\text{d}{}^{2}w}{\text{d}x{}^2}+{\rm P}w=0‚\frac{\text{d}w}{\text{d}x}|{}_{x=0}=\frac{\text{d}w}{\text{d}x}|{}_{x=l}=0。(4.11)$

显然, 它的通解是w (x) $=c\cdot \cos \mu x‚\frac{\text{d}w}{\text{d}x}=c{}_1\sin \mu x$。为使w (x) 满足两个边界条件, 必须选取μ₁l=π·n, 我们应该选最小的一个n=1。注意到$\mu {}_1=\sqrt{\frac{{\rm P}{}_1}{EJ}}$, 显然有

${\rm P}{}_1=\frac{\text{π}{}^{2}EJ}{l{}^2}‚(4.12)$

这就是两端受压自由梁的最小失稳压力。比较最小失稳推力和最小失稳压力

$\frac{\eta {}_1}{{\rm P}{}_1}=\frac{25.66}{\text{π}{}^2}=2.599。$

在本节开始时我们已经提到了这一比值, 即失稳推力是失稳压力的2.6倍。

大的轴向过载, 无论是由于大推力或者是由于大阻力引起的, 对细长体飞行器横向运动的影响, 用地面模拟试验的办法是很难测定的, 因为要造成相应的环境条件很不容易。因此, 在预先设计阶段能够给出所讨论的影响数量的估算方法是很有用处的。这一节内我们试给出一个估计式, 在小过载和小阻力条件下可望是比较准确的估计。而大过载条件下用一次逼近公式去估算, 误差可能较大。

推力和阻力对横向振动频率影响趋势从 (2.10) 式可以看得很清楚:只要有轴向力N (x) 存在, 各阶振型的频率都要降低, 轴向力越大, 降低越多。只要推力达不到失稳推力, 自由弹性体的各阶振型都可能发生, 但频率要向低频方向移动。一旦推力达到失稳推力的程度, 第一阶振型的振动或者消失, 或者变为不稳定振动。下面我们将假定推力小于失稳推力, 从而估算第一阶振型频率降低的数量。

先设推力和阻力都为零, 这时本征值问题 (2.1) 、 (2.2) 将变为通常自由梁的本征值问题,

$\begin{array}{l}\lambda {}^{2}y(x)+\frac{1}{\rho (x)}\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}=0‚\\ \left(EJ(x)\frac{\text{d}x{}^2}{\text{d}{}^{2}y}\right)|{}_{x=0,l}=0；(5.1)\\ \left(\frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}\right)|{}_{x=0,l}=0。\end{array}$

由上式决定的微分算子记为A₀,

$\begin{array}{l}A{}_0=\frac{1}{\rho (x)}\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^2}{\text{d}x{}^2}‚\\ D(A{}_0)=\{y(x)|y(x)、\frac{\text{d}y}{\text{d}x}、EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}、\\ \frac{\text{d}}{\text{d}x}EJ(x)\frac{\text{d}{}^{2}y}{\text{d}x{}^2}\in AC‚\\ A{}_{0}y(x)\in \text{L}{}_{\rho }^2(0‚l)\}。(5.2)\end{array}$

再假定 (5.1) 的第一个 (最小的) 非零本征值为λ₁₀, 相应的本征函数 (振型) 为φ₁₀ (x) , 那末依 (5.1) 式有

$\begin{array}{l}\langle \phi {}_{10}‚\phi {}_{10}\rangle \lambda {}_{10}^2=-\langle A{}_0\phi {}_{10}‚\phi {}_{10}\rangle =\\ -{\displaystyle {\int }_0^l\frac{1}{\rho (x)}}\frac{\text{d}{}^2}{\text{d}x{}^2}EJ(x)\frac{\text{d}{}^2\phi {}_{10}}{\text{d}x{}^2}\cdot \overline{\phi {}_{10}(x)}\rho (x)\text{d}x=\\ -{\displaystyle {\int }_0^{l}E}J(x)\left|\frac{\text{d}{}^2\phi {}_{10}}{\text{d}x{}^2}\right|{}^2\text{d}x。\end{array}$

所以,

$\begin{array}{l}\lambda {}_{10}^2=\frac{-{\displaystyle {\int }_0^{l}E}J(x)\left|\frac{\text{d}{}^2\phi {}_{10}}{\text{d}x{}^2}\right|{}^2\text{d}x}{{\displaystyle {\int }_0^l|}\phi {}_{10}(x)|{}^2\rho (x)\text{d}x}=\\ \begin{array}{l}\max \\ \begin{array}{l}\phi \in D(A{}_0)\\ \phi \perp {\rm Z}(A{}_0)\end{array}\end{array}[JX-*2]\frac{-{\displaystyle {\int }_0^{l}E}J(x)\left|\frac{\text{d}{}^2\phi }{\text{d}x{}^2}\right|{}^2\text{d}x}{{\displaystyle {\int }_0^l|}\phi (x)|{}^2\rho (x)\text{d}x}。(5.3)\end{array}$

上式右端是关于自伴算子本征值的极值特性的表达式。

用λ${}_1^2$和φ₁ (x) 分别表示在推力P的作用下飞行器第一个振型对应的 (2.1) 、 (2.2) 的本征值和本征函数, 用f₁和f₁₀分别表示有推力和无推力飞行器的一阶振型频率。下面证明一个主要不等式:

$\text{Δ}f=f{}_{10}-f{}_1\ge \frac{{\rm P}{\displaystyle {\int }_0^l\left(1-\frac{G(x)}{G{}_0}\right)}\left|\frac{\text{d}\phi {}_{10}}{\text{d}x}\right|{}^2\text{d}x}{8\pi {}^{2}f{}_{10}{\displaystyle {\int }_0^l|}\phi {}_{10}(x)|{}^2\rho (x)\text{d}x}。(5.4)$

我们看到, 上式右端要求事先知道φ₁₀ (x) 和f₁₀及推力P, 根据这些数据即可估算出无推力时的一阶振型频率f₁₀和有推力时的一阶振型频率f₁之差。右端所要求的这些数据都可在地面试验中准确地测定出来。

我们用算子比较法来证明 (6.4) 式确实成立 ^[11]。为使两种不同的本征值问题 (2.1) 和 (5.1) 进行比较, 必须先把 (2.1) 中的算子A和 (5.1) 中的算子A₀的定义域联系起来, 因为按定义D (A) 和D (A₀) 是不相同的, 差别在于属于D (A) 的函数在x=0附近的边界条件与D (A₀) 的函数不同。为了能够直接比较两类函数, 今对D (A₀) 中的函数加以改造。设ε>0是足够小的正数, 构造函数h (x)

$h(x)=\{\begin{array}{l}\exp \frac{(x-\frac{\epsilon }{2}){}^2}{(x-\frac{\epsilon }{2}){}^2-\frac{\epsilon {}^2}{4}},0\le x\le \epsilon ；\\ 0‚\epsilon \le x\le l。\end{array}(5.5)$

显然, h (x) 和它的所有微商都是连续的, 而且在x=0和x=ε处h (x) 和它的各阶导数均等于零。

假定EJ (x) 做为x的函数, 在区间[0, ε]上是足够光滑的。今定义

$\tilde{\phi }{}_{10}(x)=\phi {}_{10}(x)-\eta (x)‚(5.6)$

式中

$\eta (x)=\{\begin{array}{l}\phi {}_{10}(x)-\phi {}_{10}(\epsilon )(1-h(x))‚0\le x\le \epsilon ；\\ 0‚\epsilon \le x\le l.\end{array}(5.7)$

不难检查, $\tilde{\phi }$₁₀ (x) 满足边界条件 (2.2) , 因而$\tilde{\phi }{}_{10}(x)\in D(A)$。按照自伴算子本征值的极值特性, 有

$\begin{array}{l}-\lambda {}_1^2=\underset{\begin{array}{l}\phi \in D(A)\\ \phi \perp {\rm Z}(A)\end{array}}{{\displaystyle \min }}\frac{\langle A\phi ‚\phi \rangle }{\langle \phi ‚\phi \rangle }\le \frac{\langle A\tilde{\phi }{}_{10}‚\tilde{\phi }{}_{10}\rangle }{\langle \tilde{\phi }{}_{10}‚\tilde{\phi }{}_{10}\rangle }=\\ \frac{{\displaystyle {\int }_0^{l}E}\eta (x)\left|\frac{\text{d}{}^2\tilde{\phi }{}_{10}}{\text{d}x{}^2}\right|{}^2\text{d}x-{\rm P}{\displaystyle {\int }_0^l\left(1-\frac{G(x)}{G{}_0}\right)}\left|\frac{\text{d}\tilde{\phi }{}_{10}}{\text{d}x}\right|{}^2\text{d}x}{{\displaystyle {\int }_0^l|}\tilde{\phi }{}_{10}(x)|{}^2\rho (x)\text{d}x}\end{array}$