Ccond=0.000 1时小流量和设计流量工况下离心泵扬程随有效空化余量变化的非定常数值计算结果与试验结果的比较如图2.6所示。表2.2NPSHc计算值与试验结果的比较在本书的离心泵空化流动数值模拟中,质量输运空化模型的凝结项经验系数取0.000 1进行数值计算。图2.7列出了小流量下NPSHa=1.0m和设计流量下NPSHa=1.6m时,3个不同凝结项经验系数所对应的叶轮内空泡体积率分布。......
2023-06-15
Stone(1968)对全隐式差分方程提出了新的迭代解法——强隐式迭代解法,简称SIP方法。由于这种迭代方法利用解方程组对全部格点同时加以改进,迭代效率非常高,已经十分接近直接解法,因而也就更加有效。Trescott等(1977)所做的比较研究表明:对于复杂的实际问题,SIP方法明显优于其他的迭代方法。
1.强隐式差分格式
将一般的全隐式差分方程写为以下的形式(二维):
这里略去水头关于时间的上标n+1,与式(5-3-16)中每个系数相联系的各点位置见图5-3-2,共五个格点,称式(5-3-16)为五点格式。
对渗流区域内的每个格点都列出式(5-3-16)并利用给定的边界条件,就得到一线性方程组
暂时设渗流区域为一矩形,划分成N×M个差分网格(图5-3-3),格点号码按自左而右、由下而上的顺序编排,则由五点格式所确定的上式中的系数矩阵[K]具有下列形式:
图5-3-3 假设的矩形2D渗流区位置图
矩阵的第一行相当于对格点(1,1)列的方程,因为利用了边界条件,所以a1,1和c1,1都不出现,其他依此类推。可以看出,矩阵[K]的非零元素全部集中在五条对角线上。此外,它还是一个对称的矩形。矩阵[K]具有式(5-3-18)的形式是由格点号码的特殊编排方式所决定的。对于这样一个形状简单、高度稀疏的矩阵,在对它进行LU分解时,期望能分解成下列形式的[L]和[U]的乘积:
下三角矩阵[L]和上三角矩阵[U]包含的非零元素的位置类似把[K]沿主对角线分成两半。假若能做这样的分解,那么从[K]得到[L]和[U]的计算量一定是很小,因为[L]和[U]也是高度稀疏的。事实上对[K]分解得到[L]和[U]并不稀疏,使得计算量变得很大。
于是把问题倒过来,由式(5-3-19)和式(5-3-20)确定[L]和[U]的乘积,按照矩阵乘法的规则,有
式(5-3-21)的非零元素集中在七条对角线上,比[K]多了两条非零元素对角线。假若对每个格点列的差分方程不用五点格式而是在形式上用七点格式,如式(5-3-22):
与公式(5-3-22)的系数相关联的格点表示在图5-3-4中。由此形成的方程组设为
图5-3-4式(5-3-22)系数相关联的格点示意图。
图5-3-4 式(5-3-22)系数相关联的格点示意图
其系数矩阵为
2.系数矩阵的近似因子分解
比较式(5-3-21)和式(5-3-24),可直接得到[^K]的元素[L]、[U]的元素之间的下列对应关系
另一方面,根据“七点”格式与“五点”格式之间的关系可以找到[^K]与[K]的关系。事实上,利用Taylor展开式并略去高阶项得
因为
所以由式(5-3-27)可得下面的近似关系
将上式代入到式(5-3-16)中,得到
这样就把式(5-3-22)近似地化成了式(5-3-16)的形式。与式(5-3-16)比较得
式中:ω为松弛因子,目的是为了加快收敛速度。
在上式前添加松弛因子ω(0<ω<1)。
计算是按照从左而右、自下而上的顺序进行,则式(5-3-31)中的δi,j-1、δi-1,j、ηi-1,j、ηi,j-1应作为已知,因为它们已经在左、下格的计算中被求出或者可根据边界条件定出。因此,借助式(5-3-32),对每个格点(i,j)都能由“五点”格式系数ai,j、bi,j、ci,j、di,j、ei,j算出[L]和[U]的非零元素αi,j、βi,j、γi,j、ηi,j和δi,j。矩阵[L]和[U]是[^K]的因子分解,不是[K]的因子分解,但[]和[K]很近似。到这里已经实现了对系数矩阵[K]的近似因子分解。
松弛因子ω的确定通过以下公式:
其中
SIP迭代对ω的值并不敏感,ω一般大于0.92,在[0.92,0.95]之间(Peric,1995)。
3.强隐式迭代解法
解方程组[K]H=F等价于解方程组
这一方程组左端的系数矩阵[^K]已有LU分解,且[L]和[U]均为稀疏矩阵,因此很容易求解,但因式(5-3-34)的右端也含有未知的水头H,所以只能用迭代解法,将式(5-3-34)改写为下列形式
式中 H——右上角的指数表示迭代的次数,两端同减Hm,得
或者
其中
迭代的过程是:已知Hm,可由式(5-3-39)决定Rm,以此作为式(5-3-40)的右端项,解方程组求得ξm+1,再由式(5-3-38)解出Hm+1,如此重复迭代直到收敛为止。收敛的准则采用下式
其中ε是足够小的正数。
上述迭代过程主要工作是解方程式(5-3-37),[^K]已有LU分解,所以在解方程时很容易求解,只需一个前推回代的过程。
4.前推回代工程
式(5-3-43)和式(5-3-44)为中间值V和ξm+1的表示方法,这样就能顺利的解出[L]和[U],从而求解整个方程。
实际含水层区域不一定是矩形,但一定可以采用矩形区域将其覆盖,用两组平行线划分成矩形网格,在这些矩形网格中可分成内格点、靠近边界的格点和外格点三类。这样就可以运用SIP方法求解。
5.SIP算法验证
算例1
在200m×200m的区域上(浓度为0),北部和西部隔水,东部和南部边界有稳定的浓度为1mg/L,且地下水为稳定流动,孔隙平均流速为0.005m/d,流动方向平行于X轴并与X方向相反。ΔX=ΔY=20m,Δt=5d。总时段为2000d。共400个步长,DL=10,DT=1(采用陈崇希、李国敏等编著的《地下水溶质运移理论及模型》“求矩形含水层中地下水溶质运移二维ADI交替隐式差分算例”)。
区域剖分分辨率为20m×20m,共100个网格。剖分网格与初始浓度场图见图5-3-5。
图5-3-5 有限差分网格剖分图及初始浓度场图
表5-3-1和表5-3-2是时间为500天时分别用ADI算法和SIP算法算出的区域浓度场分布结果。可以看出两结果稍有不同,但是仅在小数点后第三位开始略有不同,并不影响计算精度。这些差别也主要是计算方法不同引起的。图5-3-6为图化的输出结果。
表5-3-1 ADI算法时间为500天时浓度场计算结果
表5-3-2 SIP算法时间为500天时浓度场计算结果
图5-3-6 SIP算法实例验证图
算例2
作者对三维SIP算法也进行了编程验证,作为变密度三维渗流溶质运移模型开发的基础条件,算例如下。利用SIP算法计算三维拉普拉斯方程:将边长为1的立方体剖分成20×20×20=8000的网格,Δx=Δy=Δz=0.05;网格的初始值均为0;立方体的下边、左边和后边的边界条件均为零通量边界,即U(0,Y,Z)=U(X,0,Z)=U(X,Y,0)=0;立方体的上边、前边和右边的边界条件分别为U(X,Y,1)=X×Y,U(1,Y,Z)=Y×Z,U(X,1,Z)=X×Z,(来自J.H.Ferziger和M.Peric的著作《Computational Methods for Fluid Dynamics》,2002)
方程离散格式如下:
利用SIP迭代求解结果如图5-3-7和图5-3-8所示。
图5-3-7 SIP算法验证3D切片图(从左到右从上到下z为1~21)
注:横坐标为1~21,间隔为1;纵坐标为0~1,间隔为0.1
图5-3-7中从左到右从上到下21张图分别表示z=1,21时的切片图。图5-3-8表示的是其三维立体图[顶角格点坐标为U(20,20,20)],两图能清楚地反映出算法的有效性,即在立方体的最上端边长x=1、y=1、z=1时数值最大,其余内部格点值按照离最上端远近距离均匀分布。
图5-3-8 SIP算法3D三维验证图
有关山东省水安全问题与适应对策:理论与实践的文章
Ccond=0.000 1时小流量和设计流量工况下离心泵扬程随有效空化余量变化的非定常数值计算结果与试验结果的比较如图2.6所示。表2.2NPSHc计算值与试验结果的比较在本书的离心泵空化流动数值模拟中,质量输运空化模型的凝结项经验系数取0.000 1进行数值计算。图2.7列出了小流量下NPSHa=1.0m和设计流量下NPSHa=1.6m时,3个不同凝结项经验系数所对应的叶轮内空泡体积率分布。......
2023-06-15
根据前面所介绍的串行嵌套多尺度实施方法,对一些简单算例进行分析,以验证算法的可行性及有效性,并测试两种不同方法的计算能力、计算精度及效率。表5.2不同方法计算所得等效刚度参数单位:GPa2.单轴拉伸宏观等效力学性能验证为验证嵌套多尺度方法用于材料损伤后的力学性能分析的可行性,对如图5.3所示的体元在单轴拉伸作用下损伤后的等效力学性能进行计算。图5.4考虑局部损伤演化的宏观等效应力与忽略损伤的比较......
2023-08-26
在平均化嵌套有限元多尺度计算过程中,宏观尺度应力并不是直接从宏观本构方程计算得到,需要通过细观尺度有限元计算所得细观应力进行平均化处理后得到。有关调用UMAT与ABAQUS内核语言Python并嵌入ABAQUS平台实现嵌套多尺度算法的更多阐述,有兴趣的读者可参考文献[24]。......
2023-08-26
图2.3朗格朗日网格与欧拉网格对比欧拉算法也有其不足,体现为单个循环计算时间长、材料边界不清晰、网格区域过大、冲击波耗散大、强度模拟不精确等。图2.4欧拉-拉格朗日耦合算法模型4.SPH算法SPH算法,即光滑粒子流体动力学数值算法,为固体材料大变形,尤其是存在破坏、断裂等极大变形的非线性动力学行为数值模拟提供了新的手段。......
2023-06-18
由前面的理论分析可知,在洞库储油压力为0.2MPa条件下,当洞库水封厚度为15m时,由于洞库中的气体压力大于洞库上方的水体水封能力,洞库中的气体将推动裂隙水运动向上方和两侧运移,为气体向岩体中迁移提供了的空间,并进入到上方的非饱和区,与大气相连通。当洞库水封厚度等于20m时,洞库中积聚的油气部分进入到岩体中,但由于水体厚度较大,油气压力不足以推动上方裂隙水发生大规模的迁移,因而被有效地封存于洞库上方岩体及洞库中。......
2023-06-28
计算这些序列的频率和时间平均方差形成特征向量,利用此特征向量数据进行了多方面的实验,验证其在人的行为识别方面的有效性。实验表明,当分段长度达到30帧以上时,就可获得很高的分类精度,且分段长度的变化对识别精度影响就会很小了。......
2023-06-16
回代入原通量密度式得其中Fe=,令aE=∣DeA∣,有Je=p-aEE 类似可得Jw=aWW-p 对南北方向类似地有Jn=p-aNN Js=aSS-p 将以上四项回代入原式得aPP=aEE+aWW+aNN+aSS+b 对界面上函数及其导数采取特定的构造格式,即确定了系数A和B的具体形式,就可以获得最终的离散表达式。......
2023-06-26
考虑水头、密度与浓度三者相互影响的可溶混咸淡水中溶质运移问题,需要由两个偏微分方程描述,一是考虑密度随浓度不断变化的液体流动,二是溶质运移。可用一般的描述溶质运移的对流-弥散方程描述海水入侵模型中盐分的运移。地下水溶质运移方程的离散格式。在溶质运移方程空间离散化时,为克服对流项占优时引起的过量或者数值弥散的情况,在对流扩散方程中的对流项乘以一个适当的权因子,称为上游加权法。......
2023-06-25
相关推荐