[讨论]伯努利方程是能量方程还是动量方程？

coolboy

顺便说一下偏微分方程边界条件的设置问题。我在上面的[285楼]的总结中提到：

++++++++++++++
流体的总能量沿速度场是守恒的，即可以写成沿速度场或即“特征线”的全导数。既然是全导数，则就可以对之进行积分，即沿流线积分，所得积分表达式（常数）就是伯努利方程。
++++++++++++++

这里所谓的“沿特征线进行积分”在数学上仅仅是针对双曲型方程或也即波动方程而讲的。也就是说，所谓的“沿流线积分”也即在数学上限定了或定义了推导伯努利方程必须是求解一个双曲型方程。因为双曲型方程本质上可以或必须通过“沿特征线积分”而进行求解，其边界条件则也必然由（每根）特征线或也即流线上的一点所给定。例如在本主题帖一开始的[5楼]，通流版主出了一个关于伯努利方程的应用题：

++++++++++++++
通流：一个水管连着比水管出口高十米的水箱。水管的出口连着大气，水箱也跟大气相连。计算水管出口的水的流速。
++++++++++++++

这里的水管就相当于流线或特征线。知道了水管入口处的边界条件，则特征线上其它各点（包括出口处）的值都可通过沿特征线积分而把它求出来。这里的伯努利方程则刚好对应于（全）导数为零，从而沿特征线积分刚好是一常数的特例，即流管各处的总能量是一常数。

说到这里，细心的读者应该理解关于如何设置流场入口、出口处边界条件这一问题了。从数学物理（方程）的角度来讲，首先应该搞清楚所对应的方程究竟是双曲型的波动方程，还是抛物型（或椭圆形）的热传导、扩散方程。对于前者，我们只能给出入口或出口处的一个条件。若同时给出两边界条件，则理论上应该导致问题的矛盾无解。对于后者，则我们必须具备两个边界条件才能数学求解。

这里还能看到；对于给定的某个流场，边界条件的设置同描述流场物理过程的数学模式直接相关。无黏流体的数学模式通常对应于双曲型方程。若在方程中加入黏性项（实际黏性很小，但总是存在的）则应该意识到这一过程更重要的是改变了方程的特性及相关的边界条件的设置。


[ 本帖最后由 coolboy 于 2013-11-1 11:38 编辑 ]

coolboy

推导伯努利方程的步骤是先把动量方程转化为能量方程，然后沿流线（即特征线）对能量方程积分就得到了伯努利方程。象上面通流版主给出的例题中，通过应用伯努利方程，我们可得到沿流线包括出口处的各处的速度值。尤其重要的是我们可以很明显地看出，出口处的速度值已由入口处的条件给定。我们不可能、不应该在给定入口处的条件的前提下，又再加上出口处的速度（条件）来求解沿流线各处的速度值。

我们也可以对动量方程直接沿流线积分或直接应用某数值计算软件而求得各处的速度场，结果应该是一样的。当然，这时我们并不能从理论上推出伯努利方程。此外，若对“数学物理方程”这一领域的基础知识不熟的话，可能还会错误地设置了入口处和出口处两端的边界条件之后才会考虑如何数值求解动量方程而获得各处的速度值等等。

coolboy

下面来介绍一下如何数值求解类似于伯努利方程的所谓偏微分方程中的平流（或也有称之为对流）方程。下一个帖子中的一个回帖介绍了这一问题的背景知识：

扩散方程的C-N格式计算问题
http://www.cfluid.com/bbs/viewthread.php?tid=141147

+++++++++++++++++++++++++++++

原帖由 lwd1981 于 2014-2-21 20:28 发表

                               
登录/注册后可看大图

你可以将你的数值离散格式对应的修正方程导出来，你会发现你的数值耗散项也是负的；如果你的数值耗散是正的，那么出现间断附近（数值突变的地方）的振荡将被光滑掉（这其实就是“人工粘性”的机制），从而使得你的“矩形波”会变光滑；但是你的数值耗散项是负的，它不会将会将间断附近（数值突变的地方）所出现的振荡光滑，反而会将间断附近（数值突变的地方）所出现的振荡放大。

Coolboy(2014-2-22):

这一段的叙述表明发言者对问题的数学物理内涵有着比较清晰的认识，即数值格式通常在解的“间断”处会产生虚假波以及稳定的数值格式会出现“人工粘性”来耗散掉（物理的或虚假的）数值波动。

现在我来提出一个问题，看看是否有人知道答案或知道的是什么样的答案。假如实际的物理问题确实存在类似于上述“矩形波”即具有间断的特性但方程中并不存在物理的“粘性耗散”过程，则该如何设计数值格式使得守恒的“矩形波”特性在经过长时间或长距离积分之后仍然能被很好地保持着。

至于实际的物理过程，例如下帖子中提到的守恒的总能量若初始时刻刚好是一个矩形波，则它就可以被实际流场传输很远或很久而大致不变形，至少理论上是如此。现在的问题是数值求解此问题时该如何实现？

++++++++
[讨论]伯努利方程是能量方程还是动量方程？
http://www.cfluid.com/bbs/viewthread.php?tid=114265
++++++++

数学家冯康后期有一出名的工作就是数值求解哈密顿系统并保持其保守性。但他处理的问题是常微分方程（组）而我们现在的问题是要数值求解双曲型的偏微分方程（组）并使得解的“数值粘性耗散”尽可能地小。
+++++++++++++++++++++++++++++

coolboy

开始设计平流方程（即双曲性偏微分方程）数值格式的第一步是先看一看该方程能否被改写成“通量”型式。对于目前的伯努利方程来讲，这一点刚好能自然满足。这在本主题帖前面的[188楼]就已经具体说明：

[讨论]伯努利方程是能量方程还是动量方程？ [188楼]
http://www.cfluid.com/bbs/viewth ... age%3D1&page=13

Coolboy(2012-6-11):
+++++++++++++++++++++++++++++
说到这方法三，必然要说到流体力学中的两种不同形式的守恒律。若仅是泛泛而谈地说质量守恒、能量守恒等常常会引起误解。我在本主题帖一开始的讨论中在[8楼]还刚好说到了这两种不同形式的守恒律：

http://www.cfluid.com/bbs/viewth ... page%3D1&page=1
++++++++
守恒律一般可写成两种不同的形式：（1）全导数型，（2）通量型。伯努里方程及其推导和上面[6楼]所有的讨论是建立在全导数型的表达方式上的。这时就有“沿流线的积分”及“沿流线为常数”等的说法。若采用通量型的守恒律，则我们就有（固定）流入流出面（口）的差值同体积内部源汇的守恒关系等的描述。
++++++++

所谓伯努里方程全导数型的守恒方程就是指如下的能量方程：
{1}    @H/@t+u(@H/@x)+v(@H/@y)+w(@H/@z)=0
其中的H是流体单位质量的总能量。在巴切勒的书中，其（3.5.4）式的总能量包含4项，即流体的（1）动能，（2）热能（内能），（3）压力位能和（4）（重力）位能。

流体的质量一般并不满足全导数型的守恒律。它满足通量型的守恒方程，即
{2}    @[rho]/@t+@(u[rho])/@x+@(v[rho])/@y+@(w[rho])/@z=0
或
{2}    @[rho]/@t+[DEL]（[v_bold][rho]）=0

方程{1}两边乘[rho]，方程{2}两边乘H，结果相加，则我们可得
{3}    @([rho]H)/@t+[DEL]（[v_bold][rho]H）=0

即流体的总能量在满足伯努里方程全导数型守恒律的同时，也满足通量型守恒律。但注意到全导数型守恒律指的是单位质量的总能量，而通量型守恒律指的是单位体积的总能量。我们不能仅凭也听说了众所周知的“能量守恒”等等什么的就以为掌握了所有的物理。具体的物理问题需要具体的分析，这过程中数学又是必不可少的工具。
+++++++++++++++++++++++++++++



[ 本帖最后由 coolboy 于 2014-2-24 09:00 编辑 ]

coolboy

我们现在需要数值求解的就是如下的通量型平流方程：

@([rho]H)/@t+[DEL]（[v_bold][rho]H）= S

上面右边的S表示外源项。这一项通常不是设计数值格式的重点。通量型的方程究竟好在什么地方呢？它好就好在通量型方程中所表达的量一般是指单位体积内的物理量。在这里则刚好是单位体积内的总能量。在空间划分了固定的网格之后，这“单位体积”也就完全确定而不再变化了。对于三维流场（方程），那三维空间的网格就是一个个小立方体。对于一维流场（方程），那一维空间的网格就是一段段小间隔。每一个小立方体有六个面，每一个面上都有对应的通量值（[v_bold][rho]H）。若我们能把每一个面上的通量值都能非常精确或精细地计算出来，则流入那小立方体的净通量（幅合）加上右端的外源项就给出了小立方体内单位能量随时间的增量。

至此，整个问题就归结为如何非常精确或精细地计算出流经各界面的通量值（[v_bold][rho]H）。


（待续）

lwd1981

采用守恒形式的方程和非守恒形式的方程各有各的优点，不可一概而论。一般流场中包含间断的时候有必要采用守恒形式，否则容易得到非物理的激波位置。
    对于比较光滑的流场，不必要一定要采用守恒形式，很多时候非守恒形式离散更方便。

coolboy

原帖由 lwd1981 于 2014-2-24 12:05 发表

                               
登录/注册后可看大图

采用守恒形式的方程和非守恒形式的方程各有各的优点，不可一概而论。一般流场中包含间断的时候有必要采用守恒形式，否则容易得到非物理的激波位置。
    对于比较光滑的流场，不必要一定要采用守恒形式，很多时候非守恒形式离散更方便。

若只从“离散更方便”的角度来讲，则一阶和二阶导数分别可表示为(u1-u0)/dx和(u2-2*u1+u0)/(dx*dx)等等，把它们代入任何微分方程里就可以数值求解了？但我在上面的[288楼]中的背景介绍中所提的是如下的问题：

+++++
现在我来提出一个问题，看看是否有人知道答案或知道的是什么样的答案。假如实际的物理问题确实存在类似于上述“矩形波”即具有间断的特性但方程中并不存在物理的“粘性耗散”过程，则该如何设计数值格式使得守恒的“矩形波”特性在经过长 时间或长距离积分之后仍然能被很好地保持着。
+++++

若谁要是真的知道这一问题的解答（或具体参考文献），甚至是比我预期将要给出的解答更简单的方法，不妨在此（在我给出最后解答之前）与大家共享一下。

lwd1981

原帖由 coolboy 于 2014-2-24 04:51 发表

                               
登录/注册后可看大图

若只从“离散更方便”的角度来讲，则一阶和二阶导数分别可表示为(u1-u0)/dx和(u2-2*u1+u0)/(dx*dx)等等，把它们代入任何微分方程里就可以数值求解了？但我在上面的[288楼]中的背景介绍中所提的是如下的问题：

...

关于守恒形式和非守恒形式孰优孰劣的，CFD早有定论。你的问题无非就是一个方波的纯对流输运的问题如何高精度数值求解。建议你看看有关90年代abe的CIP的文章，专门讨论这个问题，现在能够解决这个问题的数值方法太多了。
    当然，如果你有高见，愿闻其详！

coolboy

还是希望谁能共享几篇具体的参考文献。

lwd1981

http://www.mech.titech.ac.jp/~ryuutai/publist.html
Yabe 先生的Laboratory
Constrained Interpolation Profile Method
就是针对这个问题的

coolboy

原帖由 lwd1981 于 2014-2-24 13:21 发表

                               
登录/注册后可看大图

http://www.mech.titech.ac.jp/~ryuutai/publist.html
Yabe 先生的Laboratory
Constrained Interpolation Profile Method
就是针对这个问题的

就在网上很快大致看了一下这CIP方法，同我要介绍的方法应该是同源（看来也是同精度）的，都是多项式插值提高精度。区别在于是在格点上还是对立方体内插值而已。通量型的物理意义要清楚一点。最重要的区别当然是时间的先后了。我所要介绍的方法的第一篇经典论文发表于1981年。对它进行重大改进并给出原码而在我们领域内出名而成熟的第二篇经典论文发表于1986年。对它作进一步改进的论文发表于1990年，作者还是周光召以前的学生。这都早于你所说的“90年代”。

coolboy

原帖由 coolboy 于 2014-2-24 09:02 发表

                               
登录/注册后可看大图

至此，整个问题就归结为如何非常精确或精细地计算出流经各界面的通量值（[v_bold][rho]H）。

不失一般性，我们可以假定速度场是由立方体A指向立方体B。我们现在就要求出在A-B界面处的通量值。最简单的近似计算就是假定流经A-B界面的通量就由界面处速度与立方体A内的总能量值的乘积来表示。这其实就是所谓的一阶上游差分格式，其精度是很差的。在稍微好一点的二阶差分格式中，所需的通量就由A-B界面处的速度与两立方体A和B的总能量的平均值的乘积来表示。而所谓的更高阶（例如四阶）的差分格式则是采用临近更多的立方体的能量作加权平均来计算流经A-B界面的通量。

所谓一阶、二阶、四阶差分格式从微分算子到差分算子近似的角度来讲，那是对函数用泰勒级数展开到底是取多少项截断近似来说的。但其实也可以从另一个角度来讲，即各立方体之间的值是否有关系及是怎样的关系。一阶上游差分格式其实假定A和B格点的值是相互独立的。二阶差分格式所隐含的假定则是格点A到B之间的值随距离作线性变化，等等。

尽管高阶差分格式的精度要高一点，但这类根据临近点函数值进行多项式拟合的差分格式都有一些不可避免的通病，即对于类似于矩形波的有较强空间梯度的物理量（能量、痕量化学元素、雾霾浓度，等等）长时间或长距离积分之后不能保持其该有的强梯度特性，也会出现不该有的震荡或数值虚假波、不该是负值的负能量或负雾霾浓度，等等。

下面所介绍的方法则是不同于“根据临近点函数值进行多项式拟合”的一个新思路。


（待续）


注：根据上面[293楼]lwd1981回帖l所言：“现在能够解决这个问题的数值方法太多了。”其实要是明白方法的具体思路的话，大家也可在此介绍、交流一下那些方法的基本思路。

雪饮

专家啊。。。。

coolboy

原帖由 雪饮 于 2014-8-22 21:27 发表

                               
登录/注册后可看大图

专家啊。。。。

是说我吗？你是根据我帖子中的科学内容才这么认为的吗？
好象前面回帖中就有人认为我是在读硕士生，而他及他的同事们则已经硕士毕业后工作几年了，并进而以此一致认为是我错了及还不愿认错，呵呵。

coolboy

+++++++++++++
菜鸟请教特征线问题
http://www.cfluid.com/thread-143194-2-1.html
+++++++++++++

有同学在上一个帖子中提出了一个关于“数学物理方程里面的特征线对应的物理含义是什么”的问题。该问题与本主题帖所叙述的内容密切相关。我在该帖子的[44楼]作了如下的回答。现转载如下：