2 升力线理论

2 升力线理论

2.1 减阻

阻力

什么是阻力？阻力是阻止主要运动（位移向量）的力。

可以用一个简单的公式描述阻力：

\[\begin{equation} \overrightarrow{R_2}-\overrightarrow{R_1}\propto\vec{T}-\vec{D} \end{equation} \]

这里的R是反作用力（reactive force），T是推力（thrust），D是阻力（drag）。

（猜测：这个公式表明，飞行器在某种状态变化时（如加速、改变姿态），反作用力的变化量，与推力和阻力的差值成正比）

为简化考虑，我们基于黏性讨论阻力。

如下面的流程图所示，我们把阻力按粘性和非粘性分类。

​​

概念

流线和切线条件

流线（streamline）是一个跟踪流动的曲线，速度向量始终是该流线的切线。

​​

关于流线，考虑切线条件（tangency condition），可用如下公式描述：

\[\begin{equation} \frac{v}{u}=(\frac{dy}{dx})_{streamline} \end{equation} \]

迹线（pathline）总是跟踪单个粒子的流动。

角速度（旋转）与应变

注意：这一节是我猜的，如果有哪位朋侑知道Richard在课堂上是怎么讲这一部分的，欢迎评论或联系我进行修改。

我们首先引入一些符号，考虑一个三维流动。

定义：\(\lambda\)表示流体微元绕z轴的总旋转角度；\(\delta\)表示流体微元在xy平面内的总剪切位移，所导致的夹角变化（下面会对该定义进行解释）。

给出如下公式：

\[\begin{equation} \frac{D\lambda}{Dt}=\omega_z \end{equation} \]

\[\begin{equation} \frac{D\delta}{Dt}=-\epsilon_{xy} \end{equation} \]

\[\begin{equation} \omega_z=\frac{1}{2}\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right)\quad\epsilon_{xy}=\frac{1}{2}\left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right) \end{equation} \]

也就是说，总旋转角度\(\lambda\)可以视为角速度\(\omega_z\)的积分，总剪切应变\(\delta\)可以表示为应变率\(\epsilon_{xy}\)的负积分。

关于总剪切应变一阶导中的负号：

• 我们假设有一个正方形的单元体，它两条边之间的夹角是90度；
• 当流体微元开始运动，发生剪切变形时，夹角变小（被拉伸成一个菱形）；
• \(\delta\)是当前时刻下，两条边之间的夹角与初始状态的差值；
• 约定：正应变代表角度减小。

涡量

旋转是一个矢量，不仅有大小，还有旋转轴的方向。

在笛卡尔坐标系下，可将角速度矢量分解为三个分量：

\[\begin{equation} \vec{\omega}=\omega_x\vec{i}+\omega_y\vec{j}+\omega_z\vec{k} \end{equation} \]

我们引入涡量（vorticity）的定义：

\[\begin{equation} \vec{\Omega}=2\vec{\omega}=\nabla\times\vec{V} \end{equation} \]

这实际上就是速度场的旋度，即矢量做旋转运动的方向和强度。

回顾一下旋度的计算公式：

\[\begin{equation} \begin{aligned} \nabla \times \mathbf{v} &= \left( \frac{\partial}{\partial x} \mathbf{i} + \frac{\partial}{\partial y} \mathbf{j} + \frac{\partial}{\partial z} \mathbf{k} \right) \times (u \mathbf{i} + v \mathbf{k} + w \mathbf{k}) \\ &= \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ u & v & w \end{vmatrix} \\ &= \left( \frac{\partial w}{\partial y} - \frac{\partial v}{\partial z} \right) \mathbf{i} + \left( \frac{\partial u}{\partial z} - \frac{\partial w}{\partial x} \right) \mathbf{j} + \left( \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \right) \mathbf{k} \end{aligned} \end{equation} \]

回顾公式5，我们可得：

\[\begin{equation} \nabla\times\vec{V}=\vec{\Omega}=2\omega_{x}\vec{i}+2\omega_{\nu}\vec{j}+2\omega_{z}\vec{k} \end{equation} \]

下图展示了一个点涡周围，散度为0的有旋矢量场。

​​

环量

环量（circulation）的定义是：沿着一个闭合路径C的速度矢量\(\vec{V}\)的切向分量的总和。

另一种理解是：沿着给定闭合路径C的总的旋转强度。

环量的公式为：

\[\begin{equation} \Gamma=-\oint_C\vec{V}\cdot\vec{ds}=-\oint_C \left|\vec{V}\right|\left|\vec{ds}\right|\cos(\theta) \end{equation} \]

这里：\(\Gamma\)念作Gamma，\(\vec{V}=u\vec{i}+v\vec{j}+w\vec{k}\)，\(\vec{ds}=dx\vec{i}+dy\vec{j}+dz\vec{k}\)。

下图展示了一个环量的参考图象。

​​

流函数

流函数（stream function）的目的是绘制稳态流动中的粒子轨迹。

我们希望使用一个标量函数来描述速度场。

定义（考虑二维不可压缩流动）：

\[\begin{equation} \vec{V}=\psi(x,~y) \end{equation} \]

回顾一下CMT中的连续性方程：在不可压缩流动中，速度场的散度为零，流体体积保持不变，流入和流出的速度平衡。

（连续性方程就是速度场的散度）

\[\begin{equation} \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0 \end{equation} \]

定义流函数的偏导数与速度分量的关系：

\[\begin{equation} u=\frac{\partial\psi}{\partial y}\quad v=-\frac{\partial\psi}{\partial x} \end{equation} \]

代入到式12，有：

\[\begin{equation} \frac{\partial^2\psi}{\partial x\partial y}+\frac{\partial^2\psi}{\partial y\partial x}=0 \end{equation} \]

例如：给定流函数\(\psi(x,y)=\frac{1}{2}xy^2-\frac{1}{3}x^3\)，计算得到

\[\begin{equation} u=xy\quad\mathrm{and}\quad v=-\frac{1}{2}y^2+x^2 \end{equation} \]

据此可绘制图象。

​​

速度势

速度势（velocity potential）仅适用于无旋流动，即涡量=0。

回顾一下大三的数学课，任何标量函数的梯度的旋度恒为0：

\[\begin{equation} \nabla\times\nabla\phi=0 \end{equation} \]

定义一个标量函数\(\phi(x,~y)\)，使得速度矢量V等于\(\phi(x,~y)\)的梯度：

\[\begin{equation} \vec{V}=\nabla\{\phi(x,y)\} \end{equation} \]

我们把速度势与速度分量联系起来：

\[\begin{equation} u\vec{i}+v\vec{j}=\frac{\partial\phi}{\partial x}\vec{i}+\frac{\partial\phi}{\partial y}\vec{j} \end{equation} \]

这样也就是让\(u=\frac{\partial \phi}{\partial x}\)，\(v=\frac{\partial \phi}{\partial y}\)。

在无旋流动中，流体微团从一点移动到另一点所做的功，只和这两点的位置有关，而和路径有关，这就是速度势的物理意义。

库塔-茹科夫斯基理论与条件

对于理想不可压缩流体定常流动，在有势力作用下，直均流绕过任意截面形状物体，有环量的流动，将受到垂直于来流方向的升力（侧向力），如下图所示。

​​

其大小为：

\[\begin{equation} L=\rho V_{_\infty}\Gamma \end{equation} \]

其中，L为作用在绕流物体上的升力，\(V_{\infin}\)是来流速度，\(\Gamma\)是绕流物体的速度环量。

注意：该理论考虑的是二维情况，这里给出的是单位展长的升力。

我们知道，对于一个给定的翼型，在一定的迎角下，升力是唯一确定的。也就是说，此时仅存在一个确定的环量值。

确定这个环量的条件叫库塔条件。它的表述为：

对于给定的翼型和迎角，上下翼面的气流应在后缘处平顺离开，且后缘速度值保持有限。也可以说，绕翼型的环量值应使平顺离开后缘。

​​

接下来，我们深入讨论一下该理论下升力的形成机制。

​​

如上图所示。一个有环量的圆柱绕流流场，可以表示为一个没有环量的圆柱绕流的对称流场，与一个由旋转圆柱诱导出的反对称流场叠加而成。

这样，叠加的结果就是圆柱下面速度小、上面速度大的不对称绕流，然后根据伯努利方程，得到圆柱下面压强大、上面压强小，导致向上的升力。

而对于飞机翼型而言，虽然它没有直接旋转，但是通过改变迎角、翼型上下表面的不对称，也相当于给一个对称的绕流场，叠加了一个反对称的流场，导致总体流场下慢上快，进而得到机翼下面压强大、上面压强小，产生向上的升力。

翼尖涡和诱导阻力

​​

上面我们提到，飞机机翼上表面压强较低、下边压强较高。

在机翼翼尖处，高压区的空气会绕流到低压区，形成旋转的涡流，叫做翼尖涡流（wing tip trailing vortices）。

叠加考虑飞机向前运动的速度，翼尖涡流会在机翼后方，诱导（induce）产生一个向下的速度场，称为下洗（downwash）。

下洗的气流会使得机翼局部区域的相对气流方向向下偏，导致升力方向向后倾斜。

升力矢量后倾，会产生一个与飞行方向相反的分量，这就是诱导阻力（induced drag）。

​​

如上图所示，实际上，整个机翼的后缘都有向后的涡量分量，翼尖处的涡强度最大，越往机翼中间走，涡强度越小。它们共同构成了一个称为涡流片（vortex sheet）的结构。

下一步，关于计算流场中由旋涡存在产生的诱导速度，我们使用毕奥-萨瓦尔公式。

​​

如上图所示，在不可压缩流动中，强度为\(\Gamma\)（还记得之前讨论的环量吗？它是沿着给定闭合路径C的总的旋转强度）、微段长度为dL的涡线，对周围流场（这里讨论距离涡的微段距离为r的点M处）所产生的诱导速度dw为：

\[\begin{equation} d\vec{w}=\frac{\Gamma}{4\pi}\frac{d\vec{L}\times\vec{r}}{r^3}\quad dw=\frac{\Gamma}{4\pi}\frac{\sin\alpha}{r^2}dL \end{equation} \]

简化一下，我们当前考虑的机翼是有限长度的平直翼。

​​

那么，在有限长度的直线涡情况下，诱导速度为：

\[\begin{equation} w=\frac{\Gamma}{4\pi h}\int_{\alpha_1}^{\alpha_2}\sin\alpha d\alpha=\frac{\Gamma}{4\pi h}(\cos\alpha_1-\cos\alpha_2) \end{equation} \]

继续简化一下，考虑一端伸向无穷远的直线涡。

​​

此时的速度为：

\[\begin{equation} w=\frac{\Gamma}{4\pi h} \end{equation} \]

我们很快就会用到上述简化定义。

升力线模型

如何对低速机翼的气动特性进行建模呢？

对于机翼绕流，我们可以用一根具有强度\(\Gamma(z)\)的，直的附着涡线，和从附着涡向下游拖出的自由涡系代替。

对于尾流的自由涡系，假设自由涡系在机翼所处平面上，由许多根轴线平行于来流、伸向下游无穷远处的直涡线所组成。

对于大展弦比机翼，自由涡面的卷起和弯曲，主要发生在远离机翼的地方。

我们添加一个假设：自由涡面既不卷起也不耗散。那么，直匀流、绕大展弦比直机翼流动的气动模型可以建立为：

直匀流+附着涡面+自由涡面

​​

机翼上的是附着涡面，尾流是自由涡面，顶上的\(V^-\)是直匀流。

这样的“马蹄涡”的来源是，把一个无限长度的直线涡流，两端向下折90度，长得像字母\(\Pi\)（念Pi）。

这个模型的物理作用是：

1. 马蹄涡中，垂直于来流的部分是附着涡系，可以代替机翼的升力作用；
2. 沿着翼展方向（展向）各剖面上，通过的涡线数量不同：中间剖面通过的涡线最多，环量最大；翼尖/翼端剖面没有涡线通过，环量为零。
   这样就模拟了环量和升力的展向分布；
3. 展向相邻两个剖面之间，拖出的自由涡强度，等于这两个剖面上，附着涡的环量差。

进一步简化：由于大展弦比直机翼的弦长远小于展长，可以近似把机翼上的附着涡系合并，变成一条在展向变强度的附着涡线，各剖面的升力就作用在该线上。

这就是升力线假设。

此时气动模型简化为：

\[\textbf{直匀流+附着涡线+自由涡面} \]

考虑到低速翼型的升力增量约在\(\frac{1}{4}\)弦点处，因此附着涡线放在展向各剖面的\(\frac{1}{4}\)弦点的连线上，这就是升力线。

​​

​​

上述两图形象地描述了升力线模型。

2.2 升力线模型的公式

基本

下洗速度和下洗角

关于升力线公式的推导，我们需要结合下面这两个图来看。其中坐标系与变量的设定以图1为准，图2供某些特殊场景下的说明。

​​

图1

​​

图2

我们设附着涡线在展向位置y处，强度为\(\Gamma({y})\)；在y+dy处，涡强度为\(\Gamma({y})+\frac{d\Gamma}{dy} dy\)。

这里，关于y+dy处的涡强度，实际上是对\(\Gamma (y+\Delta y)\)进行一阶泰勒展开的结果。

根据涡量守恒定理，dy微段拖出的自由涡强度为\(\frac{d\Gamma}{dy} dy\)。

这里的涡量守恒定理（又称为亥姆霍兹定理）是这么说的：涡线不能在流体内部终止，必须形成闭合回路，或延伸到无穷远。

对于微元段 dy，在 y 处有一条强度为\(\Gamma({y})\)的附着涡进入该微元段，在 y + dy 处有一条强度为\(\Gamma({y})+\frac{d\Gamma}{dy} dy\)的附着涡离开该微元段。为了保持涡的闭合性，必须有一条强度为\(\frac{d\Gamma}{dy} dy\)的涡线离开这个微元段。

这里特别注意：图1中的h等于图2中的\(\zeta - z\)，而图2中的\(d\zeta\)等于我们假设的微段dy，dy的初始位置应当是从图1中的y处向翼尖方向（y指向的正方向）向外。

这样，该自由涡在附着涡线上任一点\(y_0\)处的下洗速度为：

\[dw=\frac{d\Gamma}{4\pi h}=\frac{\frac{d\Gamma}{dy}dy}{4\pi(y-y_0)} \]

对上式积分，得到整个涡系在\(y_0\)处产生的下洗速度：

\[\begin{equation} w(z)=\int_{-b/2}^{b/2}\frac{\frac{d\Gamma}{dy}}{4\pi(y-y_0)}dy \end{equation} \]

由于下洗速度的存在，机翼展向每个剖面上的实际有效速度（\(V_e\)）和有效迎角（\(\alpha_e\)）比来流速度和几何迎角小。

​​

如图所示，有效迎角：

\[\begin{equation} \alpha_\mathrm{e}=\alpha-\alpha_\mathrm{i} \end{equation} \]

这里\(\alpha_i\)叫下洗角，其公式为：

\[\begin{equation} \tan\alpha_i=-\frac{w}{V_{\infty}}\approx\alpha_i \end{equation} \]

\[\begin{equation} \alpha(i)=\frac{1}{4\pi V_{\infin}} \int_{-b/2}^{b/2}\frac{\frac{d\Gamma}{dy}}{(y-y_0)}dy \end{equation} \]

接下来，我们讨论作用于该机翼微段的升力。

使用库塔-茹科夫斯基环量定理，求升力和诱导阻力：

\[\begin{equation} \begin{aligned} &L^{\prime}=F^{\prime}\cos\alpha_i\\ &F^{\prime}=\rho V\Gamma\\ &L^{\prime}=\rho V\Gamma\cos\alpha_{i}=\rho V\Gamma\frac{V_{\infty}}{V}\\ &L^{\prime}=\rho V_{\infty}\Gamma\\ &D_i^{\prime}=L^{\prime}\tan\alpha_i\approx L^{\prime}\alpha_i \end{aligned} \end{equation} \]

升力系数与环量

下一步，我们需要求微段当地升力系数。

​​

已知\(C_L\)和\(\alpha\)之间呈线性关系，如下图所示。

​​

\[\begin{equation} \begin{aligned} C_L'(y_0)&=m_0(y_0) \alpha_{true}(y_0)\\ &=m_0(y_0)[\alpha(y_0)-\alpha_i (y_0)-\alpha_0 (y_0)\\ &=m_0(y_0)[\alpha_a(y_0)-\alpha_i (y_0)] \end{aligned} \end{equation} \]

这里，\(\alpha_a(y_0)=\alpha(y_0)-\alpha_0(y_0)\)是绝对迎角。

接下来的推导将不同于Richard上课讲的流程，他是做的数值解，建立环量和攻角之间的关系；我们这里是直接往解析解方向走，会复杂一些。

我们考虑作用于微段dy上的升力dL。回顾一下，升力的公式为\(L=\frac{1}{2}\rho V^2 S C_L\)，那么我们可以把dL表示为：

\[\begin{equation} \begin{aligned} dL&=\frac{1}{2} \rho V_{\infin}(y)^{2} [c(y) dy]C_L'\\ &=\rho V_{\infin} \Gamma(y) dy \end{aligned} \end{equation} \]

这里c(y)表示当地翼弦长度，第二个式子是从环量定理得到的。

两边消元代换一下，可得环量：

\[\begin{equation} \Gamma\left(y\right)=\frac{1}{2}\upsilon_\infty C_{L\infty}^\alpha c\left(y\right) \left[\alpha_a\left(y\right)-\frac{1}{4\pi V_{\infin}} \int_{-b/2}^{b/2}\frac{\frac{d\Gamma}{dy}}{(y-y_0)}dy\right] \end{equation} \]

三角级数展开

变换和展开

若想对本节涉及的傅里叶级数有进一步理解，推荐观看B站视频：

【【官方双语】微分方程概论-第四章：但什么是傅立叶级数呢？-从热流到画圈圈】https://www.bilibili.com/video/BV1vt411N7Ti?vd_source=7e5e34f8276fd18e94d3486b7516c07f

我们已经有了环量公式，但这个公式太复杂了。接下来的任务就是对该公式进行简化。

使用三角级数变量变换：

\[\begin{equation} y_0=-\frac{l}{2}\mathrm{cos}\theta,~y=-\frac{l}{2}\mathrm{cos}\theta_1 \end{equation} \]

回代一下，有：

\[\begin{equation} \Gamma(\theta)=\frac{1}{2}V_\infty C_{L\infty}^\alpha b(\theta)\left[\alpha_a(\theta)+\frac{1}{2\pi V_\infty l}\int_{-l/2}^{l/2}\frac{\frac{\mathrm{d}\Gamma}{\mathrm{d}\theta_1}\mathrm{d}\theta_1}{\cos\theta-\cos\theta_1}\right] \end{equation} \]

使用三角级数展开（暴力地，假设环量\(\Gamma(\theta)\)能用下述三角级数表示）：

\[\begin{equation} \Gamma(\theta)=2lV_\infty\sum_{n=1}^\infty A_n\sin(n\theta) \end{equation} \]

关于这里的常数项\(2lV_{\infin}\)，有几种角度的解释：

• 从量纲分析角度出发，环量的量纲是\(\frac{长度^2}{时间}\)，再乘以无量纲系数2；
• 从环量的公式出发，注意到在迎角项中，分子有\(2V_{\infin}l\)，把常数项设定为这个之后可以约掉。

不过，我们注意到，傅里叶级数理论针对的是定义在有限区间上的周期函数。

在升力线理论中，环量\(\Gamma(\theta)\)是定义在有限的区间\([0,~\pi]\)上的函数，并非周期函数；此外，对于典型机翼，有边界条件：环量在翼尖处趋于0，即\(\Gamma(0)=\Gamma({\pi})=0\)。

奇延拓

那么，我们是如何近似表示环量分布的呢？

我们可以用一些技巧“延拓”非周期函数，使得它在更大的区间上变成周期函数。

具体来说，对于本情境，我们隐含地假设一个在\([-\pi,~\pi]\)区间上定义的函数，它是关于原点奇对称的，把\(\Gamma(\theta)\)看做这个周期为\(2\pi\)的奇函数的半个周期，即\([0,~\pi]\)部分。这种操作叫奇延拓。

从物理上考虑，对于常规的机翼形状和攻角分布，环量分布呈现中间大、两头小的特点，这符合正弦函数一半部分的特性。

但是，为什么这里只出现了正弦函数，而没有余弦函数呢？

这是因为，奇延拓要求，函数在整个延拓区间\([-\pi,~\pi]\)内，满足奇函数性质：\(f(-x)=-f(x)\)。

由于每个\(sin(n\theta)\)都是奇函数，它们的线性组合仍然是奇函数；但是，余弦函数\(cos(n\theta)\)是偶函数，其线性组合也是偶函数，不满足奇延拓的要求。

此外，考虑边界条件，\(sin(n\times 0)=sin(n\pi)=0\)对任意整数n成立，但余弦函数不满足该条件。

求解析解

一般而言，取三角级数的前4项足以近似表示实际的环量分布。

例如，我们在\(\theta = 0\sim\frac{\pi}{2}\)范围内取4个\(\theta\)值，对应右半翼的4个剖面，建立关于系数A的4个代数方程。求解方程即可得四个A值，进而确定环量的分布。

对式34，即环量分布函数积分，乘以密度和速度得升力：

\[\begin{equation} L=\rho V_\infty\int_{-b/2}^{b/2}\Gamma(y)\mathrm{d}y \end{equation} \]

升力系数为：

\[\begin{equation} \begin{aligned} C_L&=\frac{2}{V_\infty S}\int_{-l/2}^{l/2}\Gamma(y)\mathrm{d}y\\ &=\frac{4l}{S}\int_0^\pi\sum_{n=1}^\infty A_n\sin(n\theta)\cdot\frac{l}{2}\sin\theta\mathrm{d}\theta\\ &=\pi\lambda A_1\\ &=\frac{C_{L\infty}^\alpha}{1+\frac{C_{L\infty}^\alpha}{\pi\lambda}(1+\tau)}{\alpha_a} \end{aligned} \end{equation} \]

诱导阻力系数为：

\[\begin{equation} \begin{aligned}{C_{Di}}&=\frac{2}{\upsilon_{\infty}S}\int_{-\frac{b}{2}}^{+\frac{b}{2}}\Gamma(y)\Delta\alpha_{i}(y)dy\\&=\frac{2l^2}{S}\int_0^\pi\left[\sum_{m=1}^\aleph mA_m\sin\left(m\theta\right)\times\sum_{n=1}^\infty A_n\sin\left(n\theta\right)\right]d\theta\\&=\frac{C_y^2}{\pi\lambda}(1+\delta)\end{aligned} \end{equation} \]

这里，\(\tau\)和\(\delta\)都是与机翼平面形状和展弦比有关的修正值，其计算公式为：

\[\begin{equation} \begin{aligned}&\tau=\frac{l}{2S}\int_{0}^{\pi}\frac{\sum_{n=2}^{\infty}nA_{n}\sin\left(n\theta\right)}{A_{1}}b\left(\theta\right)\mathrm{d}\theta\\&\delta=\sum_{n=2}^{\infty}\frac{nA_{n}^{2}}{A_{1}^{2}}=\frac{3A_{3}^{2}}{A_{1}^{2}}+\frac{5A_{5}^{2}}{A_{1}^{2}}+\frac{7A_{7}^{2}}{A_{1}^{2}}+\cdots\end{aligned} \end{equation} \]

升力线理论的数值解

解析解是很难求的，因此我们还是要讨论如何利用计算机求数值解。

由于时间关系，这里先贴出我手写的推导过程，如果有时间会换成markdown。

​​

​​​