2.4 瞬变流的特征线计算方法

从数学的角度来讲，特征线法在求解偏微分方程中应用非常广泛。自从文献［26］对水锤特征线计算方法进行了比较系统的阐述以后，水锤特征线的求解方法在瞬变流模拟和预测中得到了非常广泛的应用。即使如今，特征线法依旧是模拟水锤的最基本的方法。考虑到本书中的计算方法主要采用了特征线法，因此对该方法进行介绍。特别说明，该方法在文献［38］以及文献［3］都已经有详细介绍，读者若想更详细地了解可以参考这两本文献。

水锤特征线计算方法的基本思想是：通过将连续方程和运动方程进行等价变换，选择特征辅助线将偏微分方程转换成与原方程有相同解的常微分方程组；通过沿特征辅助线路径积分，将常微分方程转换成积分形式的方程组，用近似法将常微分方程组改写为差分方程组，最后通过数值计算对差分方程组进行求解。特征线法原则上可以解决任何边界和复杂管道的计算，能较方便地把阻尼的影响因素考虑到计算中来；在利用电子计算机的编程计算中，它具有较高的精度和稳定性，而且编程也较为简洁。

2.4.1 水锤特征线求解方法

简化后的瞬态水流连续方程式（2.19）和运动方程式（2.24）构成了一对准线性双曲微分方程。将这两个方程分别用L₁及L₂表示为

引入一个特征因子λ加以线性组合可得复合方程

对于组合方程，只需取任意两个不等的λ值，便可以重新组合得到一组新的方程组，而且这个方程组与原来的偏微分方程具有等价解。在原来的方程组中v和h是x和t的函数，如果在新的方程组中假设x是t的函数，根据微分方程可得

则组合方程式（2.27）可变形为常微分方程，即

只要适当选择λ的特定值，就可以满足式（2.29）的假设条件，假设λ为特定的未知数，通过求解满足假设的条件的方程可得

这两个特定的λ值，即为特征值，此时对应的特征辅助方程为

这两个辅助方程指出了x和t在特征线法计算中所应满足的特定条件，它表明了水锤波位置的变化是由时间和波速a确定的。Wylie等在文献［38］中，将这两个特征线方程描述的路径成为特征线，并且分别成为正向特征线（forward characteristic line）和反向特征线（rear characteristic line）。在这两个特征辅助方程的约束条件下，将两个特征值分别代入到常微分方程（2.30）可得

经过以上变换，两个偏微分方程转换成了一组常微分方程，而每一个常微分方程都制约于一条特征线方程所表达的特定路径，这两条路径分别称为正向特征线和反向特征线。

2.4.2 特征线方程组的差分变换

特征线方程组是常微分方程，目前还无法找到解析解，通常的做法是将其转换成差分方程的形式求解。为满足特征线的约束条件，如图2.3所示，给定管道上的两点A和B，通过这两点的正向特征线CF和反向特征线CR可以交汇到未知的P点，如果沿AP积分，那么满足正向特征线相容方程式（2.33），同理沿BP积分，可以满足反向特征相容方程式（2.34）。对于给定的管道，水锤波速a可以认为是常数，因此两条特征线为直线，在图2.3中分别用CF和CR表示。

图2.3 特征线交汇

下面以adt/g=dx/g乘以正向相容方程，并把方程中的流速写成流量的形式，通过沿正向特征线积分，方程可以转换成以下形式

由于最后一项中q和x的关系不明确，在实际中用一阶近似就可以得到足够的精度，再对以上方程式积分可得

同理，沿BP方向对反向相容方程式积分可得

式（2.36）和式（2.37）描述了管内q和h在水锤传播过程中向上下游传播规律。令B=a/（gA），R=fΔx/（2gDA²）；式（2.36）和式（2.37）可简化为以下形式

式（2.36）和式（2.37）可更进一步简化为

式中 C_F、C_R——对应正向特征线和反向特征线的系数。

如果t时刻的A、B两点的流量和水头已知，那么可以通过以上方程求得t+Δt时刻P点的流量和水头为

式（2.44）和式（2.45）即为特征线方法的基本解，在实际的应用中，需要将计算区域进行离散后建立相应的计算网格。

2.4.3 离散网格划分

以上推导过程表明，当已知P点上下游相邻节点A、B在初始时刻的水力参数时，可以通过特征线法求得P点的末了时刻的水力参数。在管道的瞬变流模拟中，为了应用特征线法，需要将管道均分成长度相等的计算单元，每个单元包含进口和出口两个断面。在一根管道中，断面的数量比单元的数量多一个，如果一根连续的管道被等分成n段，则在瞬变流模拟中，这根连续管道的单元数量为n，它的断面数目为n_s=n+1，计算单元的长度为Δx=l/n，该长度也称为空间步长（step）。如果水锤波速为a，那么对应的时间增量为Δt=Δx/a，该增量也称为时间步长（time step）。根据空间步长和时间步长，可以将一根连续的管道进行时间和空间的网格划分，该网格将连续的空间和时间离散成规则分布的网格，文献［38］中称为x-t网格。在计算中，每一个节点代表了特定时间上的特定断面。如果这些节点足够稠密，那么节点上的水力参数就能够描述连续管道上的水流条件，节点的水力参数随时间变化的过程就描述了连续管道中水流的瞬变变化过程。在计算中，只需对节点的水力参数进行计算，便可以获得连续管道中水流的整个瞬变变化过程。为了简洁地说明特征线法整个计算过程，以下对网格划分进行举例说明。

如图2.4所示为被均分成4个计算单元的连续管道，该管道形成了5个断面（节点），其中1个进口断面，1个出口断面，3个中间断面。根据特征线网格划分的条件，图中展示了4层计算网格。为了使理解变得更加容易，用不同形状对网格中的节点进行了分类，图2.4中黑色圆点表示为初始已知节点或已经求解后的节点，空心圆点表示待求解的中间节点，空心半圆点表示待求解的边界节点。这些节点的水平间距即为空间步长Δx，垂直间距为时间步长Δt。假设t时刻为当前计算的已知时刻，该时刻所有节点的水力参数都已知，在此基础上通过特征线法求解下一时刻t+Δt的各个节点的水力参数。

图2.4 特征线法计算网格划分

例如，当求解t+Δt时刻，位于x处的节点时，根据式（2.44）和式（2.45）可得到该节点的解为

式中C_F（x，t）和C_R（x，t）可分别由已知时刻的上下游相邻节点求得，即

由于特征线网格的空间步长和时间步长都是均匀等距的，因此可用序列号来标示各个节点的位置。由于在任何迭代过程中，只与计算时步的前后的已知时刻t和未知t+Δt有关，因此在计算中只需要记录一个统一的时间过程来记录各层节点的时间，则以上的解可以用更简单的形式表示为

以上过程计算了在t+Δt时刻第i点的水力参数，同样的方法，可以计算所有t+Δt时刻中间节点的参数。对于边界未知节点，由于上游只有一个反向相容方程，而下游只有一个正向相容方程，因此，无法直接通过以上特征线相交的方法求解。在瞬变模拟中，边界条件可以有多种形式，但一般都可以直接确定一个参数，如水位或者流量，也可能是一个约束方程，这样结合正向相容方程和反向相容方程便可求得边界节点的解。确定边界节点的解后，t+Δt时刻的节点全部成为已知，在下一个计算时步中可以作为新的初始值求解t+2Δt时刻的解。

以此类推，可求任意节点不同时刻的瞬变水流参数。在实际的模拟中，由于瞬变流的模拟通常是由一种恒定流状态过渡到另一种恒定流状态，因此瞬变流的起始时刻的水力参数往往可通过恒定流时的状态决定。此外，比较复杂的流动也可通过瞬变流模型进行迭代求解。在实际的瞬变流模拟中，在确定初始状态的水力参数后，也就是确定0时刻各个节点的水力参数后，可以依次求得Δt，2Δt，3Δt，…的解。只要Δt足够小，就可以看成是时间连续的瞬变水流变化过程。以上案例采用的是4个单元，通过计算机可以设定更多的单元，以使空间步长Δx足够小，这样就能在空间上反映任一时刻连续管道的水流参数的沿程分布情况。

根据以上的特征线法计算过程，通常需要给定单元的数目，然后根据单元的数目确定单元的长度Δx，此后，根据波速确定时间步长Δt=Δx/a。这种方法在计算时不方便按指定的时间步长进行计算，也不方便在复杂管网中统一时间步长，以进行同步计算。为了克服这一个问题，本书在划分计算网格时，优先选定时间步长，然后根据时间步长确定单元的长度，相应的计算方法为

由于单元数必须为整数，所以需要对波速a进行微小调整，以满足正反向特征线积分路径，波速调整的准则为

其中调整后的波速为

为了简便起见，以下调整后的波速仍然用a表示。只要划分的单元长度足够小，这样调整后的波速对计算结果的影响可以忽略不计。但采用以上调整波速的方法可以方便地设计计算时间步长，对于复杂的管网系统也可以非常方便地设置相同的时间步长，从而能更加容易地实现所有系统中各个管道之间的同步计算。

2.4.4 瞬变流的基本计算流程

根据分析对象和边界条件的具体情况，瞬变流的模拟具有很强的独特性。特别是每一项调水工程的瞬变流模拟都具有不同的规模、形式和独特的边界条件。然而，瞬变流模拟的基本原理和方法却具有相似的流程。通常情况下，瞬变流的模拟包括了以下几个主要的步骤：①分析工程的形式，选取合适的瞬变模拟模型；②概化工程的布局和基本参数；③分析工程的边界条件，建立合适的理论计算模型；④数学建模，设定瞬变流模拟的框架；⑤编制数值计算程序进行数值模拟；⑥数值计算结果验证与分析；⑦提出优化的控制方案和防护措施。

数值模拟计算程序流程如图2.5所示。

图2.5 数值模拟计算程序流程

2.4.5 特征线方法验证

为了验证水锤特征线方法的正确性，Wylie和Streeter曾采用特征线法计算与往复泵相联的简单管道中的水锤压力过程，其计算结果与实测结果的比较非常接近。如图2.6所示为Wylie和Streeter的计算结果和实测结果的比较情况。

图2.6中实线是通过传感器采集的压力变化过程线，空心点是由特征线法计算的数据点。该图所示比较说明，采用水锤特征线法计算的结果与实验结果吻合，特征线法能够较好地模拟管路中的瞬变水流问题。

图2.6 往复泵单管的水锤计算与实测比较