针对多组分两相系统的UVN-flash模型的重新构建及其在管道中二氧化碳富集混合物传输的应用

《Computers & Fluids》：A reformulation of UVN-flash for multicomponent two-phase systems with application to CO2-rich mixture transport in pipelines

　　帕德雷普·库马尔 | 帕特里西奥·I·罗森·埃斯奎维尔
荷兰阿姆斯特丹数学与计算机科学中心

摘要
在碳捕获与封存（CCS）中，富含二氧化碳的密相混合物通过管道运输是一个关键环节。准确建模需要将流体动力学和热力学相结合，特别是在减压等瞬态事件期间。在这项工作中，我们提出了一个统一的框架，用于处理管道中的两相多组分运输问题，该框架同时考虑了这两个方面。具体而言，我们采用均匀平衡模型（HEM）来模拟富含二氧化碳的两相混合物的运输过程，并通过基于亥姆霍兹能量的状态方程来提供热力学封闭。相平衡计算使用UVN-flash软件进行，并通过稳定性分析程序来检测相分离并生成相平衡计算的初始猜测值。我们引入了一种新的定制UVN-flash程序，使其与流体动力学公式相匹配。这一改进是通过引入一套更适合相平衡计算的变量实现的。所提出的框架已应用于含有富含二氧化碳混合物的储罐和管道的减压过程，证明了其在CCS相关应用中的有效性。

1. 引言
根据国际能源署[1]的报告，碳捕获与封存（CCS）将在减少温室气体排放方面发挥重要作用。在CCS中，二氧化碳从工业来源捕获后输送到储存地点（例如枯竭的气田），并在其中注入地下进行永久储存。预计大部分运输将通过管道完成。来自工业的捕获二氧化碳通常不纯，可能含有不同水平的杂质，如甲烷、氮气、胺类、硫氧化物和水等[2]。这些杂质会显著影响相行为、热物理性质，从而影响管道系统的设计和安全运行。

管道操作中的一个极端事件是减压，这可能发生在关闭操作之后[3]、[4]、[5]。在这些事件中，温度和压力可能会达到极限值，有可能超过管道壁的材料运行限制。精确预测这种行为对于安全设计管道至关重要。在本文中，我们研究了含有杂质的二氧化碳通过管道的减压过程。

一些实验研究已经探讨了富含二氧化碳的混合物的减压行为，为模型验证提供了宝贵的参考数据。Drescher等人[6]提供了相关混合物的实验数据。Cosham等人[7]进行了14次冲击管实验，包括纯二氧化碳和含有杂质的二氧化碳的实验，以研究减压行为。Botros等人[8]报告了相关混合物的实验数据，以确定减压波速。

流体混合物的两相流动由纳维-斯托克斯方程描述。关于通用两相流动模型，读者可以参考Stewart等人的工作[9]；对于二氧化碳运输 specifically，可以参考Munkejord等人的综述[10]及其中的参考文献。对于管道应用，通常会进行横截面平均处理，以获得一维双流体模型[11]。这种平均处理可能会使方程组变得不适定。为了避免这些问题，并更加强调流体动力学和热力学之间的结合，我们采用了均匀平衡模型（HEM），该模型是双曲型的，能够保证方程组的适定性。在HEM中，流体被视为两种良好混合的均匀混合物，纳维-斯托克斯方程直接用混合物密度、速度和能量来表示。对于无摩擦管道中的不可压缩流动，这种表述简化为欧拉方程，其中所有守恒量都用混合物属性表示。

除了流体动力学外，实际流体的相行为也起着关键作用，受热力学原理的支配。根据操作条件，管道中的二氧化碳运输可以是单相的，也可以是两相的。为了准确描述两相流动，需要使用实际气体状态方程（EOS）。通常使用Mie-Grüneisen加强型气体EOS；参见[12]、[13]、[14]、[15]、[16]。然而，为了准确描述实际流体行为，需要更高级的EOS。这些高级EOS是根据亥姆霍兹自由能函数定义的：给定单相的温度、体积和组分，它们可以计算所有热力学性质。需要注意的是，这样的EOS是在单相水平上定义的；在两相系统中，不能直接从总体混合物量评估热力学性质，而必须先进行相分离计算。在流体动力学模拟中，每个时间步长都已知混合物的属性，即体积（网格尺寸）、速度和能量。从这些量出发，可以计算混合物的压力。这自然会导致相分离计算，这在热力学文献中通常称为闪蒸计算。

关于闪蒸计算有大量文献。大多数研究集中在PTN-flash（压力-温度指定）问题[17]、[18]、[19]、[20]。第二常见的情况是VTN-flash（体积-温度指定）[21]、[22]、[23]、[24]。PTN和VTN都是等温闪蒸问题。然而，关于非等温闪蒸的文献较少，例如PHN-flash（压力-焓指定）、UVN-flash（内能-体积指定）或SVN-flash（熵-体积指定）。对于流体动力学模拟，UVN-flash提供了一个自然的公式，因为每个计算单元中都有内能、比体积和摩尔组成的数据。在进行闪蒸计算之前，必须确定流体混合物是单相还是两相状态。这通过稳定性分析来确定。如果稳定性分析预测存在两相状态，则进行闪蒸计算。此外，稳定性分析的结果通常为后续的闪蒸计算提供一个很好的初始猜测值。

一些研究人员研究了管道中二氧化碳运输的数值模拟。Elshahomi等人[25]使用ANSYS软件进行了减压波速的数值模拟。近年来，挪威SINTEF能源研究机构在该领域处于领先地位。他们开发了一个专用的二氧化碳流动循环，并发表了关于二氧化碳运输的一系列研究[3]、[10]、[26]。他们的工作涵盖了含有杂质和不含有杂质的二氧化碳运输的数值建模，以及对热力学行为的详细研究。大多数研究要么关注流体动力学[3]、[26]、[27]、[28]，并直接使用现有的热力学库；要么专注于热力学行为[29]、[30]、[31]。

在本文中，我们提出了一个全面的框架，用于处理管道中流体混合物的两相多组分运输，同时考虑了流体动力学和热力学方面，采用UVN-flash框架。特别是，我们关注了流体动力学求解器输入带来的挑战对闪蒸问题鲁棒性的影响。为了解决这些挑战，我们通过一组替代变量重新制定了UVN-flash问题。这种转换提高了相平衡计算的鲁棒性。本文的关键贡献是为动态管道运输模型开发并测试了定制的闪蒸程序。所提出的方法通过储罐和管道减压案例进行了验证。

本文的结构如下：第2节介绍了HEM模型的控制方程和储罐减压的动力学。第3节讨论了热力学方面，包括热力学性质的计算和状态方程的选择。第4节讨论了稳定性分析，第5节讨论了针对流体动力学的UVN-flash的重新制定。第6节讨论了时间离散化方法。最后，第7节展示了结果并与文献数据进行了比较，第8节给出了结论。

2. 控制方程
本节描述了流体流动的控制方程。

2.1. 管道中的流体流动方程
纳维-斯托克斯（NS）方程描述了流体流动的动力学。原则上，可以求解三维方程；然而，对于可能长达数公里的管道，这会在计算上变得不可行。在这种情况下，使用横截面平均的一维方程更为高效[32]。对于完全分散的两相流体流动，控制方程可以用平均混合物量表示，从而得到均匀平衡模型（HEM）。对于水平无摩擦管道中的不可压缩流动，忽略热传递效应，HEM可以表示为：
(1)
?tU + ?xF(U) = S(U)，
其中U(x,t)是守恒变量向量，F(U)是通量向量，S(U)是源项。这些项的形式如下：
(2)
U = ρuρE，F(U) = ρuρu2 + p(ρE + p)u，S(U) = 0
这里，ρu,E,p分别表示混合物的质量密度、混合物速度、特定总能量和压力。HEM模型假设瞬时热力学平衡，并且两个相以相同的速度移动。总特定能量E可以表示为：
(3)
E = e + 1/2u2，
其中e是特定内能。混合物属性可以用以下方式表示：
(4a) ρ ? αρg + (1?α)ρl，
(4b) ρe ? αρge + (1?α)ρle，
其中α表示气相的体积分数，ρg和ρl分别表示气相和液相的质量密度，eg和el分别表示相应相的特定内能（质量基础）。

方程组(1)通过加入状态方程（EOS）来封闭，该方程将压力与其他热力学量相关联。EOS在封闭系统中的作用将在第3节详细讨论。此外，为了使问题适定，必须指定适当的初始和边界条件；这些方面将在第7节讨论。

2.2. 空间离散化
第2.1节描述的HEM模型使用有限体积法（FVM）求解。首先，我们在空间上离散系统(1)，得到以下半离散形式：
(5)
?U/dt = ?1/Δx(F?i+1/2 ? F?i?1/2)，
i=1…N，
其中Ui(t)∈R3≈U(xi,t)，表示单元i中的单元平均守恒变量，Δx是空间网格间距，F?i±1/2表示单元界面i±1/2处的数值通量，N是有限体积单元的总数。整个计算域的系统(5)可以用向量表示为：
(6)
dU(t)/dt = fPipe(U(t),V(t))，
其中U(t)∈R3?是整个域(node)的守恒量全局向量，定义为：
U(t) = [U1(t),…,Ui(t),…,UN(t)]?，
Ui=[ρi,(ρu)i,(ρE)i]?，
这里，U(t)由N个单独的单元向量Ui(t)组成，ρi,(ρu)i和(ρE)i分别表示第i个单元的质量密度、动量和总能量。向量V(t)包含与每个单元的热力学状态相关的非守恒（代数）变量。其精确定义将在第6节给出。

右侧的fPipe取决于守恒变量U(t)和相关非守恒变量V(t)，表示为：
(7)
?fPipe(U(t),V(t)) ? [fPipe, 1(t),…,fPipe, i(t),…,fPipe, N(t)]?，
其中每个局部fPipe, i表示为：
(8)
fPipe, i(Ui?r:i+r,Vi?r:i+r) = ?1/ΔxF?i+1/2 ? F?i?1/2，
其中r表示模板半径（半宽度），?Ui?r:i+r?(Ui?r,…,Ui+r)，Vi?r:i+r也是如此。数值通量F?i±1/2是根据宽度为2r+1的模板数据计算的，因此依赖于守恒变量Uj和相应的非守恒变量Vj（j=i?r,…,i+r）。我们采用一阶有限体积方法和HLLC Riemann求解器来计算数值通量。之所以选择一阶方案，主要是出于鲁棒性的考虑。在涉及强减压、相变和复杂热力学模型（如Peng-Robinson状态方程）的模拟中，高阶重构可能导致非物理的中间状态（例如负密度或压力），从而导致闪蒸计算失败。对于一阶方案，我们取r=1。为了计算单元间的数值通量F?i±1/2，我们使用HLLC方案；详见附录B。

2.3. 储罐减压
储罐模型可以被视为用单个计算单元离散化的管道。许多研究人员之前已经在UVN-flash背景下考虑过储罐的动态模拟[33]、[34]、[35]、[36]。此外，[37]、[38]、[39]也在管道流动背景下使用了该模型。描述储罐摩尔组成和内能演变的控制方程为：
(9a)
dNi/dt = N?iin ? N?iout，
对于i=1,…,n，
(9b)
dU/dt = H?in ? H?out + Q?，
其中n表示混合物中组分的总数。变量Ni表示储罐中组分i的摩尔数，U是总内能。N?iin和N?iout分别表示组分i在输入和输出流中的摩尔流量。同样，H?in和H?out表示与输入和输出流相关的焓流量，Q?是供给储罐的热量速率。方程（9a）—（9b）可以简洁地用向量形式表示为：(10)dUdt=fTank(U(t),V(t))，其中状态向量U(t)∈Rn+1定义为(11)U(t)=N1(t)?Nn(t)U(t)，向量V(t)包含与热力学状态相关的变量。这将在第6节中再次讨论。右侧的fTank(U,V)由(12)fTank(U,V)=N?1in?N?1out?N?nin?N?noutH?in?H?out+Q?给出。3. 热力学方面3.1. 定义为了清晰起见，我们在总内能（U）、体积（V）和摩尔数N的指定约束下定义以下热力学概念。参考相：参考相（?）表示一个假定的单相状态，其特征是总内能U?、体积V?和总摩尔数N?=(N1?,…,Nn?)。试验相：试验相是一个初始的（非常小的量）相，用于评估系统的热力学稳定性。通过添加这个相来扰动系统，如果总熵增加，则认为系统作为单相是不稳定的，相分离是有利的。稳定性分析：稳定性分析确定在指定的U,V和N下，给定流体是保持单相还是分离成两个相，即气液混合物。如果单相状态是不稳定的，该分析提供了其温度、摩尔浓度（单位体积的摩尔数）和摩尔内能密度的估计值。更多细节请参见第4节。闪蒸：闪蒸计算确定多组分混合物的平衡相分离。在指定的热力学约束下（例如，UVN），它计算共存相之间的广延量分布，即每个相的焓、内能、熵和体积，以及满足热力学、机械和化学平衡条件的常见强度量，例如温度、压力和组分化学势。3.2. 状态方程为了闭合HEM模型和储罐模型，我们需要一个本构关系来计算管道模型中的压力和储罐模型中的焓流量。这种本构模型基于流体混合物的热力学性质，通常被称为状态方程（EOS）。流体混合物的EOS是根据任意两个热力学变量（如压力和温度或体积和温度等）以及混合物中组分的摩尔数来指定的。因此，对于一个n组分混合物，需要n+2个变量来完全指定系统的热力学状态。根据选择的两个输入热力学变量，EOS可以分为两类：基于压力的和基于体积的。通常，真实流体的EOS是根据亥姆霍兹自由能A作为体积V、温度T和摩尔组成向量N的函数来定义的：(13)A=A(T,V,N)，其中?N?={N1,…,Nn}，Ni表示流体混合物中第i个组分的摩尔数。在这项工作中，我们使用Peng–Robinson EOS [40]。各种热力学性质可以直接从亥姆霍兹能量计算得出，如下所示。(14a)p=??A?VT,N(14b)S=??A?TV,N(14c)U=A+TS(14d)H=A+TS+pV(14e)Cv=?U?TV,N(14f)μi=?A?NiV,{Nj}j≠i。这里，所有导数和函数都是作为T、V和N的函数来计算的。EOS是针对单个均匀相定义的，必须应用于相变量，即每个相的温度、体积和摩尔组成。它不能直接应用于多相系统中的混合物平均量。然而，在实际情况下，只知道广延量，如总体积和总体摩尔组成，而不知道各相之间的总体积和摩尔数分布。这就需要一个相分离计算，通过这种方法确定共存相之间的广延量和常见强度量的分布（稍后会有更多讨论）。相分离计算在热力学文献中通常被称为闪蒸计算。3.3. 闪蒸的选择根据指定的两个热力学变量不同，会出现不同类型的闪蒸问题。文献中最常见的闪蒸问题是PTN-flash [17]，其中指定了压力、温度和总体摩尔组成。PTN-flash的目标是找到平衡相分离并计算每个相的广延性质，如内能、焓和体积等。另一个重要的闪蒸问题是VTN-flash，其中指定了混合物的总体积以及温度和总体组成，例如参见[41]。PTN和VTN闪蒸通常被称为等温闪蒸，因为温度是固定的。在封闭系统中，例如过程设计中的能量平衡计算，相关的规范通常包括总内能U、总体积V和总体摩尔组成N。这种非等温规范对应于所谓的UVN-flash问题。在我们的储罐模型和HEM模型中，我们通常可以在仿真的每个时间步访问总体积（对于储罐是固定的，对于管道是单元体积）、总内能和摩尔组成。因此，UVN-flash公式为确定平衡相分离提供了自然的方法，从而确定这些情况下的压力。4. 稳定性分析在本节中，我们简要回顾了稳定性分析程序。有关此主题的详细信息，我们建议感兴趣的读者参考现有文献[18]、[24]、[42]、[43]。4.1. UVN稳定性为了评估具有固定总内能U?、体积V?和摩尔数N?=(N1?,…,Nn?)的单相系统的热力学稳定性，我们考虑一个假设的相分离成两个相：一个主体相β和一个初始（无限小）相??。结果两相系统的总熵由下式给出：(15)???Stwo-phase=S(Uβ,Vβ,Nβ)+S(U?,V?,N?)，受守恒约束：(16)?U?=Uβ+U?,(17)?V?=Vβ+V?,(18)?Ni?=Niβ+Ni?,i=1…n。参考相的熵为：(19)Sref=S(U?,V?,N?)。可以通过对参考状态（U?,V?,N?）进行一阶泰勒展开来计算主体相的熵，将初始相的贡献视为无限小扰动：(20)?????Sβ?S(U??U?,V??V?,N1??N1?,…,Nn??Nn?),(21)???=S(U?,V?,N?)?U??S?UV,N?V??S?VU,N?∑i=1nNi??S?NiU,V,N。这导致(22)??????Stwo-phase=S(U?,V?,N?)?U??S?UV,N?V??S?VU,N?∑i=1nNi??S?NiU,V,N+S(U?,V?,N?)。然后熵差ΔS=Stwo-phase?Sref为：(23)??????ΔS=?U??S?UV,N(U?,V?,N?),?V??S?VU,N(U?,V?,N?)?∑i=1nNi??S?NiU,V,N(U?,V?,N?)+S(U?,V?,N?)。使用标准热力学恒等式(24)?S?UV,N=1T,?S?VU,N=PT,?S?NiU,V,N=?μiT，代入方程(23)，我们得到(25)??????ΔS=?U?T??P?T?V?+∑i=1nNi?μi?T?+S(U?,V?,N?)。此外，由于熵是一阶齐次的，欧拉定理给出初始相的熵为：(26)???????S(U?,V?,N?)=U?T?+P?T?V??∑i=1nμi?T?Ni?。将方程(26)代入方程(25)得到：(27)??????ΔS=1T??1T?U?+P?T??P?T?V??∑i=1nNi?μi?T??μi?T?。如果ΔS>0，则单相状态是不稳定的，混合物将分离成两个相。除以?V?，我们得到著名的切平面距离（TPD）函数(28)?????????D(c1?,…,cn?,u?)=ΔSV?,=u?1T??1T?+P?T??P?T??∑i=1nci?μi?T??μi?T?，其中????u??U?/V?是初始试验相中第i个组分的摩尔内能密度，???ci??Ni?/V?是摩尔浓度。注意，独立变量是?u?和??c1?,…,cn?。初始相的温度?T?、压力?P?和化学势?μi?可以通过以下方法获得。首先，通过反转内能关系来确定温度：????u?=U(T?,1.0,c1?,…,cn?)。一旦知道了?T?，就可以计算压力：????P?=P(T?,1.0,c1?,…,cn?)，化学势类似地得到：????μi?=μi(T?,1.0,c1?,…,cn?)，i=1,…,n。TPD函数D可用于具有指定UVN的系统的稳定性测试。如果对于所有可行的???{c1?,…,cn?,u?}，D≤0，则系统作为单相是稳定的。为了找到一个使D>0的状态???{c1?,…,cn?,u?}，需要寻找函数D的局部最大值。通过对D关于其独立变量?u?和??{c1?,…,cn?}求导并应用Gibbs—Duhem关系可以得到（详细信息参见[42]、[44]）(29)?1T??1T?=0,(30)?1T?(μi??μi?)=0，i=1,…,n。方程(29)意味着初始相的温度等于参考相的温度；方程(28)简化为(31)?????D=P?T??P?T??∑i=1nci?μi?T??μi?T?=1T?(P??P?)?∑i=1nci?(μi??μi?)。有趣的是，这与VTN稳定性分析中得到的TPD函数相同。[22]。因此，原则上可以解决VTN稳定性问题而不是UVN稳定性问题。可以使用n个方程(30)来求解初始相浓度?ci?,i=1,…,n，如果这组浓度的TPD结果大于零，则系统在热力学上是不稳定的，将会分离成多个相。注意，初始相的温度?T?与参考相的温度T?相同。此外，良好的初始猜测对于获得稳定性分析测试的收敛性至关重要，这将在下一节讨论。4.2. 稳定性分析的初始猜测我们采用三种不同的初始猜测方法进行稳定性测试：•如[42]、[43]所讨论的基于单纯形的方法。•如[18]、[44]所讨论的基于饱和压力的方法。•向假设的单相浓度添加高斯随机扰动。这三组初始条件提供了多样且分布良好的解空间覆盖，从而提高了稳定性分析的鲁棒性。通过探索多个起点，我们增加了检测初始相形成并获得试验相的浓度??c1?,…,cn?和摩尔内能密度?u?的可靠估计的可能性。实施稳定性分析时的实际考虑是跳过琐碎的解，即计算出的试验相摩尔浓度与输入摩尔浓度相同的解。4.3. 从稳定性结果中获得的初始猜测获得良好的初始猜测对于闪蒸计算的收敛性至关重要。稳定性分析提供了试验相的浓度和内能密度。为了开始闪蒸计算，我们还需要体积相分离。根据Kumar等人[43]的方法，试验相的体积最初设置为总体积的一半，然后迭代减半，直到找到一个导致两相熵大于相应单相状态的相分离。然后可以通过这个试验体积来缩放内能密度和摩尔浓度，以获得试验相内能和摩尔数的初始估计值，相应的主体相量通过从总数中减去来确定。5. 为流体动力学重新表述UVN-flash如果稳定性分析表明混合物将分离成两个相，则进行闪蒸以确定平衡相分离。在计算单元的给定时间步，可用的流体变量是总质量密度ρ、比内能e和总体混合物组成，表示为摩尔分数?z?{z1,…,zn}。根据这些量，可以计算总内能Utot和第i个组分的总摩尔数Nitot，如下所示：(32)Utot=ρeV,(33)Nitot=ziρV∑jzjMj,w,i=1,…,n，其中V是单元体积，Mj,w表示组分j的分子量。闪蒸计算的流程图如图1所示。值得注意的是，闪蒸计算可能会失败或收敛到非物理解，例如负体积分数，通常当初始猜测较差时。在这种情况下，可以按照[3]、[33]的建议，在内部循环中使用更稳健但计算成本更高的PT-flash。然而，在本文考虑的测试案例中，我们没有遇到这样的困难。这可能是由于从稳定性分析结果获得了高质量的初始猜测，并且在为流体动力学应用量身定制的UVN-flash公式中使用了适当缩放的变量。定制的闪蒸的详细信息将在以下小节中介绍和讨论。下载：下载高分辨率图像（576KB）下载：下载全尺寸图像图1. 闪蒸计算流程图。5.1. TVN空间中的UVN-flash对于总内能、体积和质量（即摩尔数）固定的封闭系统，平衡时的熵趋于最大。假设系统存在于两个相中，例如气体和液体。系统的总熵Stot可以表示为(34)Stot=S(Ug,Vg,Ng)+S(Ul,Vl,Nl)。液体相的量可以用总量和气体相的量表示，即Ul=Utot?Ug,Vl=V?Vg和N1l=N1tot?N1g,…,Nnl=Nntot?Nng。方程(34)可以重写为(35)Stot=S(Ug,Vg,Ng)+S(Utot?Ug,V?Vg,Nξ)，其中(36)?Nξ?{N1tot?N1g,…,Nntot?Nng}。可以使用气体相变量作为优化变量来最大化熵函数(35)，以确定平衡相分离。这种方法已被Castier [45]和Smejkal等人[42]采用。然而，EOS允许我们将熵和其他热力学量作为T,V和N的函数来计算。因此，在每次最大化子迭代中，必须通过解决(37)U=U(T,V,N)来获得温度，对于给定的内能U、体积V和组成N。在管道或组成储罐模拟中，需要执行成千上万甚至数百万次这样的计算，这些额外的牛顿迭代用于确定温度的反演步骤可能会成为瓶颈。这激发了直接用自然状态方程变量T、V和N重新表述熵最大化问题。利用热力学关系U=A+TS，给定内能Utot的熵可以表示为（38）STVN=Utot?ATVNT，其中上标TVN表示函数参数是T、V、N1、…、Nn，且?ATVN?∑k∈{g,l}A(T,Vk,Nk)。当在TVN空间中进行UVN闪蒸时，总体积V和摩尔数Nitot的刻度可能不合适：V的范围可以从毫米到千米，Nitot可以根据组成从接近零变化到数千。为了提高数值稳定性，我们改用气相组分密度?ρg?{ρ1g、…、ρng}以及气相体积分数αg∈[0,1]作为主要变量。这些变量在广泛的组成和相分离范围内自然保持良好的刻度，从而提高了闪蒸计算的稳健性。气相体积与总体积的关系为（39）Vg=αgV，每个相中的摩尔数根据（40）Nig=ρigVgMi,w,Nil=Ni?Nig得到，其中i=1、…、n，其逆关系为（41）ρig=Mi,wNigVg。有了这些变量，熵函数可以表示为（42）STαρ=Utot?ATαρT，其中上标Tαρ表示函数参数是T、αg、ρ1g、…、ρng，且?ATαρ?Ag+Al。相亥姆霍兹能量由（43）定义为??Ag?A(T,αg,ρ1g、…、ρng)，Al?A(T,αl,ρ1l、…、ρnl)，其中ρil=(Ni?Nig)Mi,w/(V?Vg)且αl=1?αg。液相密度ρil可以用气相变量表示为（44）ρil(ρig,αg)=Mi,wNi(1?αg)V?αg1?αgρig，i=1、…、n。因此，Al可以完全表示为(T,αg,ρg)的函数：Al=A(T,αl(αg),ρ1l(ρ1g,αg)、…、ρnl(ρng,αg))。从热力学我们知道，无论使用什么变量来表示熵函数，熵都应该在给定的UVN约束下最大化。因此，方程（42）中定义的熵的临界点必须与给定UVN约束下的热力学平衡点一致，这将在下一小节中验证。熵函数的临界点由（45）获得：?STαρ?T=0，?STαρ?αg=0，?STαρ?ρig=0（i=1、…、n），这会产生n+2个非线性方程。

5.2. 热力学一致性检查
在本节中，我们将证明由方程（45）定义的熵函数的临界点与相平衡的热力学条件一致。具体来说，只需证明这些临界点对应于压力相等、各相中各组分的化学势相等、温度相等以及满足给定的内能约束。温度相等的条件是显而易見的，因为在公式制定之初就假设了一个共同的温度。我们回忆一下（46）ρig=NigMi,wVg=NigMi,wαgV??ρig?Nig=Mi,wVg，?ρig?αg=?1αgNigMi,wαgV=?ρigαg。此外，逆关系由（47）给出：?Nig?ρig=VgMi,w，?αg?ρig=?αgρig，?Vg?ρig=?Vgρig。（1）关于质量密度的稳定性条件为：（48）0=?STαρ(T,αg,ρ1g、…、ρng)?ρigT,αg,{ρjg}j≠i=?1T?ATαρ(T,αg,ρ1g、…、ρnl)?ρigT,αg,{ρjg}j≠i。因此，我们关注总亥姆霍兹自由能对ρig的导数。使用ATαρ=Ag+Al，我们得到（49）?ATαρ?ρigT,αg,{ρjg}j≠i=?Ag?ρigT,αg,{ρjg}j≠i+?Al?ρigT,αg,{ρjg}j≠i。由于ρig是Nig和Vg的函数（见方程（46））；使用链式法则得到（50）?Ag?ρig=?Vg?ρig?Ag?VgT,{Nkg}+∑k?Nkg?ρig?Ag?NkgT,Vg,{Nmg}m≠k，（51）?Al?ρig=?Vl?ρig?Al?VlT,{Nkl}+∑k?Nkl?ρig?Al?NklT,Vl,{Nml}m≠k。在方程（50）、（51）中，只有当k=i时求和项才有效。使用化学势的定义μig=?Ag?NigT,Vg,{Njg}j≠i并重新排列后，我们得到（52）?ATαρ?ρig=μig?Nig?ρig+μil?Nil?ρig+?Vg?ρig?Ag?VgT,{Nkg}︸?pg+?Vl?ρig?Al?VlT,{Nkl}︸?pl。由于Vl=V?Vg，因此?Vl?ρig=??Vg?ρig。根据方程（46），我们有?Vg?ρig=?Vg/ρig，?Nig?ρig=Vg/Mi,w。此外，使用质量平衡关系Nil=Ni?Nig，我们得到?Nil?ρig=?Vg/Mi,w。将这些表达式代入方程（52）并使用方程（49）得到（53）?ATαρ?ρig=μigVgMi,w?μilVgMi,w+pg?plVgρig=μig?μilVgMi,w+pg?plVgρig。然后稳定性条件（48）产生：（54）0=(μig?μil)ρigMi,w+(pg?pl)，?i∈{1、…、n}。注意，方程（54）是一个n个方程的系统，其展开形式为：（55a）(μ1g?μ1l)ρ1gM1,w+(pg?pl)=0，（55b）?（55c）(μng?μnl)ρngMn,w+(pg?pl)=0。（2）关于气相分数αg的稳定性条件为：（56）?ATαρ(T,αg,ρ1g、…、ρng)?αgT,ρ1g、…、ρng=0。代入ATαρ=Ag+Al并对αg求导，我们得到（57）?ATαρ?αg=?Ag?αg+?Al?αg。由于αg=Vg/V，对αg的偏导数为（58）?Ag?αg=?Ag?Vg?Vg?αg=1V?Ag?Vg。对于液相，使用Vl=V?Vg，我们得到（59）?Al?αg=?Al?Vl?Vl?αg=?1V?Al?Vl。将两者结合并回忆?A?V=?p，我们得到（60）?ATαρ?αg=pl?pgV。在稳定性点?ATαρ?αg=0，这意味着（61），对应于机械平衡的条件。将pg=pl代入方程（55）得到（62），这是化学平衡的条件。（3）关于温度T的稳定性条件为：（63）?STαρ(T,αg,ρ1g、…、ρng)?Tαg,ρ1g、…、ρng=0。回忆STαρ=(Utot?ATαρ)/T。对T求导得到：（64）?STαρ?T=?UtotT2?1T?ATαρ?T+ATαρT2。代入ATαρ=Ag+Al并对T求导得到：（65）?ATαρ?T=?Ag?T+?Al?T。使用热力学恒等式?A?T=?S，我们得到（66）?ATαρ?T=?Sg?Sl。将此代入方程（64），我们得到（67）?STαρ?T=?UtotT2+Sg+SlT+Ag+AlT2，（68）=Ag+TSg+Al+TSl?UtotT2。对每个相应用关系U=A+TS得到（69）?STαρ?T=Ug+Ul?UtotT2。因此，稳定性条件?STαρ?T=0产生（70），这是对内能的期望约束。方程（61）、（62）和（70）一起构成了在UVN规格下的热力学平衡条件。

6. 时间离散化
我们考虑由管道流动方程（6）和储罐模型（10）以及闪蒸问题（45）产生的微分-代数系统（DAE）。将状态变量分为微分（守恒）变量U和代数（非守恒）变量V，得到的DAE系统可以表示为（71）dUdt=f(U,V），（72）0=g(U,V)。下面讨论了管道和储罐模型的U和V的精确形式。

管道流动
对于空间离散化的管道流动，我们定义：U=[U1、…、UN]T，Ui=[ρi、(ρu)i、(ρE)i]T，V=[V1、…、VN]T，Vi=[ρk,ig,αig,Ti]T，其中ρkg、αg和T分别是气相中组分k的偏质量密度、气相体积分数和温度。下标i表示第i个单元格。对守恒定律进行一阶空间离散化得到半离散形式：?fi(U,V)??1ΔxF?i+1/2(Ui+1,Ui,Vi+1,Vi)?F?i?1/2(Ui,Ui?1,Vi,Vi?1)，其中F?i±1/2表示数值通量函数，fi代表系统（1）的离散化右侧。

储罐模型
对于储罐，微分变量定义为：U=[N1、…、Nn,U]T，其中Ni表示混合物中组分i的摩尔总数，U是总内能。代数变量V定义为：V=[N1g、…、Nng,Vg,T]T，其中Nig是气相中组分i的摩尔数，Vg是气相所占的体积，T是温度。储罐的函数由方程（12）定义。

时间积分方案
为了对DAE系统（71）进行时间积分，我们采用半显式前向Euler方案[46]。该方法包括对微分变量进行显式更新，然后解决非线性约束以强制代数约束。首先，计算微分更新：Un+1=Un+Δtf(Un,Vn)，然后解决非线性约束：0=g(Un+1,Vn+1)，以确定Vn+1。解决这个约束相当于按照第5.7节中方程（45）制定的闪蒸计算。

数值实验
在本节中，我们讨论结果。首先讨论储罐模型的结果，然后是管道减压的结果。我们对所有结果使用PR-EOS；详见附录D。PR的二元相互作用参数（BIP）矩阵的对角元素设置为零，非对角项在表1、表3中指定。

表1. Castier问题1和2的初始条件、流数据
数量 单位
问题1
问题2
组分（顺序）–{CH4, H2S}{CO2, C12H26, C13H28, C14H30, C15H32}
储罐初始条件
压力 P0MPa 0.10 106
0.00 106
温度 T0K 300.0 373.15
总摩尔数 N0mol {500,500}{1×10?8,0.1,0.6,0.2,0.1}
体积 Vm3 24.57 084.7
1×10?4
初始相态 –单相
两相
入口流
压力 PinMPa 520
温度 TinK 300.0 300.03 10.0
摩尔流量 finmol {4.0,6.0}{1×10?2,0.0,0.0,0.0,0.0}
热输入 Q?J/s 0.0 0.0 10.0
-100.0
相态（入口）–两相
单相
出口流
无
None
BIP ({δij}i≠j) 0.0 83.0 7.1

7.1.1. 问题1：轻组分
在这个问题中，我们考虑甲烷（CH4）和硫化氢（H2S）的混合物。储罐中的流体在t=0s时处于单相状态。有一个入口流流入储罐，该流处于两相状态。由于流的条件以温度和压力给出，我们需要在模拟开始时执行一次PTN闪蒸[17]以确定相分离，从而得到各个相的焓。总焓是通过将这些单独的相贡献相加得到的。结果表明与Castier [33]报告的结果非常吻合。最初温度急剧下降，在大约680秒时达到最低温度（约260K），之后温度开始上升，在大约1792秒时出现第二个相。此时温度斜率发生变化。另一方面，随着更多质量加入储罐，压力持续线性增加。

7.1.2. 问题2：CO2加载
在这个例子中，我们有中重碳氢化合物的混合物。最初储罐中的CO2量很少。储罐正在接收CO2。对于给定的压力、温度和进料组成，我们可以求解PR-EOS以获得体积，从而得到焓流量。最初，储罐中的流体处于两相状态。在整个过程中，以恒定速率100 J/s移除热量。图3显示了当前工作的结果，并将其与Castier [33]的结果进行了比较。同样，可以观察到非常好的匹配。直到大约341秒，储罐中还有两个相；过了这一点，只剩下一个相。然后温度相对恒定，但压力急剧增加，表明流体已经转变为液态，表现得几乎像不可压缩流体。

7.2. 管道减压
在验证了储罐模型的闪蒸求解器用于瞬态模拟之后，我们现在将其应用于管道传输的多组分流体。我们考虑了六种不同的混合物：4种两组分混合物；一种五组分混合物；以及一种单组分流体。此类数值实验已经在以前的研究中考虑过；参见参考文献[27]、[38]、[39]（针对单组分流体）和Munkejord等人[3]（针对多组分混合物）。在Munkejord等人[3]中，在管道的开放端（右侧）规定了出流边界条件。相比之下，在我们的模拟中，我们避免直接施加这一边界条件。相反，我们将计算域扩展到原始管道长度的两倍，并将问题表述为一个在中点具有初始不连续性的冲击管。采用这种方法的原因是开放端的流动条件会被阻塞，因此外部信息无法向计算域上游传播。模拟在最快行波到达管道末端之前终止。

在所有模拟中，使用了800个单元格的网格（即每侧400个单元格），而Munkejord等人[3]使用了2400个单元格。此外，我们使用了CFL值0.9，而Munkejord等人[3]使用了CFL值0.85。由于存在初始不连续性，因此使用1×10?12的初始时间步长开始模拟。后续时间步长根据最快行波确定：（73）Δt=Δxmaxi(|ui±ai|)，其中下标i对应于第i个单元格，ai表示声速。为了计算HEM声速（cHEM），需要沿饱和曲线的导数。然而，对于多组分系统，两组分包络线取代了单组分饱和曲线，使得cHEM的评估变得复杂得多。因此，我们使用Wood声速（cWood）作为实际替代方案（见附录C）。虽然HEM假设完全的热力学平衡，但Wood声速对应于“冻结”状态，即各相处于机械平衡但不是热平衡和化学平衡。然而，正如Picard等人[47]所指出的，在压力波传播期间，各相之间的质传或热传递时间不足，表明在这种情况下cWood是一个物理上合理的近似。此外，根据Fl?tten等人[48]讨论的次特征条件，平衡声速总是被冻结声速(cWood)所限制，即cHEM≤cWood。因此，使用cWood确定时间步长可以保守地估计最大信号速度，确保数值方案的稳定性，而不会违反CFL条件。另外，如2.2节所讨论的，我们使用一阶精确的HLLC方案进行空间离散化总摩尔数N是通过关系式Z=pVNRT计算得出的，其中V是计算单元的总体积，R是通用气体常数。然后通过将N乘以表3中指定的组成向量来获得摩尔向量。表2. Riemann问题初始条件摘要。案例液体L [m]压强 [MPa]温度 [K]空单元空单元空单元左右左右1{CO2, CH4 }200.028.5682.0313.65293.152{CO2, H2, N2, O2, CH4 }288.015.04.5283.15283.153{CO2 }200.010.03.0300.0300.04{CO2, N2 }283.811.992.0292.65291.555{CO2, N2 }283.812.082.0292.85292.656{CO2, N2 }283.812.02.0290.45292.35表3. 不同问题的组成。案例αLαRBIP ({δij}i≠j)CO2CH4H2N2O210.01.00.150.7260.264–––20.01.00.00.91030.01950.01150.040.018730.01.0–1.0––––40.01.0?0.0410.9––0.1–50.01.0?0.0410.8––0.2–60.01.0?0.0410.7––0.3–7.3. 验证在这里，我们重点关注使用Munkejord等人[3]报告的结果对富含CO2的混合物进行验证，具体规格见表2和表3。需要强调的是，这里考虑的设置与Munkejord等人的设置之间存在差异，这些差异可能导致结果上的不同。主要区别在于：(i)我们考虑的是一个冲击管设置，而Munkejord等人施加了一个边界条件；(ii)Peng–Robinson (PR)系数可能有所不同，因为Munkejord等人没有明确报告这些系数。尽管在建模选择上存在差异，但我们的结果与文献中的数据非常吻合，如图4所示，其中显示了二元混合物和五组分混合物在压力-温度(PT)平面上的HEM轨迹。值得注意的是，由于我们没有考虑热传递和摩擦，管道左侧部分的模拟路径对应于一个等熵过程。图5展示了t=0.1s时CO2和CH4二元混合物沿管道长度的压力和温度结果。在x≈120m处，温度剖面的斜率发生了轻微变化，这一点用橙色矩形标出。这种现象在多组分混合物中是固有的，在单组分系统中则不存在，如图10所示，其中蒸发波突然转变为接触波。这种差异可以通过相图的拓扑结构来解释：对于纯物质，饱和点在PT空间中形成一条单曲线；见图6(a)，每个压力-温度点可以对应多个相分数值。然而，对于多组分混合物，两相区域由泡线和露点线围成一个 envelop；见图6(b)。在这个 envelop 内，等质量线（相分数保持不变）从泡线（蒸汽分数0）延伸到露点线（蒸汽分数1）。因此，encapsule 内的每个压力-温度条件都对应一个唯一的蒸汽分数值。图7(a)展示了速度剖面。请注意，流体速度在接触波两侧是变化的，这与单相情况不同。这一点也在[39]、[49]中观察到。图7(b)显示了管道中（伍德的）声速的空间变化。声速是近似Riemann求解器（本例中为HLLC）中计算波速的一个重要组成部分；见附录B。与单相区域相比，可以在两相区域观察到声速的急剧下降。这种现象归因于当液体和气体相同时流体的可压缩性增加，以及气泡减弱了压力扰动。图8展示了五组分混合物的结果。整体波形与二元情况非常相似。由于初始条件、组成和EOS系数的不同，可以看到波幅和斜率（用橙色矩形标出）的轻微差异。图9(a)提供了斜坡区域的放大视图。值得注意的是，如图8(b)所示，该区域的压力没有变化。此外，如图9(b)所示，该区域对应于从两相状态到单相状态的过渡。最后，图10展示了单组分CO2情况下的空间压力-温度（图10(a)）和（图10(b)）剖面。总之，所提出的统一框架已在多个复杂度的基准问题上进行了测试，其结果与文献中的数据一致，从而证明了其稳健性。需要强调的是，稳定性分析和闪蒸计算在计算上要求很高，因此它们的相应程序是未来优化的自然候选对象。下载：下载高分辨率图像（801KB）下载：下载全尺寸图像图4. PT空间中带有800个单元的模拟路径。下载：下载高分辨率图像（256KB）下载：下载全尺寸图像图5. t=0.1s时带有800个单元的CO2和CH4混合物的冲击管结果。下载：下载高分辨率图像（327KB）下载：下载全尺寸图像图6. PT空间中带有800个单元的模拟路径。下载：下载高分辨率图像（265KB）下载：下载全尺寸图像图7. t=0.1s时带有800个单元的CO2和CH4混合物的冲击管结果。7.4. 波结构对于所有三个测试案例，解决方案都表现出与两相Riemann问题的典型波形一致的结构：1. 向左移动的稀疏波：这种平滑的波对应于管道中流体的膨胀，导致波上传压力和温度的下降。在Riemann问题的背景下，这是一种真正的非线性波[50]。2. 蒸発前沿：这一区域发生液-气相变。值得注意的是，这种波在单相情况下的Riemann问题中是不存在的。此外，这一区域与向左移动的快速稀疏波一起形成了一个复合稀疏波。具体来说，对于(a)纯组分：蒸发前沿突然结束，并以接触不连续性过渡到单相区域。(b)多组分：蒸发前沿表现出一个斜率区域，对应于PT空间中的不同蒸汽相分数值。之后，它以接触不连续性过渡到单相区域。3. 向右移动的接触不连续性：这种波跟随冲击波，并标志着物质的分离。在这一波中不发生流体混合。与单相流体不同，它伴随着流体速度的变化。4. 向右移动的冲击波：这种波对应于流体压缩，伴随着熵的产生，并导致温度和压力的增加。下载：下载高分辨率图像（268KB）下载：下载全尺寸图像图8. t=0.22s时带有800个单元的五组分混合物的冲击管结果。下载：下载高分辨率图像（255KB）下载：下载全尺寸图像图9. t=0.22s时带有800个网格单元的五组分混合物的冲击管结果。下载：下载高分辨率图像（249KB）下载：下载全尺寸图像图10. t=0.2s时带有800个网格单元的单组分（CO2）的冲击管结果。8. 结论在本文中，我们开发了一个统一的框架，该框架集成了流体动力学和热力学，用于模拟两相条件下的多组分流体传输。作为测试案例，我们考虑了与CCS应用相关的富含CO2的混合物。我们采用了均匀平衡模型(HEM)来描述流动，并使用基于亥姆霍兹自由能的状态方程来模拟热力学性质。关键贡献是开发并测试了一个专为动态管道传输模型设计的UVN-flash框架。UVN-flash的重新表述依赖于引入了一组新的变量，即ρ1g,...,ρng,α,T。这种选择与流体动力学求解器的输入无缝对接，与标准的UVN/TVN变量相比，提供了更好的缩放变量。此外，我们证明了当用这些变量表达熵函数时，其临界点对应于经典的热力学平衡条件。此外，我们还讨论了UVN-flash的稳定性分析公式，这使我们能够可靠地检测单相和两相状态。这是进行稳健的多相瞬态管道模拟的关键组成部分。该方法在一系列减压场景上进行了测试，包括两个罐问题（一个二元混合物和一个五组分混合物）以及纯CO2、四种二元混合物和一个五组分混合物的管道案例。数值实验表明，所提出的框架能够一致且准确地描述耦合的热力学和流体动力学问题，并能够捕捉到对安全管道设计相关的极端条件。还提供了算法的流程图以及所得DAE系统的时间和空间离散化细节。由于热力学计算在计算上非常耗时，探索加速这些计算的策略是一个有趣的研究方向。另一个有前景的研究方向可能是研究SVN-flash在管道减压中的应用，因为它可以提供一种高效且准确的方法来模拟等熵流的热力学。cRediT作者贡献声明Pardeep Kumar：写作——原始草稿，可视化，验证，软件，方法论，概念化。Patricio I. Rosen Esquivel：写作——审阅与编辑，监督，项目管理，资金获取，概念化。在撰写过程中使用生成式AI和AI辅助技术的声明在准备这项工作时，作者使用了GitHub Copilot来提出措辞和数学排版。使用该工具/服务后，作者根据需要审阅和编辑了内容。作者对出版物的内容负全责。