本文概述了D.B.Spalding教授迄今在科学和工程方面的一些贡献。从Spalding早期关于燃烧的工作和他在传质理论方面的独特工作出发,简要介绍了Spalding未发表的“统一理论”。在此之后,帝国理工学院小组在算法方面的发展导致了现代计算流体动力学的诞生,包括著名的SIMPLE算法。还介绍了燃烧、多相流和湍流模型的发展。最后,提到了计算流体动力学的一些学术和工业应用以及随后几年考虑的传热应用。
在这种性质的论文中,不需要使用通常的命名法部分。Spalding对科学和诗歌的热情都是广为人知的。
也许正是这些看似不相干的追求之间的相互作用,赋予了他发明奇妙的描述性缩写词的才能,这些缩写词描述了他的许多科学贡献。这个清单很长;这里只提供了几个。其他的则散布在整个论文中。
a有时也被称为相间滑移分析仪
1960年,Spalding在本刊第1卷发表了一篇开创性的文章[10],在这篇文章中,他假设了质量传递的“欧姆定律”,即
=gB (1)
右边的第一项代表层流的贡献,其余的是湍流的贡献。一旦用这种方式表示有效粘度,就有可能扩展这一思想来得到有效导热系数的表达式
对于中等普朗特数,这个表达式有时简化为。该函数至今仍被广泛用于紊流计算中提供壁面边界条件。
1964年,Spalding发表了一份航空研究委员会报告[21],标题为“湍流边界层和壁面射流中摩擦、传热和传质的统一理论”。1964年12月,这份报告被重新打印并做了注释[22]:有87页打印的傻瓜页和29张图。这是Spalding对紊流剪切流的预测作出非凡贡献的真正开端。似乎标题的范围不够宏大,论文的图1(此处转载)揭示了他的视野的范围,尽管正如他在论文的第1页所明确指出的那样,并非所有所代表的特征通常同时出现。值得注意的是,这篇论文从未在期刊上发表过,即使是以缩写形式发表。一种可能的解释是,无论是在他自己的小组内部还是在其他地方,发展都非常迅速,以至于Spalding认为统一理论很快就会被取代,不值得进一步研究。
① Spalding自己最近评论说,他觉得这篇论文“对一本书来说太短了,对一篇期刊论文来说也太长了”。
其中为无量纲速度,为自由流速度,壁面函数u+ (y+)的适当形式包含传质效应。这个方程的关键特征,与Coles提出的著名关系非常接近,是参数的出现,据说,“将在后面发挥重要作用”。引入了守恒性质的一个类似方程,它可以表示混合物中化学惰性成分的组成、停滞焓等。在平面自由湍流混合层信息的指导下,引入了初步的夹带规律。Head曾试图将夹带率与形状因子(定义为联系起来,yG是速度边界层厚度,而Spalding则展示了他通常的物理洞察力,意识到夹带率更可能与(1-)成正比,即速度剖面的射流强度或尾流分量。事实上,在最初的统一理论[21]中,他发展了一对更复杂的方程
其中为无因次夹带速率。这两个方程后来都被Spalding简化了[22],因为他们被认为过于复杂,并且发现高估了夹带速率,并在不久之后被他进一步修改了[26]。统一理论被证明与湍流边界层理论的许多其他方面是一致的,事实上,它表明了如何发展一个统一的方法来涵盖诸如传热、传质、粗糙度甚至性质变化等影响。
发展恒定性质边界层统一理论的纯流体力学方面的责任交给了两个研究生,Escudier和Nicoll。湍流边界层(<1)的实验数据由前者进行分析,而壁面射流(>1)的实验数据较少因此需要大量的实验工作。夹带率是用公式推导出来的
本文进行了大量的计算和讨论,最后得出结论:“边界层的预测与实验基本一致,除了接近分离条件的地方”和“壁面射流的预测是好的,无论是对零压力梯度还是对逆压力梯度”。另一方面,Escudier和Nicoll对逆压力梯度边界层性能差给出了通常的解释,即将实验数据不满足二维条件归咎于实验数据。这里有些道理;参见Fellouah and Pollard[28]。前面已经提到,在《统一理论》发表的第二年,Spalding修改了无量纲带射与尾迹参数的关系方程[26]。实际上,正是在后一篇论文中,Spalding首次考虑了一种基于求解湍流边界层动量积分方程和动能亏缺方程的边界层预测方法。他表明,在Rotta[29]和Truckenbrodt[30]的成熟理论中,存在带射函数和耗散积分之间的形式关系。他继续推荐改进的,但暂时的表达式,对<1有效,包括近分离的边界层,以及对>1有效。Spalding说:“作者目前的观点是,最终使用‘动能’法比使用‘质量守恒’法更可取。”他的理由是“更直接地将与密度变化、曲率半径和其他效应联系起来”。很明显,在Spalding的眼中,到1965年,夹带理论实际上已经成为历史:他说“不应该忘记,基于夹带理论的作用,尽管可能很短暂,但很有价值”,但没有迹象表明,消散积分方法也注定是短命的。
基于动量方程和动能方程解的边界层流动流体力学方面的预测,后来通过使用温度剖面的两个参数:热边界层厚度和尾迹参数,扩展到温度方程。这些是由能量方程和“T-平方”方程得到的,这是动能方程的对应项。该研究发表在1966年芝加哥国际传热会议上[32]。这时,Spalding开始有效地使用加权函数的思想。动能方程和“T-平方”方程不需要物理上可理解的基础;它们可以被认为是动量和能量方程的积分,乘以一个加权函数(可以是因变量本身)。为了寻求更大的灵活性和更广泛的适用性,考虑了许多加权函数族。例如,利用u(或T)的幂得到了动能和相关方程。y的幂(作为加权函数)导致了各种“动量矩”方程。另一个选择是使用子域加权函数。有了这么多的积分方程,轮廓现在可以有大量的参数;它们可以是多项式,或者更好的是,分段线性。在这段时间的探索中,Spalding还提出了边界层中横流坐标的各种自变量。距离y或无量纲距离是明显的候选;但他提出了流函数(及其无量纲对应物),它简化了横流速度或流量的表示。
Spalding当时很有信心,他的统一理论和分段剖面图可以表示许多工程兴趣流。尽管该理论取得了许多显著的进展,但具有强大的有利或不利压力梯度的流动仍然超出了UT的范围。这些流动包括重要的工业流动类别,如壁面射流和分离流动,包括管道台阶后面或腔体上方的流动。统一理论的主要问题是需要使用一组代数方程来计算轮廓参数。这些方程可以产生奇点,例如,在壁面射流到边界层的过渡中。
剪切在分离流动中起着关键作用,包括那些边界层被不利压力梯度破坏的流动,或者导致分离的几何形状。由于UT的核心依赖于Navier-Stokes方程的抛物线版本,因此这种流动不在UT的范围之内。Spalding的想法是,可以找到能够进化并代表这些流动的概况,否则这些流动就没有“自相似”行为。因此,他指派了两名新学生来解决将这种流动纳入统一理论范围的问题:沃尔夫施因已经在1964年8月加入了这个小组,Spalding让他解决平板上撞击射流的问题;在一个通道中向后台阶后面的流动和在一个驱动的方形腔内的流动是Runchal关注的焦点。这些扩展将牢固地建立统一理论,不仅适用于边界层等抛物流,而且适用于具有强压力梯度、再循环和碰撞的椭圆流。到1965年底,从Runchal和Wolfshtein的初步工作来看,UT在分离或高压梯度流动方面明显存在两个严重缺陷。UT的基本结构是抛物型Navier-Stokes方程,并包含“自相似性”存在的固有假设。Spalding意识到UT和剖面法有一个致命的缺陷。轴向扩散破坏了自相似性,并在分离和反向流动中起着关键作用,因此没有简单的方法来表示轴向扩散的作用。UT不能扩大到这种流动。在尝试将UT扩展到强压力梯度流动的过程中,很明显,当时唯一可用的解决Navier-Stokes方程的完整(椭圆)形式的一般方法是有限差分法。在某种程度上,有限差分法(FDM)与统一理论的中心思想相似:两者都使用分段多项式轮廓。但也有一个本质的区别。UT依赖于求解完成配置文件的匹配常数,而FDM使用简单的通用分段配置文件直接求解控制变量的值。本质上,在计算域中,除了感兴趣点(或节点)的近邻之外,没有明确的轮廓。这可能被认为是一个局部发展的概要,而不是一个通用的概要。
Navier-Stokes方程的有限差分方法已经存在很长时间了:Thom[33]在电子计算机出现之前就已经很好地使用了它们,Burggraf[34]最近将它们应用于方形腔体中的分离流动。很明显,使用FDM方法比修改UT方法更简单、更通用。1966年1月,Runchal和Wolfshtein结合了他们各自的研究目标,开始编写一个通用的FDM计算机程序。由于他们的直接兴趣是解决二维问题,纳维-斯托克斯方程的流函数涡度(φ-ω)形式是他们关注的焦点。φ-ω方法消除了方程中的压力,也产生了一组两个而不是三个耦合的偏微分方程。在早期,很明显存在标准的局限性,然后首选使用FDM的中心差异选项。对于高雷诺数流,解不能收敛。稳定性分析表明,这是因为在高雷诺数下,矩阵系数不再是正定的。Spalding接着做了一个类比:猪圈里吹来的风总是臭烘烘的——或者(北半球的)北风通常带来寒冷。这些讨论导致了“逆风”概念的发展和使用。后来发现,Courant等[35]等人之前也使用了单侧差异。Spalding做出的重要贡献是采用了有限体积方法(FVM)的概念作为数学发展的核心,而不是使用传统的fdm。③Spalding将计算网格比作一系列的“罐”(容纳节点)和“管”(沿着网格输送流体)。这些改进迅速带来了成功。到1966年6月,Runchal和Wolfshtein报告了壁面射流和后台阶后分离流动以及驱动盖方形腔内低和高雷诺数流动的结果。这在帝国理工学院的一系列论文中有报道[36-39]。到1966年末,这种方法已经建立起来,并开始被世界各地的其他研究人员广泛使用。它随后由 Gosman等人以书的形式出版。[40]
③ 直到一段时间后,“有限体积”一词才被采用。
φ-ω方法的一个缺点是,它基本上局限于二维流动,不能很容易地以任何明显的方式扩展到三维。使用矢量涡度和矢量速度势(与流函数对应)的可能扩展在当时被认为过于复杂;人们担心二维领域令人愉快的简单性会消失,因此进一步的进展(如下所述)是基于所谓的原始变量,即速度分量和压力。然而,涡度-矢量势公式实际上后来被Mallinson和de Vahl Davis[41]以及其他许多人成功地使用。
随着FVM的发展和成功,以及无法将UT扩展到具有强压力梯度的流动,UT的进展(基于积分方程)在1966年夏天突然停止,当时Spalding与帕坦卡尔一起前往芝加哥参加1966年国际传热会议,发表了一篇关于UT工作的论文[32],并在暑假期间进行了一个咨询项目。在此期间,随着椭圆流的FDM已经证明的成功,他们决定基本上放弃UT的“分段剖面”方法,转而使用FDM来处理二维边界层的抛物形式的Navier-Stokes方程。他们引入了一项创新,用一个基于变量的无维流函数(本质上是von-Mises变换)代替与墙的法向距离,作为横流坐标、速度(和温度)的分段线性剖面和子域加权函数。在自由边界处使用合适的夹带公式来确定计算域向周围自由流的增长。该方法具有二维边界层、自由剪切流、壁面射流和管道流的一般程序所需的速度、灵活性、准确性和简洁性。
因此,到1966年底,椭圆型和抛物线型FDMs都牢固地建立起来,整个集团的重点转移到该技术的应用和改进上。随着论文[42,43]的发表,以及后来Patankar和Spalding[44,45]和Gosman等人[40]的著作的出版,椭圆和抛物线方法都获得了相当多的追随者,并被帝国理工学院、斯坦福大学和其他地方的研究人员广泛使用。许多博士论文将该方法应用于各种抛物线和椭圆流。这也是对湍流数学模型的研究如火如荼的时候。CFD代码为湍流和燃烧模型的测试提供了方便的工具。辐射的二通量和六通量模型[46]以及不同粒径范围的粒子浓度公式也在这些方法的框架内实现。
到1968年,越来越多的证据表明,在许多应用中,具有逆风差异的FDM和FMV会导致过度的假(数值)扩散。因此,努力转向寻找更好的方法来表示对流项。研究了许多高阶格式,包括沿流线的矢量差分概念,但没有一个证明是完全令人满意的。Spalding[47]随后提出了一种指数方法来代替上风向,但最终帝国理工学院的研究小组最终确定了一种“混合”差分方案[48],该方案基于局部Peclet数混合了中心和上风向差分方法。在相当长的一段时间内,这[40,43]成为求解高雷诺数流动的标准方法,但后来被更高阶上旋、QUICK[49]和矢量上旋方案所扩展,并通过使用通量而不是变量的剖面。
人们意识到,“罐管”方法可以被高斯定理所取代。这避免了泰勒级数的使用,并导致在任意形状的控制体积中导出离散方程组的优雅形式。
这种广义的FVM通过产生计算流体力学(CFD)的新领域而彻底改变了工程实践。通用和灵活的计算机代码可以用来解决任何复杂几何形状和任何雷诺数下的流体流动问题。在一个重要的方面,FVM也改变了计算流体动力学背后的基本思想。与FDM和有限元法(FEM)使用的数学方法不同,对FVM的理解是通过物理方法来增强的。感兴趣的焦点不是变量及其与其他节点的相互关系参数;现在的重点是控制体积和跨越它们边界的通量,以及其中包含的源或汇。
进一步的研究和Spalding的持续创新,导致了抛物边界层方法的一个新的改进版本,由Patankar和Spalding在该书的1970年版中发表[45]。其中包含一个名为GENMIX的新计算机程序。1977年,Spalding出版了一本书[50],本质上是一本GENMIX的用户手册。该方法现在包括一种改进的计算夹带的方法和一种确定受限流动中未知压力梯度的程序。新方法现在包含了先前为高雷诺数流动椭圆模型开发的混合格式。
Patankar-Spalding抛物边界层过程基于一种纯“行进”过程的非迭代格式。Spalding通过“汽车转向类比”解释了在受限流动中压力梯度的非迭代确定。在该方法的许多应用中,自由边界夹带的计算仍然是振荡和不稳定的来源。从形式上讲,这些方程是非线性的,迭代解将以一种直接的方式解决这个问题。(稍后在通用CFD代码PHOENICS中,Spalding[51]解决了这些问题,并在前进的步骤中使用了迭代)夹带速率(以及在受限流动中未知的压力梯度)应该被认为是方程中的特征值,由所得解的某些总体性质确定。根据我们现在所知道的,现在可以构建一个更健壮和有效的二维边界层程序;但在当今时代,几乎没有动力去开发这样的方法。
1967年,Spalding做出了一项看似微不足道的贡献,但对后来CFD工具的思考和发展产生了深远的影响。他引入了“一般输运方程”的概念来表达所有二阶对流-扩散输运方程。他建议把它们写成一个方程,用一个广义变量u来表示速度分量、焓或化学物质。这个广义输运方程有四个组成部分:(1)瞬态或累积,(2)对流,(3)扩散,(4)源项。因此,Navier - Stokes、能量和质量(物种)输运方程都可以用这个方程来表示:
这种表述的早期例子出现在Gosman等人1969年的书中[40]。这是他对数学表达“统一”的努力。这种方法导致了CFD方法的快速发展,因为人们可以集中精力解决单个方程,而不是查看每个变量的细节。
Spalding在1966年左右引入的另一项创新是用“罗盘”符号来表示FD或FV模板。这简化了通量和的可视化和表达,与他对FV的“坦克和管道”方法相吻合。围绕一个中心点(或CV)的六个点(或控制体)被指定为“北”、“南”、“东”、“西”、“高”和“低”。因此,线性代数方程的“标准”形式,即所谓的零残差形式为
到1969年,帝国理工学院团队的兴趣开始转向三维流体。流函数涡度变量在二维流动中具有明显的优势,但如前所述,在三维流动中使用它们是有问题的。在对这种流动进行了一些初步的(φ-ω)扩展之后,Spalding开始寻找直接基于Navier-Stokes方程原始形式的方法。他们尝试了许多备选办法,但没有一种完全令人满意。Caretto等人[52]的耦合SIVA方案带来了第一个成功的迹象。然而,耦合方案被认为是计算昂贵的,重点转移到寻找一种有效的顺序或分离的方法来求解控制方程。此时,Harlow及其同事[53]已经证明了交错网格算法的优势,Chorin[54]发表了关于压力分裂(或压力投影)算法的研究成果。Patankar和Spalding[55]将这些思想与早先开发的FVM技术和混合差分相结合,得出了一种鲁棒且通用的三维Navier-Stokes求解器,称为SIMPLE(压力链接方程的半隐式方法)。SIMPLE算法取得了经典的地位,并引领了新型CFD技术的商业化。
1970年左右,有人努力开发三维边界层(或三维抛物线流)的程序。人们很快意识到,这不是二维边界层过程的简单或直接的扩展。三维边界层在顺流方向上表现为抛物型,在横流平面上表现为椭圆型。对于这种“椭圆”问题,不能使用流函数-涡度过程,因为流动的三维性排除了使用流函数作为变量。显而易见的选择是直接使用速度分量。
Spalding的突破性贡献在于发现了压力梯度的作用。为了保持该抛物型问题求解过程的行进性,Spalding要求在给定的x点上,由横跨整个边界层的单一压力梯度dp/dx驱动流向动量方程。该压力梯度可以从外部自由流的压力梯度(对于无约束流动)或总体质量守恒(对于受限流动)中获得。然后,横流速度分量由其动量方程和横流平面上满足局部连续性方程所需的压力变化来控制。虽然计算了横流平面上的压力变化,但它没有用于修改用于流向动量方程的流向压力梯度dp/dx。因此,在处理压力时故意接受了一种不一致。Spalding坚持认为,这种不一致是使用推进过程的抛物线流的好处所付出的一个小代价。在实际应用中,横流压力变化很小,不会产生明显的误差。
如果这是一个重要的问题,为什么我们在研究二维边界层时不知道呢?二维边界层的解是否也包含关于压力的不一致?Spalding指出,这个问题一直存在于二维边界层中,但仍然是不可见的;正是三维边界层把这个问题带到了风口浪尖。
在二维边界层中,我们通过求解流向动量方程获得流向速度,其中我们使用单一的dp/dx值。然后由连续性方程求出横流速度。横流动量方程根本没用!如果我们使用它,我们将能够计算横流方向的压力变化,从而暴露与使用的dp/dx值的不一致。
在三维边界层中,有两个横流速度。连续性方程本身不能决定这两个速度分量。水流在两个横流方向上的分布是由横流平面上的压力变化决定的。我们必须求解横流动量方程,并计算横流平面上的压力变化。压力的不一致是藏不住的!Patankar和Spalding[55]在文中描述了压力链接方程的半隐式方法(SIMPLE)。该方法的主要成分是:
通过求解沿横流平面dp/dx值不变的沿流动量方程来计算沿流速度分量。
通过求解横流动量方程和连续性计算横流速度分量和横流压力变化。
后一步需要开发一种新的算法,这是SIMPLE的核心。人们认识到我们总能解出任何特定压力场的动量方程;但是得到的速度场通常不满足连续性方程。然后,我们可以对指定的压力场提出“修正”,使得到的(修正的)速度满足连续性。这导致了所谓的“压力修正方程”的诞生。它实际上是连续性方程的一种重写形式,它提供了一种直接计算压力修正(以及压力)的方法。
尽管在1972年的论文中,Patankar和Spalding在三维抛物流的背景下提出了SIMPLE算法,但很明显,该方法中包含了二维椭圆流的程序,因为这就是处理横向流速度的方法。此外,该方法可以很容易地扩展到三维椭圆流,因为它不依赖于任何特定于二维情况的概念(如流函数)。因此,SIMPLE为求解任何维度的耦合动量和连续性方程提供了一个通用的计算程序。再一次,这种方法在世界范围内被广泛使用。该方法的改进变体,如SIMPLEC[56]、simple[57]和simple[58]后来由不同的研究人员开发,部分是为了加强速度和压力之间的耦合,并提供更快收敛的更鲁棒的算法。
在20世纪70年代中后期,Spalding开始将类似simple的方法扩展到多相流。InterPhase Slip Algorithm (IPSA)是一种两相欧拉-欧拉算法[59,60]。它是单相SIMPLE算法的自然而优雅的扩展。虽然在其最初的公式和下面的描述中,只确定了两个阶段,但所有的主要思想都很容易扩展到多相(即两个以上)流程。它首次出现于1977年在杜布罗夫尼克的ICHMT研讨会上,见Spalding[60]。
Spalding认为两相流是由两个相互穿透和混合的连续体组成的。每个位置(因此,每个有限体积电池)都可能存在两种相(例如液体和气体)。各相的量以其体积分数表征。这种两相流模型没有描述两相在局部的分布情况。体积分数表征不能区分大块液体和细小液滴。正是这种特征导致了这个公式的简单性,也导致了它的主要弱点。Spalding后来提出了消除这一弱点的扩展[61]。
对于两个互穿连续体,他提出了两组相的守恒方程。除了多了一对相间输运项外,这些方程与单相方程相似。因此,两相在某一点处的速度差会导致由于一相对另一相施加剪切力而产生的力项,两相之间的温差会导致相间传热。相间传输在很大程度上取决于可用于剪切或热传递的表面积。这反过来又要求颗粒相液滴的大小。体积分数表征不能提供这些信息。事实上,多相流问题是封闭问题的一个经典例子,其中湍流可能是工科学生最熟悉的例子。
用分离的方法求解两相的守恒方程似乎是可行的;也就是说,先求出一个相的速度和温度,然后求出另一个相的速度和温度。因此,在求解一个相位时,假设另一个相位的速度和温度值是已知的。当相间输运占主导地位时,这一过程将非常缓慢地收敛,因为另一相的假设值将主导当前相的行为。Spalding已经在另一种情况下看到了这个问题。在本文其他地方讨论的热交换器分析中,壳侧流体和管侧流体的温度以同样的方式相互影响。因此,Spalding发明了部分消除算法(PEA)来精确地处理这种情况。诀窍是用代数变换的形式来解温度方程,消除了“其他”温度的直接影响。(仍然有一些项间接包含“其他”温度的相邻值,这就是为什么转换是部分而不是全部。)减少了“其他”温度(或速度)的影响,迭代过程得到了令人印象深刻的加速。Spalding立即将两相流问题确定为PEA的候选问题。结果表明,PEA的使用对于大多数两相流计算的成功是必不可少的。还施加了“求和统一”要求,以确保在溶液过程中始终保持相体积分数之间的耦合。
如果将整个连续性方程写成质量守恒方程,则密度较大的流体(液体)的项过多地支配了密度较小的流体(气体)的项。为了克服这一困难,Spalding写了一个“体积连续性方程”[62],该方程赋予了两相相当的权重。气液分析仪(GALA)采用体积连续性而不是质量连续性,有利于计算含有密度不连续的流体的运动,特别是那些大量级的流体,如气体和液体之间的运动。这种整体的连续性方程是建立两相流压力或压力修正方程的基础。虽然单一“共享”压力的概念是该程序的基础,但IPSA允许存在两种不同的流体压力,只要已知两种单独压力之间的适当代数关系。因此,相位连续性方程用于获得单个相的体积分数,并与SIMPLE类似,将整体连续性转换为压力校正方程。
与SIMPLE一样,IPSA迅速被世界各地的研究人员采用,IPSA的变体构成了当今几乎所有多相有限体积CFD代码中两相欧拉-欧拉公式的基础。两相流的基本框架后来使用“代数-滑移法”(有时称为混合模型或漂移-通量模型),从概念上扩展到各种多相流。本质上,滑移代数模型由混合物的连续性和动量方程、分散相的连续性方程和滑移速度的代数关系组成。我们可以考虑不同大小范围的液滴,不同浓度的气体包,等等。随后在PHOENICS代码中实现IPSA时,相位连续性方程中也加入了相位扩散项。Spalding用梯度扩散假设对这一项进行了建模,以表示与波动速度和体积分数之间的相关性相关的湍流通量。后来一些研究者主张在动量方程中排除相扩散项,而采用湍流色散力,如Lahey等[63]。在Spalding的双流体湍流模型Malin和Spalding[64]的一些变体中,相扩散项被证明是一个重要的元素。后来(即1980年后),Spalding利用两相概念设计了先进的湍流模型[64-67](考虑湍流是层流和湍流块的两相混合物)和先进的燃烧模型[68-70](考虑反应气体和未反应气体的混合物)。
随着二维薄剪切流的Patankar和Spalding[42,44,45]有限体积求解器的发展,需要(就湍流建模而言)转移到提供在整个计算域中确定有效湍流应力和热通量的途径,Spalding将这一问题简化为获取局部湍流粘度的问题。Patankar-Spalding方法中的初始模型是Prandtl的混合长度假设,该方案在早期的几篇论文中使用,包括他的小组提交给斯坦福会议的Kline等[71]。
然而,Spalding的设想已经在为循环流动提供相应的计算程序,在这里他认识到混合长度方案是相当不充分的。在文献[72]中,他着手解释了为什么在分离流中,壁面热流密度通常在再附着点达到最大值。正如他所说的那样[72]
“回想起来,在再附着点,(时间平均)剪应力为零,这一事实就显得特别有趣;大多数计算传热率的公式都是基于边界层中流动的物理理论,它们会预测,当剪切应力为零时,传热也必然消失。”
他解释说,分离流动的明显矛盾行为是由于湍流能量k向壁面的扩散输运。他通过对湍流能量方程的近似解析解,用经验系数证明了这一点,为这一论断提供了证据,从完全不同的流动类型中选择的扩散和耗散项,正确地再现了观察到的斯坦顿数和雷诺数之间的幂律依赖关系。在过去的二十年里,湍流能量方程已经被几个湍流建模者所采用,如Prandtl[73],但Spalding的论文是第一个展示了将该方程应用于分离流的理想化,以及它如何解决了为什么分离流中的传热遵循不同于边界层中的行为规则的问题。大概正是从这项研究中,人们决定将Prandtl-Glushko湍流能量模型嵌入到稍晚的再循环流数值方案中,Gosman等[40])。
第二个主要贡献,与上述相关,涉及改进近壁对数速度剖面(外部区域的数值解固定在其上,以避免在当时不可能实现的数值要求-更不用说在椭圆求解器中求解粘性子层的建模问题)。如上所述,Spalding认识到通常的公式
以及相应的温度分布
在驻点附近有几面的不足,因为在那里,摩擦速度是一个不合适的速度标在证明了分离流计算中包含湍流能的必要性后,他提出将上述壁面定律公式中出现的无量纲参数改为:
星号量引入了湍流能量k(在平衡状态下相同,相同,但并非相反)。
方程的改变:(12)、(13)至式。(14)和(15)只是简单地用代替。事实上,系数Cμ被选为0.09,这样,在壁面附近的平衡流动中(壁面剪切应力与“对数律”区域的湍流能量成正比),公式给出的结果与更传统的形式相同,方程(12)、(13)。
Eq.(14)的另一个重要特征是壁面剪切应力与壁面平行速度成正比,这在流动接近分离时很重要。Spalding版本的墙律从20世纪70年代初开始在帝国理工学院的软件中使用,并最终在Launder和Spalding的评论论文中发表[74],下文进一步参考。在过去的30年里,这种形式一直是大多数商业CFD规范的标准墙体处理方式。与原来的壁面定律不同,它总是给出至少在定性上正确的循环流动中的传热模式,尽管目前在研究级CFD中采用了更可靠的壁面函数方案。
尽管引入湍流能量来计算再循环或其他远离平衡的流动具有非常实际的好处,但Spalding意识到,由于prandtl- Glushko模型(像混合长度假设一样)需要湍流长度尺度的经验规范,因此它的适用范围极其有限。他觉得这更普遍,更完整的模型除了需要一个k的方程外,还需要一个进一步的输运方程来计算流场上湍流长度尺度的分布。受Kolmogorov[75]和Rotta[76]的工作启发,他们分别提出了波动频率ω和k与长度尺度乘积L的长度尺度方程(即kl方程),Spalding推动了使用该方程的模型的发展,即双方程模型的发展。在早期的报告[77]中,他使用Rotta的积分形式的长度尺度方程表明,可以用同一组常数计算各种自由湍流剪切流(特别是混合层、平面射流和风扇射流),而混合长度模型则需要不同的混合长度与流宽之比来给出正确的传播速率。然后,为了通过数值计算来支持这些发现,并为他的CFD团队创建的计算机代码开发一个通用模型,他让两名博士生在Patankar-Spalding边界层代码的帮助下,开发和测试一个使用Rotta长度尺度方程、k的输运方程和prandtel - kolmogorov涡流粘度假设的模型。Rodi负责自由流动模型的开发和测试,Ng负责近壁流动模型的开发和测试。这一研究导致了k-kL模型的建立。在Rodi和Spalding[78]中,对于自由流动,证实了该双方程模型优于混合长度模型。这项工作还提出了所谓的平面/圆射流异常,即在停滞环境中平面和圆射流不能用相同常数的模型来计算。Ng和Spalding[79]的近壁计算表明,如果在kL长度尺度方程中加入额外的一项,以产生正确的近壁长度尺度衰减,则该模型在这些流动中也表现得相当好。
Brian Spalding在紊流燃烧模拟领域做出了许多非常重要的贡献。的确,许多当前流行的关键思想在燃烧模型中已经出现在他的开创性发展中,包括混合在分子水平上的经常控制作用,以及火焰前沿表面现象所起的关键作用的问题的多尺度性质。用Spalding自己的话来说, [92]:
“很少有画家发明新的色彩,它们在画布上的排列构成了新奇。”
Spalding在紊流燃烧方面被引用最多的贡献可能是他1971年在第13届燃烧研讨会上发表的论文[93]。根据当时可用的湍流、预混、密闭、崖体稳定火焰的实验观察,他反思了火焰蔓延角度与实验条件(如燃料类型、等效比、接近流的平均速度或温度和压力水平)的相对独立性。根据这一证据,他得出结论,在很大范围的火焰条件下,化学动力学确实是辅助于空气动力学过程的。在1967年发表的一篇更早的前cfd论文中[94],Spalding分析了假设新鲜混合物在被夹带进入剪切层时立即燃烧的含义(当时称为“夹带-燃烧”模型)。当时提供的假设导致了正确的方向;但是,从多年后第一个CFD计算提供的有利条件来看,它被认为是粗糙的。因此,在1971年的论文中,夹带燃烧模型被“混合燃烧”模型所取代,该模型后来被称为“涡流分解”(EBU)模型。EBU的基本特征是,湍流燃烧的速度是大块新鲜气体分解成小块气体的速度,从而为分子过程的作用创造了界面区域;由于燃烧主要发生在这个界面区域,在足够大的Damköhler数下,假设燃烧速率可以通过求助于涡流破裂速率来模拟。
1971年的开创性论文[92]讨论了湍流剪切流中的火焰,其中,使用平衡假设,可以估计动能衰减率等于产量。因此,在混合控制下,单位体积的反应速率设为:
式中τ为混合物的局部反应度(或反应进度)。为了适应化学动力学是速率控制过程的情况,我们将阿伦尼乌斯速率与这个混合速率结合在一个谐波平均值中。
该模型在1975年发表的一篇论文中得到进一步完善[95],这是对燃烧室进行全面、多维建模的首次尝试之一。计算完全是椭圆的(与之前实验室火焰的主要抛物线模拟相反);紊流采用双方程k-ε模型;并建立了辐射的四通量模型。利用k和ε预测的优势,将涡流分解模型重新塑造为现在广泛使用的一种形式,即:
其中g是标量(混合分数)方差[81]。
如上所述,许多研究人员使用了几十年的涡流破裂模型(有一些次要的、非基本的变化),Spalding在很大程度上认为它是第一个粗糙的尝试[93,95],用于解释火焰空气动力学对燃烧速度的经常控制的影响,以及界面-表面形成和分子混合微观尺度现象的相关性。原始EBU模型在随后的研究中被Spalding和他的同事完善。因此,在[96,97]中,他概述了紊流燃烧的“滑动-拉伸”模型,其中“连贯的气体体在穿过火焰时被挤压和拉伸”。这一思想在[98]中得到了进一步阐述,并随后被赋予了ESCIMO模型的特点。
ESCIMO代表吞噬、拉伸、相干、互扩散和运动观察者,即模型的主要组成部分。Spalding借鉴了许多先前作者的思想,制定了一个模型[100],该模型可被视为适合在CFD代码中体现的湍流燃烧的第一个多尺度解释。对于EBU,假定化学反应发生在分离反应物(扩散火焰)或反应物与生成物(预混火焰)的界面区域。界面区域被排列成“褶皱”、“包裹”、“三明治”或“层”(图2);它们形成于流体的主要剪切区;它们在对流时伸展;然后当它们被其他漩涡重新吞没时“死亡”。因此,亚太经社会的任务是双重的:第一,计算褶皱的统计数字(即褶皱大小或宽度的局部分布),这被称为“人口部分”[100];其次,计算每个褶皱内的(一维)属性分布(这被称为“传记部分”)。第一种是欧拉的观点;第二种是拉格朗日,因为瞬态的一维轮廓依赖于(排除简化)流过流时的褶皱历史(这是模型名称的移动观察者部分)。另外,前者描述的是宏观尺度上的流动,对流和拉伸发生在宏观尺度上,而后者关注的是微观尺度上的流动,分子扩散和化学反应发生在微观尺度上。Spalding及其同事将ESCIMO模型应用于几种典型构型,如搅拌反应器、挡板火焰和自由射流[100,101],可以通过适当的简化假设计算出折叠种群的解析解。后来,它成功地应用于湍流,非预混氢-空气火焰[100]。
尽管ESCIMO在一个框架中体现了后来在广泛接受的模型中发现的关键概念,但整个科学界对ESCIMO的接受程度并不高。这个框架使得用20世纪70年代末的计算机可以解决这个问题。因此,非预混燃烧的小火焰模型使用了薄反应层分布的相同概念,其参数化使用的不是长度尺度(折叠宽度),而是时间尺度的倒数(耗散率)。褶皱的传记部分的时间依赖性(例如,层内的物种一维轮廓,图2)在所谓的稳态小片中被免除。这是局部计算的,与包裹(或小火焰)历史无关。Kerstein[102]最近独立导出的线性涡模型(LEM)显著地显示了ESCIMO的一些共同建模原则。值得注意的是,LEM还采用了嵌入三维计算的一维计算,以解决三维网格无法解决的小尺度混合和化学反应之间的耦合问题。与Spalding的ESCIMO不同,小尺度是在欧拉框架中计算的;通过一维场的随机搅拌来模拟湍流混合的效果,而不是采用局部的“涡流分布”。
欧拉多相方法的发展,如上文讨论的IPSA[60,103],为湍流燃烧建模增加物理真实感开辟了新的机会。Spalding很快认识到,在火焰前缘,燃烧(燃烧)和未燃烧的物质混合在一起,就像“陆地和泻湖景观”一样[68],可以建模为两相流,其中一相代表燃烧流体,另一相代表未燃烧流体。因此,在这个框架中,燃烧主要被视为发生在包裹边界上的界面现象,并由新鲜流体夹带到燃烧的流体中引发。燃烧速率主要由夹带速率决定,这需要建模。夹带模型通常涉及每单位体积的碎片表面积,这是长度尺度的倒数,可以解释为特征碎片尺寸;对于这种碎片大小,已经提出了几种代数和传输模型。燃烧流体内部也提供了燃烧的动力学控制。有趣的是,Magnussen[104]的单流体涡流耗散概念也有类似的概念,其中的“流体”(使用Spalding的术语)是“周围环境”或平均混合物和“精细结构”;
三维抛物线过程的早期应用之一是有曲率的管道。1972年,在抛物计算程序发表后不久,Pratap5[120]在Spalding的监督下,尝试计算螺旋盘管中的流动和传热。这项工作包括对发展中地区和完全发达地区的研究,完全发达地区是空间渐近状态。在热交换器、河流弯道、弯曲动脉、管道弯头等各种应用中,盘管中的流动具有重要意义。当流体流过弯曲的管道时,管道的曲率诱导二次流,在中心区域向外流动,在近壁区域向内流动。这些二次流被称为迪安涡[121,122]。在发展中区域,流向由抛物线型演变为外侧有峰的畸变型。这是动量从弯曲的内部区域转移到外部区域的结果。
④ V.S.Pratap目前被称为S.P万卡。
以上两个工作考虑了大的半径比。然而,当管道曲率较大时,抛物线假设就会失效。Spalding提出了一种新的过程,称为“半椭圆”或“部分抛物线”过程,它解释了加速流向流动的横向压力变化。Pratap和Spalding[128]跟进了这一想法,他们证明了捕捉这种效应的可能性。将该方法应用于强曲率导管内的流动,并采用Pratap[129]的实验数据进行比较。
绕垂直于管道轴的轴旋转的管道中的流动呈现出另一种情况,物体对流体施加力。Majumdar等[130]应用三维抛物线计算程序求解控制三维流动方程的抛物线形式。动量方程包含由于离心力和科里奥利力而产生的额外项。在随后的研究中,Majumdar和Spalding[131]将部分抛物线过程应用于旋转管道中的流动。与数据的一致性总体良好;然而,观察到某些差异。该算法认识到三维标量压力场确实是椭圆形的,并且在流动的横平面上存在流体再循环,这就需要将稀缺的计算机内存分配给压力和速度分量,但由于在主流方向上没有从下游到上游的流动,因此只需要将压力存储在内存中,因为流动基本上是抛物线流动。Pollard和Spalding[132]采用了这一思想,发展了部分椭圆方法。
Majumdar等[133]研究了圆形风管和围绕平行于风管轴线的水平轴旋转的方形风管完全发达区域内湍流流动的局部和平均换热特性。结果表明,旋转增加了换热。结果与实验数据吻合较好。Majumdar先前的一项研究[134]也考虑了绕垂直轴旋转的径向扩散器。
部分抛物线算法的发展证明了Spalding对流体力学和一般输运现象的掌握。上述关于具有体力的流动——受旋转和曲率影响的管道流动——的工作是开创性的。新开发的部分抛物线法(以及抛物线法)已成功应用。这导致了一套新的计算程序来预测复杂的三维流动,到1975年,帝国理工学院的传热部门正在考虑这种流动(“复杂流动”当时被定义为具有外部施加的应变速率,或者两个或多个基本湍流相互作用[135],而不是“基本流动”:射流、尾迹、通道、边界层)。在此期间的研究重点之一是建立使用湍流模型预测三维湍流的能力,这在本文的其他地方进行了讨论。在许多情况下,还需要实验数据,如果设备在世界其他地方存在,就把学生送到那里。回顾1970-1980年期间该小组发表的约75篇论文,可以发现大量涉及Spalding的主题,例如:自由和强制湍流对流,Abdelmeguid和Spalding[136];通过和绕过各种形状和障碍物的流动,Caretto等[137]燃烧,Serageldin和Spalding [138];低普朗特数传热,ElHadidy等[139],t形结流动,Pollard[140];Majumdar等[141]。许多这类研究的一个显著特点是将计算方法和模型的发展结合起来,这些方法和模型能更好地捕捉所研究现象的主要物理特征。
工程师通常使用什么数值方法来计算固体物体的应力和应变?通常的答案是“有限元法”(FEM),它通常应用于结构分析问题的模拟,并形成许多计算机代码的基础。然而,有许多多物理场问题,在这些问题中,固体应力必须同时计算,例如,湍流流体的传热、传质、化学反应等等。这种流固相互作用(FSI)问题具有很大的科学意义和实际必要性。例如,旋转涡轮机中的热应力和离心应力,人体动脉中的血液流动等弹性管中流体流动中的压力波传播,或包含压电膜的打印机腔中的流动,等等。在这些和其他情况下,方法的选择不太清楚,因为与有限体积方法(FVM)相比,FEM算法对于流体动力方程耦合系统的模拟可能没有那么复杂,有限体积方法(FVM)发展得很好,目前几乎所有主要的CFD代码都采用了有限体积方法。此外,此类CFD代码已经包含了许多物理现象(湍流、辐射、化学反应、多相流等)的子组件模型,这些模型是模拟复杂FSI问题的必要组成部分。
Spalding是最早提出FVM可以很容易地用于求解弹性问题的人之一,并且,建议在基于FVM的代码中开发方程可能比将现有的大量物理流体力学和传热模型从已建立的基于FVM的CFD代码移植到基于FEM的代码中更简单。Beale和Elias[148,149]展示了如何将FVM代码PHOENICS应用于应力分析问题,利用蠕变流体流动和固体力学中的流函数和应力函数之间的类比。该方法在理论上类似于Runchal等人[43]的φ-ω方法,因此仅限于二维平面应变/应力问题。Spalding随后考虑了位移,而不是力的类比,后者在三维空间中具有相当普遍的主要优势。
线性热弹性的控制方程以平衡方程为基础,结合广义Hooke定律和应变-位移关系,得到Navier位移方程。Spalding最初将位移方程考虑为如下形式:
毫无疑问,Spalding和他的同事们在1969年创立了一家名为燃烧热能和动量有限公司(CHAM)的公司,开创了商业CFD行业。该公司于1974年更名为ConcentrationHeat and Momentum Ltd.,提供商业CFD服务,该服务是基于Spalding在帝国理工学院的研究小组提出的开创性技术。虽然CHAM的活动和CFD代码PHOENICS的起源(见下文)不可避免地交织在一起,但早在PHOENICS出现之前,CHAM基本上是世界上唯一的商业CFD公司,直到20世纪80年代初,Creare公司开始销售FLUENT代码,其早期起源可以追溯到Swithenbank等人[152];参见Runchal[153]。在1974年至1980年间,CHAM为广泛的工业领域的工业客户开发了许多特定应用的CFD计算机代码,包括航空航天,汽车,国防,化学,环境,消防和安全,船舶,制造和工艺,核电和化石燃料发电行业。这些CFD规范(其中一些进入了公共领域)也为中国商会承接CFD咨询合同提供了手段,旨在为行业解决实际问题。CHAM的客户大多是英国、欧洲和美国的大型工程公司。在此期间,在Spalding的指导和监督下,中国商会取得了许多“第一”。例如,它为核工业的蒸汽发生器创建了第一个CFD代码:为此,CHAM必须设计能够模拟同一空间内蒸汽和水的分离流动的CFD程序。另一个例子是,CHAM为铝生产行业编写了第一个电解冶炼厂(ESTER) CFD代码,用于模拟霍尔电解槽型的多阳极电解冶炼厂,其中熔融金属/电解质界面上的电磁效应和重力波可以相互作用,从而限制设备的性能。后来,ESTER被基于PHOENICS的产品Rosten所取代[154],该产品也被行业广泛使用[155-157]。煤和燃气炉、热交换器、柴油和汽油(汽油)发动机、动力冷凝器、蒸汽发生器和自然通风冷却塔都是CHAM生产的第一批CFD软件包的设备项目。
在20世纪70年代,前面描述的部分抛物线方法为一些3D工业应用提供了比全椭圆代码更实用的路线。以部分抛物型CFD代码FLASH为例,采用曲线坐标系计算船舶和潜艇壳体的三维绕流。这项工作主要由Ishikawa-HarimaHeavy Industries (IHI)(日本)、NationalMaritime Institute (NMI) (Teddington)和AdmiraltyExperimental Works (AEW) (Haslar)资助,Spalding及其同事发表了大量出版物,例如Abdelmeguid等人[158]和Markatos等人[159]。随着计算机技术的日益强大,FLASH被内部的椭圆型CFD代码所取代,然后最终被基于PHOENICS的CFD模型所取代,该模型采用一般的非正交体拟合坐标系,Malin等[160]。部分抛物线概念在涡轮机械应用中也被证明是有用的,1976年,Singhal和Spalding[161]开发了SEMFT部分抛物线CFD代码,用于预测涡轮机械叶栅中的二维稳定湍流。这个代码是1977年提供给苏尔泽兄弟公司的,它能够处理不可压缩、亚音速、超音速或跨音速的气流。CATHY3代码是另一种部分抛物线式涡轮机械代码,提供给许多公司,包括1978年的西屋电气公司。采用曲线坐标系对旋转涡轮机道内的三维定常湍流进行了预测。该代码来源于Pratap和Spalding[128]以及Majumdar和Spalding[162]的早期工作。1979年,它被同名的全椭圆代码Malin等人[163]所取代。该椭圆型CFD代码包含一个空化模型,该软件被提供给NASA Marshall,用于离心叶轮内稳定的三维湍流预测。
如前所述,Patankar和Spalding[142]率先将分布阻力概念应用于壳管式换热器三维流动和传热的CFD建模中。自20世纪70年代初以来,该方法已被CHAM广泛用于进行热交换器,蒸汽发生器和冷凝器在国防,船舶,核能和电力工业中的应用的CFD模拟。这些应用是由许多公司资助的,包括燃烧工程 (美国) 劳斯莱斯(英国)、国家核能公司(英国)、电力研究所(美国)、欧洲原子能共同体(意大利)、北方工程工业帕森斯(英国)、通用电气公司(英国)和大瀑布(瑞典)。Spalding及其同事通过这些研究和其他研究发表了大量论文,例如Singhal等人[164,165]、Al-sanea等人[166]、Jureindini等人[167]和Malin[168]。一些应用包括通过IPSA进行均匀、代数滑移或双流体建模的两相流模拟。最终,PHOENICS成为了CHAM所有此类计算的载体,但著名的前PHOENICS CFD代码是由Spalding和其他CHAM人员开发的ATHOS,CALIPSOS, FOCS, LOBSTER, STELLA和URSULA。
许多工业多相计算的起源可以追溯到GALA[61,62]和IPSA[60,61],两者都是由Spalding在70年代中期设计的。这些在之前的手稿中有描述,但这里应该提到的是,Spalding的杰出论文[61]也描述了不同的界面跟踪技术,作为使用GALA模拟自由表面流动的替代方案。一种方法涉及跟踪粒子,而另一种方法涉及求解液体高度的守恒方程。在20世纪70年代末和80年代初,CHAM和帝国理工学院的团队都成功地采用了界面跟踪技术,例如,参见Maxwell[169]、Awn[170]和Castrejon和Andrews[171]。IPSA对工业多相流的数值计算产生了巨大的影响。它构成了CHAM及其客户使用PHOENICS和早期CHAM代码(Moult等人[172]、Baghdadi等人[173]、Markatos[174,175], Marchand等[176],Al- sanea等[177],Boisson和Malin[178]等人)进行的大部分两相流模拟的基础。不可避免地,这种方法后来在商业CFD产品中找到了竞争的途径。
Spalding满足于在许多合适的多相流应用中使用代数滑移模型。该模型在20世纪70年代后期在CHAM用于蒸汽发生器应用,Singhal等[179,180],但在20世纪80年代,它的使用变得更加广泛,用于化学和材料加工行业。在此期间,使用了几种模型的变体,主要限于考虑重力和/或离心力场,例如:Pericleous和Drake[181](旋风中的颗粒分离);Tacke和Ludwig[182](包体分离);Pericleous和Patel[183](搅拌槽式反应器中的固体悬浮物)。最后,Spalding提出了一个更一般的公式,即使用GALA来满足体积连续性。该方法于1992年在PHOENICS实施。近年来,使用代数滑移建模的多相能力开始出现在其他通用CFD软件中。
在20世纪70年代中期到80年代初,另一个新的应用领域是建筑物中火灾和烟雾运动的预测。当时,英国Borehamwood消防研究站(FRS)资助了CHAM在2D和3D椭圆MOSIE CFD代码开发和应用方面的开创性工作,Markatos等[184]。这些特殊用途的火灾和烟雾模拟代码被基于PHOENICS的模型所取代,该模型随后被FRS和建筑研究机构(BRE)使用,以创建著名的CFD火灾模拟代码JASMINE, Markatos & Cox [185], Kumar等[186],该代码使用PHOENICS代码的早期版本作为其CFD引擎。通过这种方式,Spalding的贡献可以看作是为FRS和BRE提供了基本的CFD能力,使他们能够开发和验证JASMINE中包含的CFD技术,使CFD现在成为消防安全工程中公认的组成部分。
在CHAM的早期,燃烧建模是业务的重要组成部分,到1973年,该公司为Rocket Propulsion Establishment(Westcott, UK)创建了一个非常先进的椭圆CFD燃烧代码,称为BAFL,参见Tatchell和Spalding[187]和Jensen等人[188],用于预测以超音速或亚音速移动的火箭下游区域的流动、燃烧和传热。BAFL走在时代的前面,它结合了许多复杂的数学模型,使其能够处理可压缩的、湍流的、气体之间双向化学反应的再循环流动,温度和浓度的波动,特定波长波段的辐射传热(Lockwood和Spalding[46]),以及与不同尺寸范围的颗粒物质的反应(Spalding[189])。事实上,CHAM出售给燃气轮机制造商的第一批CFD代码中,有一些是专门用于燃烧室模拟的。这些燃烧器代码中有几个已经使用了十多年,当然也得到了用户的进一步发展;并且其中至少有一个通过美国陆军进入公共领域,使竞争的CFDcode供应商能够开展业务;他们很快就照做了。1976年,CHAM向NASA Lewis研究中心提供了一个燃烧CFD代码,该代码使用化学平衡模型来表征碳氢化合物反应,并使用化学动力学模型来预测氮的氧化物,Elghobashi等[190]。这可能是预测污染物形成的第一个CFD模型。20世纪70年代另一个值得注意的CHAM燃烧规范是CORA3,它被交付给了许多工业客户。除其他外,该CFD燃烧代码利用了Spalding的“eddybreak”[93191]和“假设PDF”方法[81191],这些方法考虑了湍流波动。Fiveland提到了Spalding的工作对CFD在发电行业的影响,参见Hanna[192],他于20世纪70年代在美国的Babcock和Wilcox开始了CFD在锅炉制造业的研究工作:
“当时我知道Spalding教授在英国开创性的CFD工作,我们看到了使用这种‘新’CFD技术改进我们的锅炉和燃烧器设计的价值。”
20世纪70年代的计算机内存小,执行速度慢,因此需要巧妙的编程来解决工业问题。此外,在一些圈子里,对于计算机模拟复杂的湍流和多相流是否会有任何实际用途,人们持怀疑态度。1978年的某个时候,Spalding设想了一个能够处理所有流体流动过程的单一CFD代码的想法。因此,CHAM放弃了开发单个应用特定CFD代码的政策,并在1978年底开始创建世界上第一个通用CFD代码:PHOENICS。
许多科学家和工程师在他们职业生涯的某个阶段,通过撰写他们专业领域的专著而脱颖而出。然而,Spalding将他的时间和精力投入到一个实用的CFD软件上,这反过来又帮助催生了一个完整的行业领域,并创造了前所未有的就业机会和财富。因此,CFD作为计算机辅助工程的一个子集可以直接追溯到抛物线、双曲或椭圆型数值积分代码系列PHOENICS的发展[51,193-195]。
PHOENICS于1981年作为FORTRAN-input程序商业化推出,后来在1984年作为交互式输入CFD代码推出。这是第一次对所有热流体问题使用单一的代码。这代表了一个重大的哲学变化。实际上,在那个时候,对于通用的CFD代码是可能存在的想法存在着很大的阻力。PHOENICS的成功也要归功于其他一些人,尤其是H.I. Rosten,如果不提及J.C. Ludwig自代码开始以来所做的重要和持续的贡献,那就太失职了。然而,PHOENICS本质上是Spalding的心血结晶,直到今天仍然是他关注的主题。PHOENICS代表了CFD发展的范式转变;工程师们现在能够将之前在帝国理工学院和其他地方开发的算法和封闭假设应用于大量现实世界的工业问题。一旦达到了这一步,CFD就不可逆转地走出了学术界的象牙塔,进入了工程产品和工艺的现实世界。在软件设计、内存管理/优化、可移植性、文档和培训方面,向客户提供工作软件产品的需求将专业水平提升到了一个新的水平,而这在大学环境中是不常见的。由于代码的通用性,CFD超越了热流体问题,发展到了一系列多物理场问题,例如,化学气相沉积。
在20世纪80年代早期,PHOENICS是唯一可用的通用CFD代码,它代表了Spalding小组先前发表的许多科学原理的体现,其中许多原理被精心地复 制到一个案例库中,供用户执行。为层流或湍流编写单相和多相代码,具有复杂的化学动力学,传热和传质(包括热辐射)和可变性质,这是非常独特的概念。该代码被设计为一个被称为“地球”的核心。这包含了所有已建立的硬编码和用户无法访问的实用程序(为了避免没有经验的人损坏,从而便于维护和支持)。还提供了开源组件,其中对物理模型(如湍流闭包、壁函数和求解算法)进行了编码,并且可以由用户检查、更改或替换。“钩子”(线程)在软件的不同位置提供,供用户附加自己的代码,从而允许扩展代码功能,并使其能够用于许多不是专门为其设计的问题。这种“开放源代码”的概念最早出现于1981年,后来被其他一些cfd代码供应商所采用。PHOENICS是基于有限体积法的。代码的核心是单相最简单算法[58],SIMPLE[55,137]的一种变体。使用IPSA[60](欧拉)和GENTRA(拉格朗日)来调节多相流。用户可以从许多线性方程“求解器”中进行选择,例如逐点,平板和整个场(或提供自己的),并在每个变量的基础上使用惯性或线性方法来控制松弛。最初采用的是交错速度方案[53],尽管后来在直线坐标和贴体坐标中添加了并列方案[196]作为选项。边界条件和源项本身被编码为线性源项,在式(11)中,其中G是几何因子,C是源项“系数”,V是源项“值”,用户可以对其进行实质性控制;因此它们是可编程的。代码基于隔离方案,对每个字段变量执行多次“迭代”,并对整个计算域进行多次“扫描”,直到达到收敛,求解器停止。点值和残差最初是打印出来的,在后来的代码版本中,由图形用户界面(GUI)动态绘制,为用户提供对结果的信心度量。
在PHOENICS的早期版本中,用户可以根据需要使用FORTRAN在预处理器' SATELLITE '中或在运行时' GROUND '中定制代码。通过GROUND,用户实际上可以通过执行各种操作的函数和子例程访问存储在单个“F数组”中的所有EARTH变量和数组。代码的性质允许广泛的可编程性,因此如果用户愿意,可以完全绕过所有默认的求解器等等,而支持他们自己的实用程序。需要强调的是,大量的学术和工业CFD计算都是基于PHOENICS代码进行的。1985年第一届国际PHOENICS用户会议记录[197]显示,从一开始就有大量的PHOENICS应用于非常广泛的行业。此后每两年举行一次的用户会议,连同PHOENICS计算流体动力学杂志和公开文献的出版物,将工业应用的数量扩大到真正巨大的比例。Spalding的学生在帝国理工学院(Imperial College)的计算流体动力学部门(Computational Fluid Dynamics Unit)[1988 - 215]和其他地方撰写了许多基于PHOENICS的博士论文。虽然现在研究生在他们的研究工作中使用商业代码的概念几乎没有争议,但在当时这是相当新颖的。
一些代码的发展包括:贴体网格(1983年);生物流动——膀胱、肺等(1987年);可编程输入语言和共轭传热(1989);以专家为基础的系统以改善趋同(1991);固体应力/应变(1992年);并行执行(1993年);细网格嵌入(1997)。后来,引入了新的编程功能,允许用户创建自己的物理模型,而不需要诉诸FORTRAN。随着用户基础的增长,案例库也在增长。多年来,工业界和学术界用户的无数发展都被采用到代码中。此外,特殊版本的PHOENICS被开发用于广泛的应用:电解冶炼厂,化学气相沉积,电子冷却,高炉,冷却塔,外部空气动力学。中国商会成立了PHOENICS专用用户组。
PHOENICS是世界上第一个通用CFD代码;到1981年出版时,基本的科学原理已经被广泛地发表了。因此,一旦CHAM证明了通用CFD代码可以被创建并在市场上盈利,竞争对手很快就会紧随其后。PHOENICS的应用和原理被大量模仿;现在市场上有许多这样的CFD软件包,由竞争厂商提供。大多数这些代码的血统可以追溯到中国商会本身或帝国理工学院团队,中国商会和PHOENICS的影响有时受到员工流动的影响。当然,可以说,中国商会在创建和营销PHOENICS方面的创新开创了一个基于CFD软件和服务销售的全球行业,尽管没有确切的数据,但目前该行业的年营业额已超过六位数。
在第一版PHOENICS发布的几年内,Spalding在CHAM发起了基于PHOENICS的特殊用途产品(SPPs)的创建和营销工作,这些产品采用了特定于应用程序的GUIs。其中第一种是ESTER,它是专门为铝工业开发的。其他spp在20世纪80年代末和90年代初迅速跟进,例如,HOTBOX(用于电子冷却),FLAIR(用于建筑物中的加热,通风和火灾和烟雾运动),TACT(用于冷却塔),PHOENICS-VWT(用于赛车空气动力学)和PHOENICS-CVD(用于化学气相沉积)。Spalding早在1979年就主张将CFD应用于CVD反应器[216],其中包括使用嵌入式精细网格来更精确地解析晶圆附近流体流动和沉积过程的精细结构的建议。
PHOENICS-CVD, Heritage[217]仍然是SPPs中技术最先进的,能够通过模拟多组分气体中的流体流动和传热来模拟各种CVD反应器的行为,包括气相(均相)和表面(非均相)化学反应、热辐射和等离子体效应。在20世纪80年代和90年代,PHOENICS被CHAM和其他组织大量用于冶金和材料加工应用。此类应用通常具有挑战性,可能涉及流体流动、传热传质、相变、电磁学、自由表面流动、多相和化学的组合。已经提到了在铝工业中成功使用PHOENICS-ESTER。此外,Alcan International与CHAM合作,在PHOENICS中使用IPSA模拟浸入式液态金属射流的混合和凝固,参见Enright等人[218]。CHAM还使用IPSA为NASA Lewis研究中心创建了基于PHOENICS的CFD模型来模拟二元材料的固液相变化,参见Prakash[219,220]。在20世纪80年代中期,CHAM还开始开发基于PHOENICS的HIsmelt过程CFD模型,该模型涉及通过将铁和煤粉直接注入熔融浴来直接熔炼铁矿石,以及随后在顶部空间释放和燃烧熔炼气体。20世纪90年代初,中国商会的员工流动促进了Rio Tinto(澳大利亚)、Davis等人[221]对该模型的持续改进。许多其他公司也开发了基于PHOENICS的材料加工CFD模型,包括Arbed Recherche(卢森堡)、Centro Sviluppo Materiali(意大利)和Hoogevens(荷兰)。在20世纪80年代和90年代初,PHOENICS在麻省理工学院也被Szekely及其同事广泛用于材料加工操作的CFD模型的开发,在此期间发表了大量论文,例如:Ilegbusi和Szekely[222,223](中间包和连铸),Cartwright等[223](Czochralski系统),Choo和szkley[224](弧焊)。Saluja等[225](搅拌熔体)。直到今天,PHOENICS还在瑞典冶金研究所(Metallurgical Research Institute AB)和墨西哥国立理工学院(National Polytechnic Institute)高级研究和研究中心(Centerfor Advanced Study and Research)使用了多年,用于进行金属加工研究,例如,参见Solhed和Jonsson[226]和Rodriguez等人[227]。
在20世纪80年代早期,用于表示复杂几何形状的笛卡尔切割细胞方法在CHAM用于特定应用,但直到十年后,Spalding才开始启动一项工作计划,为PHOENICS配备一种普遍适用的设备,称为PARSOL(部分固体)。其动机是隐式地处理导入的CAD几何图形,并大大简化与使用体拟合非正交网格相关的网格生成。尽管在Spalding[228]中进行了一些讨论,但PARSOL从未正式发表。该方法对流域的大部分区域采用了背景笛卡尔网格,并对被实体切割的单元进行了特殊处理,从而保留了符合边界的网格。PARSOL也是部分非结构化的,因为它允许通过嵌入的细网格区域对局部网格进行细化。该方法计算分数面积和体积,并采用了一系列特殊算法来计算界面面积、计算壁面应力、计算固体边界附近的平流和扩散等。笛卡尔切割单元法现在已经成为结构化和非结构化体一致性网格的公认替代方案,并且现在至少在其他三种商业CFD代码中使用。
在1993年,CHAM在PHOENICS中实现了一些采用正态壁面距离阻尼函数的低雷诺数湍流模型。人们很快认识到,使用传统的基于搜索的几何程序来确定壁距分布,对于处理PHOENICS工业应用中使用的任意复杂几何图形来说,计算代价太大。就在那时,Spalding提出了一种巧妙而又极其经济的方法,即通过求解标量变量的泊松型输运方程来计算壁面距离的近代值。Spalding的建议主要是直观的,基于对两块板之间热辐射的类比。这种基于微分方程的壁距算法立即在PHOENICS中用于所有相关的湍流闭包。该方法仅在第10届布莱顿国际传热会议的海报会议上由Spalding本人正式提出,但此后它被Tucker[229,230]等人利用,改进并发表,最终该方法或其某些变体进入了大多数商业CFD代码。
壁面距离法的发明与Spalding提出的零方程低雷诺数湍流模型LVEL同时出现,该模型用于确定流体流过具有许多固体物体的空间时的湍流运动粘度[231]。在这种情况下,附近固体之间的网格密度通常太粗,任何更先进的湍流模型,如双方程k-ε模型,都无法有意义地使用。紊流运动粘度是通过Spalding壁面定律计算的[232],该定律涵盖了整个层流和紊流状态。该模型的一个关键要素是前面提到的基于泊松的壁面距离微分方程。LVEL模型在PHOENICS中实现,并且在电子冷却应用中具有很大的实际优势,因此几乎不可避免地,该模型开始出现在其他商业CFD代码中,特别是那些致力于电子元件和系统的热设计的代码中。Spalding在他的一篇论文中写道[231]
“LVEL模型应该被视为为传热工程师的问题提供了一个实用的解决方案,而不是提供了一个新的科学见解。”
当然,该模型的实际用途已被许多工作者证明,包括Rodgers等人[233]、Dhinsa等人[234]和Choi等人[235]。后者报告说,LVEL模型与更先进的双方程模型一样有效,用于模拟机架式服务器,其最大的优势是“可以节省大量的计算时间(三倍或更高),特别是在进行瞬态CFD模拟或在不同机架设置下进行稳态运行时”。
Spalding认识到工业工程师需要在合理的计算时间内获得适当的解决方案,从而开发了IMMERSOL辐射模型。该模型由Spalding于1994年提出,主要用于经济地表示拥挤空间中的辐射,如电子冷却应用,以及在流体中分布许多导电固体的任何其他问题。IMMERSOL是浸入式固体(immersion SOLids)的缩写,CHAM将其应用于PHOENICS中,作为一种经济的选择,可以在光学厚度和厚度介质中表示复杂几何形状的辐射。Spalding尚未在公开文献中发表该模型,但Rasmussen[236]和Yang等人对其进行了描述。[237]。该模型涉及辐射系数的扩散方程的解,其中扩散系数与介质的吸收和散射系数以及固体壁之间的距离成正比。壁间距离是由泊松方程自动计算的,正如前面在湍流模型的壁距计算中所解释的那样。辐射方程是通过简化Hamakar和Schuster的六通量模型而设计的,参见Spalding[238],忽略了辐射的方向效应,但将其扩展到透明介质是基于物理直觉。IMMERSOL提供了一个经济可行的近似辐射传递方程,给出了两个无限长的平行板之间的热辐射的精确解决方案;在更复杂的情况下,它的预测总是可信的,而且是正确的数量级。Yang等人[237,239]和Rasmussen[236]都报道了该模型在燃烧应用中的成功应用,Hien和Istiadji[240]也报道了该模型在太阳辐射下的建筑通风研究中的成功应用。到目前为止,除了那些使用PHOENICS的组织外,IMMERSOL似乎对工业CFD影响不大。在Rasmussen[236]的工作中发现了一个例外,他通过用户编程将IMMERSOL实现到另一个商业CFD代码STAR-CD中,然后扩展模型以处理波长依赖性。
“Spalding能够将数学天赋与丰富的物理直觉相结合,这导致了一种科学方法,这种方法得到了在实际工业问题的“前沿”工程师的普遍赞赏,”和“商业CFD软件开发人员并不总是以他们的慷慨认可而闻名,因此有必要指出,实际上所有今天销售的商业代码都要归功于Spalding和同事的开创性工作。”
本文作者:
V. Artemov , S.B. Beale , G. de Vahl Davis, M.P. Escudier, N. Fueyo, B.E. Launder, E. Leonardi* , M.R. Malin, W.J. Minkowycz, S.V. Patankar, A. Pollard, W. Rodi, A. Runchal, S.P. Vanka(注:Eddie Leonardi 于2008年12月14日去世。)
来源:多相流在线