流体力学读书笔记

高等流体力学读书笔记

论文题目: 特征线法读书笔记

姓名: 杨志伟

学号: 113108000839

专业:兵器发射理论与技术

指导教师: 周建伟

日期: 2013年12月

1 特征线法

1.1 理论的引出

在考虑了两对面管壁都外折使得两束膨胀波相交，以及膨胀波束在自由边界上反射等问题时，单有绕外钝角流动的公式就不够使用了，需要一种使用于解复杂问题的方法，这就是特征线法。在概况性上说，特征线法确实比绕外钝角的解法进了一步，只要是两个自变数的双曲线型偏微分方程都能用。定常超声速流（包括平面及轴对称的无旋和有旋流）与一维费定常流（不论亚声速还是超声速）的运动方程都是双曲型的。这几种流动能在数学上归在一起，正是反映了在物理上这几种流动都是以波的形式进行变化这样一个事实。

定常亚声速流场上，流动的变化不是以波的形式进行的，任何扰动都没有界线可言，扰动遍及全场，变化都是连续的，任何流动参数（速度、密度和压强等）不仅本身连续，而且它对空间坐标的导数也都连续。与此相反，在定常超声速流场上，扰动都是有界的，像激波在流场中是以突跃面的形式存在的，流动参数本身在突跃面上有突跃的变化，称为强突跃；另一种扰动也是有界的，例如膨胀波（或微弱压缩），界线是马赫波，流动参数本身在波上是连续的，但它的导数在波上可以不连续。如图1所以，定常超声速气流流过外钝角，在第一道膨胀波

O L

11的上游，各流动参数都是均一的，对 的导数到处都是零，但一到

O L上便开始

11

变化了，虽然流动参数本身在

O L还是连续的，但无变化的直匀流区突然在这条

11

线上有变化的扇形膨胀地带相接，诸流动参数的导数在

O L上必是突然从零变为

11

某一定值。在最后一道波

O L上，导数从一定值跃变为零。在中间各道波上，流

12

动参数的导数取了特殊的突跃值——零。

图1 外钝角绕流

流动参数的导数在其上发生突跃的线（或面）称为弱突跃线（或面），以区别与流动参数本身发生突跃的强突跃面（如激波）。马赫线（膨胀波或微弱压缩波）就是弱突跃线。当然，突跃之中可以包括突跃值为零（即不发生突跃）的情况。这样就可以把绕外钝角流动的扇形地带中的每一道膨胀波都包括在内，都是弱突跃线。

上述几种流动的运动方程都是双曲型二阶偏微分方程。解双曲型二阶偏微分方程时，存在有一类特殊的曲线，即函数的二阶导数在这种曲线上所以不连续，因而知道了因数在这种曲线一侧的数值时，不可能靠级数展开去推得函数在曲线另一侧的数值来。双曲型二阶偏微分方程所独有的这种函数导数可以不连续的曲线，在数学上称为特征线。流场上的弱突跃线就是数学上的特征线。根据数学上的特征线理论，用数值解法解决定常超声速流动的问题，就是气体动力学中的所谓特征线法。

特征线理论除了定义特征线、并根据定义找出特征线方程之外，还确定；沿特征线上的各参数的变化之间必须遵守一定的关系，这是特征线所具有的另一方面的性质。所以说，特征线有两种性质：一是跨过这种曲线，函数(例如位函数)的二阶导数可以不连续；二是沿着这种曲线，函数的一阶导数(例如流速)的变化之间又必须遵守一定的规律。就是说，既有突跃的一面、又有规律性的一面。在做计算时，要同时利用这两种性质，尤其是利用沿特征线的规律性。这个规律性在绕外钝角流动的问题中已经用过了，那就是气流折角与流速有一一对应的关系。不过，这个关系在绕外钝角流动的问题里，没有从数学上强调它是沿特征线的关

系，表面上看来反而像是跨特征线似的。原因是，在那个简单问题上，没有全面讨论流场，只是画了要用的马赫线。事实上，与流线成μ角的马赫线还有另一族存在，如图2上的PQR 线。后面可以证明，在这个具体问题上，所谓沿特征线的变化关系正是沿PQR 那一族马赫线的变化关系。

图2 外钝角绕流的马赫线

本节先把定常平面超声速流和轴对称超声速流的运动方程列一下，以说明数学上的共同点，然后导出这些方程列的特征线来。接下去详细讲解特征线法在乎面无旋定常超声速流问题上的应用及解题法。后面再讲超声流流过圆锥体的稻确解及轴对称特征线法。

1.2 两个自变数的运动方程

平面不可压位流的运动方程是

22220x y

φφ

??+=?? 而在平面定常无旋可压流的条件下，此式府写为如下形式

2

2222222221111120c x x c y y

c x y x y φφφφφφφ

?????????????????-?+-?-???

=???? ? ? ??????????????????????? 通乘以2c ，并引用符号

x x φφ?=?，y y φ

φ?=?，22xx x φφ?=?，22yy y φφ?=?，2xy x y

φφ?=??

则上式化为

()()2

22220x xx x y xy y yy c

c φφφφφφφ-?-??+-?= （1-1）

定常无旋轴对称流的运动方程直接写为

()()2

2

22

220x

xx x r xr r

rr r c c c r φφφφφφφφ-?-??+-?+= （1-2）

这两个方程是同一类型的，都是具有两个独立变数的二阶非线性偏微分方程。它们有时也称为拟线性方程，因为对最高阶导数而言是线性的。此二式可写成一个共同的形式

20xx xy yy A B C D φφφ+++= （1-3）

一般说来，式中的系数A ，B ，C ，D 是x ，y ，x φ，y φ的函数。这样一个二阶偏微分方程，究竟属于哪种类型，要看判别式()2B AC -大于、等于、还是小于零。现在式（1-1）中的

()()()()22222222222

222x y x y x y B AC c c c c c V c φφφφφφ??-=?--?-=+-=-??

式（1-2）中的

()()22222

222x y B AC c c c V c φφ??-=+-=-??

在超声速流(V c >)中，这两个判别式都大于零，所以式（1-1）与（1-2）都是双曲型的方程。现在总概括成式（1-3）来研究，目的是要得到解决实际问题的具体办法。不过，并不能单纯从数学角度来研究式（1-3）的解法，而是要与气体流动的物理情况密切结合起来。函数φ就是速度位，x φ就是x 分速，y φ就是y 分速。也常把x φ和y φ组成的平面称为“速度面”。

1.3 哥西问题及特征线方程

在一定的边界条件下直接积分式（1-3）是很困难的，所以就希望用数值解的办法求出需要的答案。以平面流动来说，假设知道了待求函数在某一条曲线上的函数值，如果能够设法求出函数在该曲线附近的数值，就能把全流场上的解答一步步地找出来，这在数学上叫做哥西问题。哥西问题的数学提法是这样：给定函数的偏微分方程，并在xy 平面上沿某一条曲线L 给定该函数的数值，问能否推出L 曲线附近的函数值。

在所讨论的问题中，未知函数是φ。此处限于讨论介质属性(各流动参数)是连续的流场。这样，x φ和y φ(诸分速)就是x 和y 的连续函数。如果从曲线L 出发，往附近任意走一小步，x φ和y φ的改变量应为

x x x xx xy d dx dy dx dy x y

φφ

φφφ??=

+=?+??? （1-4） y y y xy yy d dx dy dx dy x

y

φφφφφ??=

+

=?+??? （1-5）

因此，如果能够以曲线L 为基地逐步向外开拓，位函数φ除了应该满足式（1-3）之外，还应满足式（1-4）及式（1-5）。而要做到这一点，各个二阶偏导数的数值必须是能够确定的。难道说有不确定的情况吗？有的，下面就谈这个问题。

从几何上想一下，微分方程式（1-3）的每一个解可以看作是x ，y ，φ空间中的一个三维曲面，叫做积分曲面，而且每个解都定义了一个满足式（1-3）的函数(),x y φφ=。在积分曲面上可以画出很多条不同的空间曲线，每一条曲线在

xy 平面上都有一定的投影。在指定点处，从这些被投影的曲线向外每走一小步，

就可按式（1-4）与式（1-5）确定出相应的x d φ与y d φ值。注意，积分曲面上的任意一点处，式（1-3）都是能满足的。此外，式（1-4）及式（1-5）是适用于位于该曲面上的任何曲线的无限小微段所对应的改变量。因此，位函数的二阶导数应同时满足下列三式

2xx xy yy A B C D φφφ++=-

0xx yy x dx dy d φφφ?+?+= （1-6） 0xy yy y dx dy d φφφ+?+?= （1-7）

这一组方程可以看作是xx φ，xy φ及yy φ的线性代数方程，解得

22

002200

x

y x y xy A D C dx d d dy A dy d C dx d D dx dy

A B C A dy B dxdy C dx dx dy dx

dy

φφφφφ-??+??+??=

=

?-?+? （1-8）

xx φ以及yy φ的形式与xy φ是类似的，分母行列式都一样，只是分子行列式有所不同罢了。

由式（1-8）可以看出，如果分母行列式等于零，那么，诸二阶偏导数就不能确定，也就无法根据给定曲线上的函数值来推算在曲线附近的函数值。这样一类具有特殊性质的曲线，它们虽然是积分曲面上的线条，但是x φ及y φ的导数在这些线上可以不连续，这类曲线，如果存在的话，称为解的特征曲线（或简称为特征线），特征线在xy 平面上的投影，称为物理面特征线。注意到x φ及y φ就是u 和v ，位函数的二阶导数不连续就意味着速度的导数不连续，因而所有流动参数的导数也都不连续。这正符合在本节开头处说的弱突跃的定义。所以特征线的物理意义很清楚，就是弱突跃线。

令式（1-8）的分母等于零，即得特征线在物理平面上的投影的微分方程，是

2

20dy dy A B C dx dx ????

-+= ? ?????特特

（1-9） 由此，得

dy dx ??

=

???特

（1-10） 式（1-10）就是物理面特征线的微分方程，它规定厂物理特征线斜率的变化规律。该式只有当20B AC ->时才有意义。可见，只有双曲型的偏微分方程才具有特征线。椭圆型方程，由于20B AC ->而不存在特征线。再看，该式中有正负号、表示在同一点可以有两个斜率，也就是说在xy 平面上有两族特征线存在。负号定为第Ⅰ族特征线，正号定为第Ⅱ族特征线。而在式（1-10）中的3个系数A ，

B ，

C 是与x φ与y φ有关系的。为了作出物理面特征线，就必须确定x φ及y φ沿着

特征线是怎样变化的。

1.4 函数的导数沿特征线的变化

在特征线上，既然分母行列式等于零，那么，在用式（1-3）、式（1-6）及式（1-7）求φ的3个二阶偏导数时，分子行列式也必须为零，否则会出现无限大的答案；而在物理问题里，参数总是有限值。分子行列式为零便规定了x φ与y φ的变化之间有一定的关系。以求xy φ这个二阶偏导数约分子行列式来说（xx φ和yy φ也一样），它是式（1-8）的分子，即

0x y A dy d C dx d D dx dy φφ??+??+??= （1-11）

由此解得

y x x d A dy D dy d C dx C d φφφ????

??=-- ? ? ???????特特特

（1-12） 把式（1-10）代入式（1-12），得导数沿特征线的变化规律（也称为“相容性条件”）为

y x x d D dy d C d φφφ????

=- ? ?????特特

（1-13）

式（1-13）给定了x y φφ-平面上的特征线斜率与x φ，y φ，x ，y 之间的关系。因x u φ=，y v φ=，故又把式（1-13）称为速度面上的特征线方程。与物理面特征线相类似，速度面特征线也有两族，第Ⅰ族取负号，第Ⅱ族取正号。

为了明确起见，常把物理面和速度面上的两族特征线方程分别写为

dy dx ??=

???Ⅰ

（1-14a ）

dy dx ??

=

???Ⅱ

（1-14b ）

y x x d D dy d C d φφφ????

= ? ?????Ⅰ

Ⅰ （1-15a ）

y x x d D dy d C d φφφ????

= ? ?????Ⅱ

Ⅱ （1-15b ）

需要强调一下，式（1-10）虽然是根据函数的导数发生突跃的条件导出来的，然而决不是说在每一个具体的流动问题里都必定发生弱突跃，只是说凡发生弱突跃必在特征线上罢了。同理，式（1-13）也只是规定了沿特征线的变化规律，并没有说沿每一条特征线都必须有变化，只是说沿着特征线如果有变化，其变化规律必是式（1-13）。在具体问题中．究竞在四条特征线上发生弱突跃，究竞沿哪条特征线仑变化，那是由具体的边界条件所规定的。在绕外钝角流动的例子中，第Ⅰ族特征线是发生弱突跃的线，但沿着第Ⅰ族特征线没有变化；第Ⅱ族特征线是参数沿之起变化的线，但不是发生弱突跃的线。

1.5 利用特征线关系式作数值解

可以利用式（1-14）及（1-15），根据给定的导数值进行数值解去求其他点的导数值。因为位函数的导数就是流速，正是所要求的参数。具体求解的时候，耍把式（1-14）及式（1-15）中的微分用差分代替，并根据上游的已知数据结合边界条件一步步作下去。参看图3，假定有十分靠近的、且不在同一条特征线上的两点1N ，2N ，并知道了这两点上的一切流动参数。这时用式（1-14）可以算出过1N ，2N 两点的特征线方向来。实际存在的特征线一般是曲线，但当1N ，2N 十分靠近时，过1N 的第Ⅱ族特征线与过2N 的第Ⅰ族特征线的交点也不会离1N ，2N 太远，

做初步近似，可以用宜线段代替曲线。这样只要根据算出来的斜率值，并从1N 和2N 分别作两直线段，便得一个交点11N ，这个新点的几何伦置很容易确定，即

111,1

111y y dy x x dx -??

= ?-??Ⅱ （1-16）

112,2

112y y dy x x dx -??

= ?-??Ⅰ （1-17）

,1dy dx ?? ?

??Ⅱ和,2

dy dx ?? ???Ⅰ分别是过2N 点的第Ⅱ族特征线斜率与过2N 点的第Ⅰ族特征线

斜率。解此二代数方程，得坐标()1111,x y 。导数值()11x φ和()11y φ可用式（1-17）求出。因1N 至11N 是沿第Ⅱ族特征线变化的，故按式（1-15b ）得

()()()()()111111111111

y y x x D y y C φφφφ?-=---? （1-18） 同理，因2N 至11N 是按第Ⅰ族特征线变化的，故按式（1-15a ）有关系式

()()

()()()2

11211

211

2

22

y y x x D y y C φφφφ?-=---? （1-19） 脚注1，2，11分别表示1N ，2N ，11N 点的数值。联解式（1-18）及式（1-19），即得11N 点的()11x φ及()11y φ之值。

图3 特征线网示意图

若需要精确一些的话，可在求出初步近似的()11x φ及()11y φ值之后，把11N 点的A ，B ，C ，D 值算出来，然后用1N ，2N 与11N 的平均A ，B ，C ，D 值再重列式（1-18）及式（1-19），解出更精确一步的()11x φ及()11y φ值来。

如果原给的不止1N ，2N 两点，还有3N ，4N ，…一系列不在同一条特征线上的点，则重复上述步骤，可得12N ，13N ，…等新点。再以11N ，12N ，13N ，…为新的起点，即又可求出21N ，22N ，23N ，…等点来。这样就可以根据给定的情况算出下游一定区域中的流动情况或反推出上游一定区域中的流动情况来。

如果已知数据的点正好给在一条特征线上，就无法推算新的点。这时需要附加其他条件一起进行计算，这在后面再讲。

2 平面无旋流的特征线法

2.1 物理面与速度面上的特征线方程

参看式（1-1），平面定常无旋流的

2222x A c c u φ=-=-， x y B u v φφ=-?=-?

22

22y C c c v φ=-=-， 0D =

把这些值代人式（1-14a ）及式（1-14b ），得物理面上两族特征线方程为

dy dx ??

== ?

??Ⅰ,Ⅱ

（1-20） 式（1-20）表示的是特征线和x 轴之间夹角的正切关系式。如果算一下特征线与速度矢V 量之间的夹角的正切，即可证明特征线就是马赫线。因为特征线与

x 轴之间的夹角是()arctan dy dx ，速度矢量V 与x 轴之间的夹角是()arctan v ，二者之差的正切是

2

tan 1

1dy v dx u u v dy Ma u dx

-

==

=-+

（1-21）

该式说明特征线与速度矢量之间的夹角为马赫角，所以特征线就是马赫线。下方符号是第Ⅰ族特征线，相应地，在式（1-20）中取下方符号；上方符号是第二族。参看图4，顺着V 正的指向看过去，第Ⅰ族特征线向左伸，第Ⅱ族则向右伸。

图4 特征线与流线的关系

如记V 与x 轴之间的夹角为δ，则把式（1-20）作适当变换或直接按图4写下物理特征线方程

()tan dy u dx δ??

=+ ???Ⅰ,Ⅱ

（1-22）

再把A ，B ，C ，D 代人式（1-15a ）及式（1-15b ），得速度面特征线方程（即位函数的导数沿特征线的变化关系）

2222

dv B uv V c du C c v ±-??

=== ?-??Ⅰ,Ⅱ

1-23） 若用合速V 及角度δ来表示上式，形式会简洁得多。换算如下

222V u v =+，arctan v u δ??

= ???

微分得

V dV udu vdv ?=+

2

udv vdu

d V

δ-=

合并此2式，得

dv u v dV du v Vd u v u δ??+ ???=??

- ???

将此式与式（1-23）联立消去v u ??

???

，最后得

2

2

tan dV

u Vd V c

δ

==-

在均能流场中，cr c 是常数，将上式分子分母通除以cr c 后得

tan d d λ

μδλ

=? （1-24）

此式又称为“相容性条件”。式巾的正号是沿第Ⅰ族特征线变化的关系式，相应地在式（1-23）中取下方符号。上方符号则是第Ⅱ族。该式规定：沿特征线，δ如有变化，λ必有相应的变化；当然也可以不变。

总之，两族物理面特征线及速度面特征线方程是

()tan tan dy dx d d δμλμδλ???

=+ ??????

???=? ?????Ⅰ

Ⅰ

（1-25） ()tan tan dy dx d d δμλμδλ???

=- ?????????=-? ?????

ⅡⅡ （1-26） 2.2 平面无旋流的特征线网图

注意式（1-24），该式中不包含物理面上的坐标（即无dx 项）。因此，可以

直接积分。因为tan μ也是λ的函数，代

tan μ=

入式（1-24），得

2

2d d λδλ=

正负号的上下次序换了一下，把第Ⅰ族特征线的符号写在上面了。积分得

C δ=±+

?

（1-27）若规定1

λ=时，气流方向0

δ=，则0

C=。

2.3 四种边值问题

图5 第一边值问题

第一边值问题：已知定常超声速流场上某一条普通曲线AB上所有各点的坐标及λ值和δ值，求以AB线及过A，B两点的两条特征线AC和BC所围区域内的流动情况。

注意给定的只是AB曲线，特征线AC和BC是未知的，等问题解完了才能确

定下来。解法是这样的：参看图5，在物理面的给定曲线AB上选取一系列的点

1

P，2

P，…，

n

P，间距一艇说来应很小，具体取决于所需的精确度。各点的坐标(),x y 及λ和δ值皆已知，具体执行计算时采用如下步骤。

引用记号

()

tan

Aδμ

=+

()

tan

Bδμ

=-

tan

Dμ

=

tan

Eμ

=-

将特征线方程和相容件条件式（1-25）和(7—N)改写成以下形式

y A x ?=?? 沿C + y B x ?=?? 沿C -

D λλδ?=??? 沿C +

E λλδ?=??? 沿C -

此处，C +表示第Ⅰ族特征线，C -表示第Ⅱ族特征线。在图5止，假设1AQ ，12PQ ，

23P Q 等是第Ⅱ族特征线C -；11PQ ，22P Q ，33PQ 等是第Ⅰ族特征线C +。因此，从1P 和2P 点出发求2Q 点时，可得下列方程

()

2121Q P Q P y y B x x -=- ()

2222Q P Q P y y A x x -=-

()2

1

1

2

1

Q P P Q P E λλλδδ-=-

()2

2

2

2

2

Q P P Q P D λλλδδ-=-

式中，A 和D 是指1P 与2Q 点的A 和D 的平均值；而B 和E 则是指2P 与2Q 点的B 和E 的平均值。由此解得

()(

)2

121

2

P P P P Q y y Ax Bx x B A

---=

- （1-28a ）

()

2121Q P Q P y y B x x =+- （1-28b ）

()(

)2

12211

2

1

2

P P P P P P Q

P P D E E D

λλλδλδδλλ---=

- （1-28c ）

()2

1

1

2

1

Q P P Q P E λλλδδ=+- （1-28d ）

其中的

()1212P Q A A A =

+，()

221

2P Q B B B =+ ()1212P Q D D D =+，()

221

2P Q E E E =+

计算过程要作迭代。开始时，为了求平均值A ，B ，D ，E ，需要估计2Q 处的流动参数λ和δ，然后按式（1-28）算出点2Q 处的参数。得到新值2Q λ和2Q δ以

后，即可计算新的平均值A ，B ，D ，E 。按此步骤作迭代，直到点2Q 处的第()1n +次的1n λ+值与第n 次的n λ值之差小于一个规定值（例如610n λ-)，就认为迭代收敛了。更多的细节不在此讨论了。

然后从A 起，每两个邻点决定一个新点的位置及流动参数。例如A 和1P 决定

1Q 。1Q 的几何位置由过A 的第Ⅱ族特征线及过1P 的第Ⅰ族特征线的交点所决定。点2Q 的x ，y ，λ，δ是由点1P 和2P 处的参数确定。依此类推。这样一步步做下去可以把1Q ，2Q ，…，n Q 各点全定出来。同时也定出了过A 的第Ⅱ族特征线上的第一个小段1AQ ，及过B 的第Ⅰ族特征线上的第一个小段n BQ 。依此类推，终于可以把曲线三角形ACB 中的流动情况都定出来。C 点及AC 和BC 曲线也在迫近过程中定出。

第二边值问题：在物理面上给定过某个A 点的两条特征曲线AB 和AC ，及沿线的流动参数λ和δ，求曲线四边形ABCD 区域内的流动情况。

参看图6，本问题中的BD 及CD 线是待画的。

图6 第二边值问题

解法是在AB 上取分点1P ，2P ，…，在AC 上取分点1Q ，2Q ，…。1P 和1Q 不是同一条特征线上的点，通过尸1P 和1Q 作两条不同族的特征线交于1R 点，按第一边值问题的解法，从1P 和1Q 可以定出1R 来。1R 再和2P 或2Q 去决定下一步的点。一步步做下去，总可以把ABCD 区域中的流动情况处处都定出来。

第三边值问题：给定流场上的一段特征线AB 及线上各点的流动情况λ和δ，

另给定一段过A 点的固壁AD ，求曲线三角形ABC 内的流动情况。

参看图7，BC 是待定的特征线。

图7 第三边值问题

解法是先把A ，B 两点及其间各分点1P ，2P ，…的λ和δ定出来。过1P 作第Ⅱ族特征线，交固壁于1Q 点，1Q 的流速指向必与该点处固壁的切线一致。在固体边界上，流速大小是其指向决定的。往下都是类推的做法，不必多说。

第四边值问题：给定流场上—条特征线AB ，及AB 上各点的流动情况λ和δ，且A 点在自由边界上，自由边界以外是已知压强为a p 值的静止气体，求曲线二角形ABC 内的流动情况。参看图8，AC 是自由边界，BC 是特征线，都未知，有待做出。

图8 第四边值问题

与前一问题不同处在于，自由边界广气流的压强必须等于外界的a p 。这个压强规定了自由边界上的流速值，然后流速规定指向，从而得出边界流线AC 来。

解法是先按外界压强算出边界上的λ值：()0a p p πλ=求得λ值，其余步骤与第三边值问题相同。