第六章变几何问题连续和移动边界.docx
《第六章变几何问题连续和移动边界.docx》由会员分享，可在线阅读，更多相关《第六章变几何问题连续和移动边界.docx（27页珍藏版）》请在冰点文库上搜索。
第六章变几何问题连续和移动边界
第六章变几何问题：
连续和移动边界
V.R.GUNDABALA,W.B.JZIMMERMAN,A.F.ROUTH
1DepartmentofChemicalandProcessEngineering,UniversityofSheffield,
NewcastleStreet,SheffieldS13JDUnitedKingdom
2DepartmentofChemicalEngineering,CambridgeUniversity,PembrokeStreet,
Cambridge,CB23RA
E-mail:
求解过程中，由于几何模型发生变化而导致区域网格改变时会出现几何结构的连续性问题。
本章中，我们举两个例子——流道中由于不同尺寸孔板引起额外压力损失问题是稳态几何连续性的一个例子。
从概念来看，这个问题与第5章贝纳尔问题中通过Rayleigh数进行的变量连续性没什么不同。
第二个例子是液体中悬浮的乳胶颗粒的干膜。
在这个问题中，考虑了两个移动前缘；一个通过坐标变换，另一个通过光滑界面模型。
该技术对于之前模拟该问题采用的移动弱项方法有所改进。
1.引言
1.1几何连续性
我们已经有了一些参数连续性的例子——通过参数在一定范围内一系列的较小变化，以相邻参数值的解作为新参数求解初值。
只要参数不通过歧点且参数步长足够小，即可保证从旧解到新解的平滑过渡。
即使存在歧点，旧的分支仍然是一个可行解，参见我们在第5章中对贝纳尔对流问题的讨论。
几何连续性与参数连续性在一个重要方面存在性质上面的不同。
在几何连续性中，区域几何结构的改变使得需要对网格重新划分。
我们应该小心的区分几何连续性的改变。
例如：
在管流中，众所周知流动可由雷诺数表征：
（1）
该无量纲数包含了流体密度ρ，入口速度U，直径D，和粘度μ，可用于描述流体的动力相似性。
因此，管道直径改变引起充分发展流并不属于几何连续性的范围，而更接近于常规的参数改变问题。
与上一章类似，模拟的例子给出了稳态和瞬态的模型，本章我们将给出稳态几何连续性与瞬态几何连续性的例子。
对于前者，不同的模型所用的几何结构不同，因此网格会略有不同。
因此，一个参数到下一个参数的解是不兼容的（不同的网格自由度）。
如果旧的解作为新几何结构的初始值，则需要将旧的解在保持一致性的情况下映射到新的区域。
在这里的例子中，稳态模型PDE方程组是线性的，因此解可以直接由1FEM步获得。
因此，将旧的解映射到新几何结构没有附加值。
对于瞬态问题，随着移动的前缘，涉及到区域减小的问题。
该区域在每一时间步后均发生改变，因此将上一时间步的解映射到新的区域是必要的。
这导致在每一时间步后必须重新划分网格。
对于自由边界问题这是关键的一步。
例如膜流动和喷气流，边界的位置和速度场的解密切相关，边界应位于能够满足压力平衡的位置。
二维不可压、层流纳维-斯托克斯方程能以一系列标准方法求解（有限差分，有限元，谱元，格子-波尔兹曼以及多重网格技术），对于固体边界条件以及复杂几何形状已能够由商用标准模拟引擎进行计算。
标准计算流体动力学包有两个标准引擎：
（1）能够满足复杂几何形状的网格生成器，
（2）能够求解包含压力项（如纳维-斯托克斯方程中对于连续性方程的拉格朗日乘子）通用传递方程组的偏微分方程求解器。
这两步通常分别进行，先生成网格，后进行求解，FMELAB软件在这方面并无不同。
1.2自由边界问题和任意拉格朗日-欧拉变换
计算流体动力学对于自由边界问题处理的并非很好。
可以考虑耦合流动解和网格生成的迭代算法，但是自动化的标准计算包很难实现。
Ruschak[1]描述了使用网格自适应方法实现压力边界条件的标准方法。
Goodwin和Schowalter[2]成功应用于对毛细管粘度喷射的处理之中对网格位置与流体方程和边界条件的同时求解，采用有限元方法对模型方程进行离散并采用牛顿迭代法求解。
FIDAP，通常用于处理自由边界流体，使用迭代流体求解/椭圆网格重生成方法，而不是同时牛顿迭代。
在缩小了的区域中的瞬态模型，我们采用坐标变换和平滑界面的方法以改变整个时间域上的几何域及两个移动边界条件。
通常对于考虑多相流的模拟，由于平滑界面模型允许拓扑结构的改变，因此对于移动边界模型较为适用。
本书第8章以水平集合方法求解液滴结合的问题。
原则上，COMSOLMultiphysics软件也能够采用ALE模式（任意拉格朗日-欧拉变换）求解后者问题。
ALE模式可通过COMSOLMultiphysics中集成的方程来实现，随着移动网格位置剩余而增加。
在所有的应用模式下，包含这些附加项的标准化方法均为“激活”网格。
例如，允许给定自由边界是对上一版本FEMLAB的改进，但涉及到自由边界计算的问题较少，只有少数如考虑到表面张力引起的流动等问题。
实际上，对FEMLAB软件包的这一改动并非主要针对这个问题。
考虑机械动作（压气机与透平）与传递现象联系的需要是这个改动的一个目的，另一个目的是为了考虑结构力学与传递现象的联系。
我本希望用一个简单的模型来说明ALE移动边界条件耦合的方法，但这里建议读者阅读COMSOLMultiphysics模型库中蠕动泵和挑梁模型。
2.稳态几何连续性：
孔板引起的流道压降
在本节中，我们将考虑两个需要考虑几何连续性的模型，分别是孔板和粘性流管道中的片状物。
这两个模型是相互联系的，事实上仅存在微小差别。
图1为孔板几何结构，我们可估计当板阻塞率为ε时，保持恒定体积流率所需的额外压降。
图中绘制了当板阻塞率为40%时所生成的典型网格。
由于有必要分辨孔板周围的微小特征，孔板周围的网格必须保持足够致密，网格加密成为一个重要问题。
我们将从几何参数改变来外推所需的附加压降，而不是直接求解高阻塞率问题。
图1.粘性流管道中孔板问题生成的网格。
用于表示阻塞百分数的参数为ε＝0.4
尽管有可能考虑孔板周围任意雷诺数的计算，层流的主要影响与人为消去的雷诺数-斯托克斯方程相似。
该现象的根由在于大多数的耗散发生在孔板缝隙之间，在该位置，流体加速，然而小的缝隙导致了较强的粘性摩擦，可近似根据斯托克斯方程模拟动量传递：
（2）
其中，μ为流体粘度，ρ为密度，所有的均代表它们在流体力学中的常用意义。
方程
（2）有量纲，准稳态，无惯性项。
由于该方程是线性的，对此方程有过详尽的分析。
H.Ockendon和J.K.Ockendon[3]等人的文章可提供较好参考。
Homsy等[4]提供了一系列很好的对小雷诺数粘流消失可视化分析。
我对孔板问题的注意来自于Dugdale教授，他通过约束孔板自身能量耗散最小获得了锋利孔板附近的解，求解过程中忽略了容器其他边界上的能量耗散。
他的观点是因为孔板非常小，所有的流体强制通过孔板，几乎所有的能量必须通过孔板耗散，对于缝隙长度为a的三维孔板给出量纲分析，能量耗散速率E须满足：
（3）
其中，Q为体积流率，W为功流。
c为未知比例常数，该常数数值由Dugdale根据吸附能极值从理论上进行了计算，也能够通过测量Q和压降从实验上获得。
在二维体系中，通过相似量纲分析可获得单位长度耗散损失
和截面面积流率
，从而可得：
（4）
Dugdale报道了通过糖浆实验确定常数c范围在3.17至3.30之间，他的理论计算结果为3.30。
Bond[6]给出了孔板问题与管长为2ka的哈根-伯叔叶管流问题的相似性，其中a为孔板半径，压降等值为k＝0.631，表明c=3.21.
我们曾对管道中紧合粒子的拖曳较感兴趣。
对于与（3）（4）相同的理论基础，管道中的颗粒由狭缝宽度决定。
Zimmerman对于薄片[7]（宽边方向的运动）以及球体[8]在圆柱管中的沉积做了相关研究，表明当颗粒半径较大时曳力增大较快（狭缝宽度为a）。
对小颗粒半径（1-a）采用摄动法并将级数展开相加，有可能确定当颗粒接近管壁时问题的奇异性。
通过对薄片宽边方向运动问题的量纲分析，（4）可用于考虑二维或非中心对称狭缝的二阶奇异特性，O（a-2）。
由于缝隙宽度随着与球体中心相对的角度不断发生变化，球体问题不适合采用量纲分析方法。
Bungay和Brenner[9]计算了球体曳力的奇异性阶数为O（a-5/2）。
Harlen[10]采用有限元方法,发现该问题较难收敛，也揭示了解析大尺度范围差分问题在数值上存在的困难。
即使对于线性模型，当小长度尺度决定了流体的动力学特性时，求解仍非常困难。
我猜测较小狭缝宽度下的动力学问题可以通过对较大狭缝宽度问题解的外推来获得。
本章中，我们首先计算阻塞率为ε的孔板存在引起的流道附加压损。
缝隙半径与阻塞率相关a=1-ε.这个问题与颗粒沉积问题在概念上的不同较小。
例如Shail和Norton[11]计算了薄片宽度方向在圆柱流道中的运动，也计算了与之相对应的问题－保持薄片在流道中静止。
这些问题的定量关系由于方程
（2）的线性特性也成线性相关，有望通过一个问题来反映另一个问题。
表1孔板插入流道的二维模型。
文件名orifice05.m
模型导航栏
选择二维空间
选择COMSOLMultiphysics|FluidDynamics|IncompressibleNavier-Stokes（modens）
完成
Draw菜单
定义物体
定义直线，类型：
多段线
坐标x：
00222.052.05550
y：
0110.950.951100
完成
Draw菜单|Coerceto|Solid
以上生成复合结构CO1
Options菜单|
Constants
常数名称：
rho0，表达式：
0
常数名称：
mu0，表达式：
1
常数名称：
Umean，表达式：
1，完成
Physics菜单|
Boundarysettings
选择边界1，设定进/出口边界条件。
u0=6*Umean*s*（1-s）;v0=0
选择边界8，设定常规流动/压力边界条件，p=0（默认）
接受其它边界默认无滑移边界条件。
完成
Physics菜单|
Subdomainsettings
选择子域1
设定ρ=rho0;η=mu0;Fx=0;Fy=0
选择init选项卡，设定u（t0）=Umean,v（t0）=0;p（t0）=0
完成
Mesh
点击工具栏上三角形按钮生成默认网格
加密网格两次
Solve菜单
求解管理器，选择稳态非线性求解器
Advanced选项卡：
设定变量缩放为空
点击工具栏上求解按钮（＝）
保存为orifice05.m文件
2.1孔板插入2-D流道的模型
表1为流道中孔板问题模型的设置说明。
在绘图模式中，通过矩形相减获得复合几何对象，并设置孔板阻塞率ε。
较繁琐的一个方法是建立多段线并将其转换为可编辑的几何对象。
我们的目的是建立能够自动改变几何体参数的MATLAB脚本，这是因为在COMSOLMultiphysicsGUI中实现自动改变几何体参数是很困难的。
入口边界条件是二维流道中充分发展的哈根－伯庶叶流，Umean为表征入口条件的单参数。
弧长参数s用于表征竖直方向的依从关系。
s在每个边界均可用，取值范围在0到总弧长范围之内。
由于抛物线速度场与边界信号方向无关，因此没有必要考虑信号方向。
用标准网格参数，点击重新划分网格按钮两次，可得到如图1所示网格。
注意输出给出了压力基准，因此我们可以得到相应的压力分布。
输出环境可有多种选择，标记为正常流体/压力，在流道末端重生成哈根-伯庶叶流体的压降。
应用哪一个边界条件是一个模拟方面的问题，涉及到将所模拟问题进行理想化及估算。
在这种情况下，管道上游有一段流体发展段，出口条件为“straightout”，表明出口与入口也类似。
在求解器参数下，“scalingofvariables”设置为无选择，适用于平均流速已经确定的问题。
此外，我们选定密度为零，使得这个问题成为线性问题以改进模型的收敛性。
线性问题已有成熟的收敛性判断方法－单矩阵逆变换。
无论如何，我们采用稳态非线性求解器－因为该求解器是不可压缩纳维－斯托克斯方程的默认求解器，并非一个可选的求解器。
经过两次迭代，误差可以较精确的满足COMSOLMultiphysics的要求，初始误差10-12是一个较好的迭代起点。
图2给出了当ε＝0.05时速度场的箭头图。
图2ε＝0.05凹槽的速度箭头图和等压线图
图3ε＝0.05凹槽的中心线压力分布图
图3为中心线方向的压力分布，由后处理菜单中的横截面参数绘图生成。
因边界为单位长度，在边界1上的压力积分给出了边界上的平均压力。
该数值在消息框中显示为60.462。
计算哈根－伯庶叶流动压降△p=60。
充分发展u速度分布为：
（5）
代入方程
（2）求得常压力梯度为-12Umean。
在下游五个单位长度距离上，可计算
以保证出口压力为p=0。
因此，附加压降为0.462（由于粘性力与速度的比例关系而无量纲化）。
练习6.1
细分网格并计算附加压降。
使用工具栏中的标准细化按钮，使用旧的解作为计算初值。
对网格的统一性和附加压降的变化作一评论，值得再一次细化网格吗？
现在回到绘图模式，双击凹槽底部的竖直部分。
编辑孔板阻塞率为40%,保持孔板宽度仍为0.05.求解。
图4为速度适量的箭头图及压力的等值图。
很明显，速度分布必须在角落转弯，由此导致更多的流体中断，也因此造成了更多的能量耗散。
边界积分给出了压力损失为△p=84.85。
注意沿出口流体x-速度的边界积分为1，是Umean的数值。
图4给出了等压线，明显看出压力在孔板附近的迅速耗散。
由于需要强制保证流体通过孔板，因此在孔板附近的上游位置出现最大压力。
我们如何实现几何连续性？
在这种情况下，又有问题是线性的，所有的模型都利用相邻几何参数（阻塞率）获得较好的收敛性能。
尽管如此，还需要认真研究网格细化问题以保证分辨率。
首先，输出一个模型m文件。
然后编辑该文件实现几何结构参数化。
本章用到的MATLABm文件的第一部分如下：
%%%%%%%%%%%%%%%%%orifice.m%%%%%%%%%%%%%%%%%%%%%%%%%%%%
slot=[0.05:
0.05:
0.75];
%WZ:
Setupstorage
output=zeros（length（slot）,5）;
%WZ:
NowlooparoundthewholeFEMLABmodelm-filewithj
forj=1:
length（slot）eps=slot（j）;
%Constants
={’rho0’,’0’,’mu0’,’1’,’Umean’,’1’};
%Geometry
carr={curve2（[0,0],[0,1]）,curve2（[0,2],[1,1]）,
curve2（[2,2],[1,1-eps]）,...
curve2（[2,2.05],[1-eps,1-eps]）,curve2（[2.05,2.05],[1-eps,1]）,...
curve2（[2.05,5],[1,1]）,curve2（[5,5],[1,0]）,curve2（[5,0],[0,0]）};
g1=geomcoerce（’curve’,carr）;g2=geomcoerce（’solid’,{g1}）;clears
s.objs={g2};s.name={’CO1’};s.tags={’g2’};
=struct（’s’,s）;=geomcsg（fem）;
%Initializemesh
=meshinit（fem,’report’,’off’）;
%Refinemesh
=meshrefine（fem,’mcase’,0,’rmethod’,’regular’）;
=meshrefine（fem,’mcase’,0,’rmethod’,’regular’）;
=meshrefine（fem,’mcase’,0,’rmethod’,’regular’）;
%Applicationmode1
clearappl.class=’FlNavierStokes’;er={4,2};
er={2,1};nsuffix=’_ns’;clearprop
sis=’static’;=prop;clearbndbnd.u0=
{0,’6*Umean*s*（1-s）’,0};={’noslip’,’uv’,’strout’};
=[2,1,1,1,1,1,1,3];=bnd;clearequ=
{{’Umean’;0;0}};=’rho0’;=’mu0’;er=
{{1;1;2}};er={{1;1;2}};=[1];=equ;
{1}=appl;={’ref’};r=1;=
’SI’;
%Multiphysics
fem=multiphysics（fem）;
%Extendmesh
=meshextend（fem）;
%Solveproblem
=femnlin（fem,’solcomp’,{’u’,’p’,’v’},...
’outcomp’,{’u’,’p’,’v’},’uscale’,’none’,’report’,’off’）;
fem0=fem;
%Solveproblem
fem=adaption（fem,’init’,,’solcomp’,{’u’,’p’,’v’},...
’outcomp’,{’u’,’p’,’v’},’ntol’,1e-006,’nonlin’,’on’,...
’solver’,’stationary’,’l2scale’,[1],’l2staborder’,[2],...
’eigselect’,[1],’maxt’,,’ngen’,2,’resorder’,[0],...
’rmethod’,’longest’,’tppar’,1.7,’uscale’,’none’,...
’geomnum’,1,’report’,’off’）;
%Integrate
I1=postint（fem,’u’,’dl’,[1],’edim’,1）;
I2=postint（fem,’K_x_ns’,’dl’,[4,5,6],’edim’,1）;
I3=postint（fem,’p’,’dl’,[1],’edim’,1）;
%WZ:
writeoutoutputtothearrayoutput
output（j,1）=slot（j）;output（j,2）=I1;output（j,3）=I2;output（j,4）=I3;
output（j,5）=length（.p）;
%WZ:
endourj-loop
end
%WZ:
recordresultstofile
savefem
dlmwrite（’’,output,’,’）;
图4.速度场箭头矢量图及等压线（ε=0.05）
该m文件脚本可能是本书中最长的，因此保留了部分注释。
值得注意的第一点是几何结构是由脚本通过调整阻塞率生成的。
通过建立练习中ε=0.4的几何对象。
大多数的程序语句由COMSOLMultiphysicsGUI生成，选择“另存为m文件”进行保存。
自行添加的主要是控制程序以及结果保存部分的代码。
注意这里采用了自适应求解器，可以由GUI来实现（将自适应求解器参数面板上的对话框打勾即可）。
对整个脚本的求解在一个常规的linux工作站下大约需要10CPU分的时间。
表2中对自由度的估计通过.u长度来给出。
表2孔板模型的主要输出结果。
平均出口x方向速度，
积分粘性应力，自由度数目随孔板阻塞率的变化。
ε
出口速度
粘性剪切力
△p
自由度
0.05
1
-0.78522
60.463
90511
0.1
1
-1.0336
61.551
94235
0.15
1
-1.2105
63.248
100340
0.2
1
-1.4211
65.605
99151
0.25
1
-1.5548
68.729
103260
0.3
1
-1.8241
72.788
100920
0.35
1
-1.9963
78.043
96579
0.4
1
-2.3204
84.864
106450
0.45
1
-2.6338
93.817
103000
0.5
1
-3.1446
105.78
104110
0.55
1
-3.6672
122.21
103880
0.6
1
-4.383
145.46
105340
0.65
1
-5.2545
179.9
99807
0.7
1
-6.933
233.9
99231
0.75
1
-9.2868
325.69
92935
表2由几何连续性研究的结果组成。
明显的，附加压损随着阻塞因子的增加而快速增加。
粘性应力从ε=0.05到ε=0.75也存在快速的增加。
（6）
尽管（6）式在从ε=0.05到ε=0.75的范围内较好的拟合了所得数据，但如果在一个较小的范围ε∈[0,0.5]，对1－ε的劳伦特级数可给出对数据更好的拟合。
（7）
图5和图6绘制了根据方程（5）和（7）得到的曲线以及计算所得数据点。
由（7）式得出的附加压损为△p=371，由（6）则为214,模型给出309。
（7）式的关键特点在于如果考虑到式中分子中系数项几乎相等，该结构与量纲分析的结果基本一致。
图5△p随ε的变化（孔板厚度0.05）,曲线为ε在[0,0.5]之间的三次曲线拟合
图6由1－ε逆幂的劳伦特级数进行拟合
图7Re=1500,ε=0.4。
注意背流面粘滞流区域。
流体的惯性导致流线跨过该粘滞区域。
在粘性流中，流体“转过拐角”并没有分离
练习6.2
Dugdale的孔板是锋利的。
我们的孔板厚度则为整个流道厚度的5%。
试着将孔板变得更薄：
4％，3％，2%。
对附加压损会有怎样的影响呢？
根据文献[12]，当缝隙宽度较小时，孔板厚度对曳力有一定的影响。
如果缝隙是扁平的，则Dugdle的量纲分析是正确的，如方程（3），但是如果颗粒是有限弧度的，则Bungay和Brenner的结果O（a-5/2）更为合理。
练习6.3
（a）改变孔板的上边界条件为对称边界条件。
这个模型是二维板状物的粘性绕流问题。
尝试几何连续性。
（b）改变m文件，不令上一次的几何结构作为下一次计算的初始条件。
计算快一些了吗？
（c）这个例子并不是真正存在的物理现象，试着增加流函数-涡fangcheng用于考虑浮力对流的影响。
平板问题由Kim[13]进行了研究，在缝隙宽度较小时采用了摄动法，将展开的级数项加和得到曳力的奇异特性。
练习6.4
这些例子没有精确的采用几何连续性，这是因为由于Re=0几何形状改变成为线性问题，使得由于惯性力造成的非线性项随之消失。
假设Re=1500.在s中给出常数，设置rho0=1500。
图7中，注意背流面粘滞流区域。
流体惯性使得流体转过拐角，未造成分离。
当Re增大时，该流型结构与线性流有显著不同，求解较为困难。
有两个有效的策略：
对于Re的参数连续性或者几何结构角度的连续性。
后者方法涉及到使用asseminit命令将旧的解内插到新的网格上：
＝assemini
升级会员 