变何问题连续和移动边界

1、第六章 变几何问题：连续和移动边界V.R.GUNDABALA, W.B.J ZIMMERMAN, A.F.ROUTH1Department of Chemical and Process Engineering, University of Sheffield,Newcastle Street, Sheffield S1 3JD United Kingdom2Department of Chemical Engineering, Cambridge University, Pembroke Street,Cambridge, CB2 3RAE-mail: w.zimmermanshef.ac.

2、uk求解过程中，由于几何模型发生变化而导致区域网格改变时会出现几何结构 的连续性问题。本章中，我们举两个例子一一流道中由于不同尺寸孔板引起额外 压力损失问题是稳态几何连续性的一个例子。从概念来看，这个问题与第5章贝纳尔问题中通过Rayleigh数进行的变量连续性没什么不同。第二个例子是液体 中悬浮的乳胶颗粒的干膜。在这个问题中，考虑了两个移动前缘；一个通过坐标 变换，另一个通过光滑界面模型。该技术对于之前模拟该问题采用的移动弱项方 法有所改进。1 .引言1.1 几何连续性我们已经有了一些参数连续性的例子一一通过参数在一定范围内一系列的 较小变化，以相邻参数值的解作为新参数求解初值。只要参数不通

3、过歧点且参数 步长足够小，即可保证从旧解到新解的平滑过渡。 即使存在歧点，旧的分支仍然 是一个可行解，参见我们在第 5章中对贝纳尔对流问题的讨论。几何连续性与参数连续性在一个重要方面存在性质上面的不同。 在几何连续 性中，区域几何结构的改变使得需要对网格重新划分。 我们应该小心的区分几何 连续性的改变。例如：在管流中，众所周知流动可由雷诺数表征：UDRe 该无量纲数包含了流体密度 p,入口速度U,直径D,和粘度也可用于描述流 体的动力相似性。因此，管道直径改变引起充分发展流并不属于几何连续性的范 围，而更接近于常规的参数改变问题。与上一章类似，模拟的例子给出了稳态和瞬态的模型， 本章我们将给出

4、稳态 几何连续性与瞬态几何连续性的例子。 对于前者，不同的模型所用的几何结构不 同，因此网格会略有不同。因此，一个参数到下一个参数的解是不兼容的（不同 的网格自由度）。如果旧的解作为新几何结构的初始值，则需要将旧的解在保持 一致性的情况下映射到新的区域。在这里的例子中，稳态模型PDE方程组是线性的，因此解可以直接由1FEM步获得。因此，将旧的解映射到新几何结构没有 附加值。对于瞬态问题，随着移动的前缘，涉及到区域减小的问题。该区域在每 一时间步后均发生改变，因此将上一时间步的解映射到新的区域是必要的。这导 致在每一时间步后必须重新划分网格。 对于自由边界问题这是关键的一步。 例如 膜流动和喷气

5、流，边界的位置和速度场的解密切相关，边界应位于能够满足压力 平衡的位置。二维不可压、层流纳维-斯托克斯方程能以一系列标准方法求解（有限差分，有限元，谱元，格子-波尔兹曼以及多重网格技术），对于固体边界条件以及复杂 几何形状已能够由商用标准模拟引擎进行计算。标准计算流体动力学包有两个标 准引擎：（1）能够满足复杂几何形状的网格生成器，（2）能够求解包含压力项（如 纳维-斯托克斯方程中对于连续性方程的拉格朗日乘子）通用传递方程组的偏微 分方程求解器。这两步通常分别进行，先生成网格，后进行求解，FMELAB软件在这方面并无不同。1.2 自由边界问题和任意拉格朗日-欧拉变换计算流体动力学对于自由边界问

6、题处理的并非很好。可以考虑耦合流动解和网格生成的迭代算法，但是自动化的标准计算包很难实现。Ruschak1描述了使用网格自适应方法实现压力边界条件的标准方法。Goodwin和Schowalter2成功应用于对毛细管粘度喷射的处理之中对网格位置与流体方程和边界条件的 同时求解，采用有限元方法对模型方程进行离散并采用牛顿迭代法求解。FIDAP,通常用于处理自由边界流体，使用迭代流体求解 /椭圆网格重生成方法，而不是 同时牛顿迭代。在缩小了的区域中的瞬态模型，我们采用坐标变换和平滑界面的 方法以改变整个时间域上的几何域及两个移动边界条件。通常对于考虑多相流的模拟，由于平滑界面模型允许拓扑结构的改变，

7、因此对于移动边界模型较为适用。 本书第8章以水平集合方法求解液滴结合的问题。原则上，COMSOL Multiphysics软件也能够采用ALE模式（任意拉格朗日- 欧拉变换）求解后者问题。 ALE模式可通过COMSOL Multiphysics中集成的方 程来实现，随着移动网格位置剩余而增加。在所有的应用模式下，包含这些附加 项的标准化方法均为“激活”网格。例如，允许给定自由边界是对上一版本 FEMLAB的改进，但涉及到自由边界计算的问题较少，只有少数如考虑到表面 张力引起的流动等问题。实际上，对 FEMLAB软件包的这一改动并非主要针对 这个问题。考虑机械动作（压气机与透平）与传递现象联系的

8、需要是这个改动的 一个目的，另一个目的是为了考虑结构力学与传递现象的联系。我本希望用一个简单的模型来说明ALE移动边界条件耦合的方法，但这里建议读者阅读COMSOL Multiphysics模型库中蠕动泵和挑梁模型。2 .稳态几何连续性：孔板引起的流道压降在本节中，我们将考虑两个需要考虑几何连续性的模型，分别是孔板和粘性流管道中的片状物。这两个模型是相互联系的，事实上仅存在微小差别。图 1 为孔板几何结构，我们可估计当板阻塞率为e时，保持恒定体积流率所需的额外 压降。图中绘制了当板阻塞率为40%时所生成的典型网格。由于有必要分辨孔板 周围的微小特征，孔板周围的网格必须保持足够致密，网格加密成为

9、一个重要问 题。我们将从几何参数改变来外推所需的附加压降，而不是直接求解高阻塞率问题。图1,粘性流管道中孔板问题生成的网格。用于表示阻塞百分数的参数为6 0.4尽管有可能考虑孔板周围任意雷诺数的计算， 层流的主要影响与人为消去的 雷诺数-斯托克斯方程相似。该现象的根由在于大多数的耗散发生在孔板缝隙之 间，在该位置，流体加速，然而小的缝隙导致了较强的粘性摩擦，可近似根据斯 托克斯方程模拟动量传递：(2)其中，以为流体粘度，p为密度，所有的均代表它们在流体力学中的常用意义。方程(2)有量纲，准稳态，无惯性项。由于该方程是线性的，对此方程有过详尽的分析。H.Ockendon和J.K.Ockendon

10、3等人的文章可提供较好参考。Homsy等4提供了一系列很好的对小雷诺数粘流消失可视化分析。我对孔板问题的注意 来自于Dugdale教授，他通过约束孔板自身能量耗散最小获得了锋利孔板附近的 解，求解过程中忽略了容器其他边界上的能量耗散。 他的观点是因为孔板非常小， 所有的流体强制通过孔板，几乎所有的能量必须通过孔板耗散，对于缝隙长度为 a的三维孔板给出量纲分析，能量耗散速率 E须满足：Q2-Ec3- W Q p3 a其中，Q为体积流率，W为功流。c为未知比例常数，该常数数值由 Dugdale根 据吸附能极值从理论上进行了计算，也能够通过测量Q和压降从实验上获得。在二维体系中，通过相似量纲分析可获

11、得单位长度耗散损失E和截面面积流率Q,从而可得:QE c W Q p（4）aDugdale报道了通过糖浆实验确定常数 c范围在3.17至3.30之间，他的理论 计算结果为3.30。Bond6给出了孔板问题与管长为2ka的哈根-伯叔叶管流问题 的相似性，其中a为孔板半径，压降等值为k= 0.631,表明c=3.21.我们曾对管道中紧合粒子的拖曳较感兴趣。对于与（3）（4）相同的理论基础，管道中的颗粒由狭缝宽度决定。Zimmerman对于薄片7（宽边方向的运动）以 及球体8在圆柱管中的沉积做了相关研究，表明当颗粒半径较大时曳力增大较 快（狭缝宽度为a）。对小颗粒半径（1-a）采用摄动法并将级数展开

12、相加，有可能 确定当颗粒接近管壁时问题的奇异性。通过对薄片宽边方向运动问题的量纲分 析，（4）可用于考虑二维或非中心对称狭缝的二阶奇异特性，O（a-2）。由于缝隙宽度随着与球体中心相对的角度不断发生变化，球体问题不适合采用量纲分析方 法。Bungay和Brenner9计算了球体曳力的奇异性阶数为 O（a-5/2）。Harlen10采 用有限元方法，发现该问题较难收敛，也揭示了解析大尺度范围差分问题在数值 上存在的困难。即使对于线性模型，当小长度尺度决定了流体的动力学特性时， 求解仍非常困难。我猜测较小狭缝宽度下的动力学问题可以通过对较大狭缝宽度 问题解的外推来获得。本章中，我们首先计算阻塞率为

13、e的孔板存在引起的流道附加压损。缝隙半 径与阻塞率相关 a=1- &这个问题与颗粒沉积问题在概念上的不同较小。例如 Shail和Norton11计算了薄片宽度方向在圆柱流道中的运动，也计算了与之相对应的问题一保持薄片在流道中静止。这些问题的定量关系由于方程（2）的线性特性也成线性相关，有望通过一个问题来反映另一个问题。表1孔板插入流道的二维模型。文件名orfice05.m模型导航栏选择二维空间选择 COMSOL Multiphysics| Fluid Dynamics|Incompressible Navier-Stokes (mode ns)完成Draw菜单定义物体定义直线，类型：多段线坐标

14、 x: 0 0 2 2 2.05 2.05 5 5 0V： 0 1 1 0.95 0.95 1 1 0 0完成Draw 菜单 |Coerce to| Solid以上生成复合结构CO1Options 菜单 | Constants常数名称：rho0,表式：0常数名称： mu0,表式：1常数名称：Umean,表辽式：1,元成Physics 菜单 |Boundary settings选择边界1,设定进/出口边界条件。u0=6*Umean*s*（1-s）； V0=0选择边界8,设定常规流动/压力边界条件，p=0 （默认）接受其它边界默认无滑移边界条件。完成Physics 菜单 |Subdomain se

15、ttings选择子域1设定 f=rho0; r=mu0; Fx=0; Fy=0选择 init 选项卡，设定 u(t0)=Umean, v(t0)=0； p(t0)=0 完成Mesh点击工具栏上三角形按钮生成默认网格 加密网格两次Solve菜单求解管理器，选择稳态非线性求解器Advanced选项卡：设定义量缩放为空 点击工具栏上求解按钮（=）保存为orfice05.m文件2,1孔板插入2-D流道的模型表1为流道中孔板问题模型的设置说明。在绘图模式中，通过矩形相减获 得复合几何对象，并设置孔板阻塞率e o较繁琐的一个方法是建立多段线并将其 转换为可编辑的几何对象。我们的目的是建立能够自动改变几何体

16、参数的 MATLAB脚本，这是因为在 COMSOL Multiphysics GUI中实现自动改变几何体 参数是很困难的。入口边界条件是二维流道中充分发展的哈根伯庶叶流， Umean为表征入口 条件的单参数。弧长参数S用于表征竖直方向的依从关系。S在每个边界均可用， 取值范围在0到总弧长范围之内。由于抛物线速度场与边界信号方向无关， 因此 没有必要考虑信号方向。用标准网格参数，点击重新划分网格按钮两次，可得到如图1所示网格。注 意输出给出了压力基准，因此我们可以得到相应的压力分布。 输出环境可有多种 选择，标记为正常流体/压力，在流道末端重生成哈根-伯庶叶流体的压降。应用 哪一个边界条件是一个

17、模拟方面的问题，涉及到将所模拟问题进行理想化及估 算。在这种情况下，管道上游有一段流体发展段，出口条件为“ straight out”，表 明出口与入口也类似。在求解器参数下，“scaling of variables”设置为无选择，适用于平均流速已 经确定的问题。此外，我们选定密度为零，使得这个问题成为线性问题以改进模 型的收敛性。线性问题已有成熟的收敛性判断方法-单矩阵逆变换。无论如何， 我们采用稳态非线性求解器-因为该求解器是不可压缩纳维-斯托克斯方程的 默认求解器，并非一个可选的求解器。经过两次迭代，误差可以较精确的满足 COMSOL Multiphysics的要求，初始误差10-12

18、是一个较好的迭代起点。图 2给 出了当 0.05时速度场的箭头图。diE.i! p期MI 弱 Anew flrfd |15三rW岁.三 tttmwtt.一二 rtTIWW11TTST. 三mwrrt 二 rtWtrm rTfrnwtrTT m rnWW一 三 fwm, 三 mrrtm / ftnftrtn hwfTTrd第k更-C.3I7 .晨八W3j三,二-3.54455工用可1械K但K)5=i7即UHL5UMh;U12图2 e =0.05凹槽的速度箭头图和等压线图图3 e =0.05凹槽的中心线压力分布图图3为中心线方向的压力分布，由后处理菜单中的横截面参数绘图生成。 因 边界为单位长度，

19、在边界1上的压力积分给出了边界上的平均压力。 该数值在消 息框中显示为60.462。计算哈根伯庶叶流动压降 p=60。充分发展u速度分布 为：U 6Umeany(1 Y)代入方程(2)求得常压力梯度为-12Umean。在下游五个单位长度距离上，可计算Pinlet 60Umean以保证出口压力为P=0。因此，附加压降为0.462 (由于粘性力与速度的比例关系而无量纲化)。练习6.1细分网格并计算附加压降。使用工具栏中的标准细化按钮，使用旧的解作为 计算初值。对网格的统一性和附加压降的变化作一评论，值得再一次细化网格 吗？现在回到绘图模式，双击凹槽底部的竖直部分。编辑孔板阻塞率为40%,保持孔板宽

20、度仍为0.05.求解。图4为速度适量的箭头图及压力的等值图。很明显， 速度分布必须在角落转弯，由此导致更多的流体中断，也因此造成了更多的能量 耗散。边界积分给出了压力损失为 p=84.85。注意沿出口流体x-速度的边界积分 为1,是Umean的数值。图4给出了等压线，明显看出压力在孔板附近的迅速耗 散。由于需要强制保证流体通过孔板，因此在孔板附近的上游位置出现最大压力。我们如何实现几何连续性？在这种情况下，又有问题是线性的，所有的模型都利用相邻几何参数(阻塞率)获得较好的收敛性能。尽管如此，还需要认真研 究网格细化问题以保证分辨率。首先，输出一个模型 m文件。然后编辑该文件 实现几何结构参数化

21、。本章用到的MATLAB m文件的第一部分如下：%orifice.m%slot=0.05:0.05:0.75;% WZ: Set up storageoutput=zeros(length(slot),5);% WZ: Now loop around the whole FEMLAB model m-file with j for j=1:length(slot) eps=slot(j);% Constantsfem.const = rho0 , 0 , mu0 , 1 , Umean , 1 ;% Geometrycarr=curve2(0,0,0,1), curve2(0,2,1,1),c

22、urve2(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 ,g2i=geomcoerce( solid ,g1); clear sfem.draw=struct(s.objs=g2; = CO1 ; s.tags= g2 ;s ,s); fem.geom=geomcsg(fem);reportoff )；% Initialize meshfe

23、m.mesh=meshinit(fem,fem.mesh=meshrefine(fem,mcase,0,rmethod);, regularfem.mesh=meshrefine(fem,mcase,0,rmethod, regularfem.mesh=meshrefine(fem,mcase,0,rmethod, regular% Refine mesh）；）；% Application mode 1clear appl appl.mode.class = FlNavierStokes ; appl.gporder = 4,2;appl.cporder = 2,1; appl.assigns

24、uffix =lear prdp _ns ; cprop.analysis= static ; p = prop; clear bnd bnd.u0 = 0, 6*Umean*s*(-s) ,0; bnd.type = noslip , uv , strout ;bnd.ind = 2,1,1,1,1,1,1,3; appl.bnd = bnd; clear equ equ.init = Umean ;0;0; equ.rho = rho0 ; equ.eta =porder =mu0 ; equ.c1;1;2; equ.gporder = 1;1;2; equ.ind = 1

25、; appl.equ = equ;fem.appl1 = appl; fem.frame = ref ; fem.border = 1; fem.units =SI ；% Multiphysics fem=multiphysics(fem);% Extend meshfem.xmesh=meshextend(fem);% Solve problemv ,none , reportoff );fem.sol =femnlin(fem, solcomp , u , p outcomp , u , p , v , uscalefem0=fem;fem=adaption(fem,% Solve pro

26、blem,1e oninit ,fem0.sol, solcomp , u , p , v , .solver , stationary eigselect ,1, maxtrmethod, longest ,geomnum ,1, report outcomp , u , p , -00v6, , nonltionll2scale ,1, l2staborder ,2, .,10000000, ngen ,2, resorder ,0, .tppar ,1.7, uscale , none,.off );dl ,1, edim ,1);,dl ,4,5,6, edim ,1);dl ,1,

27、edim ,1);% Integrate I1=postint(fem,u,I2=posti nt(fem, K_x_ns I3=postint(fem,p, % WZ: write out output to the array outputoutput(j,1)=slot(j); output(j,2)=I1; output(j,3)=I2; output(j,4)=I3;output(j,5)=length(fem.mesh.p);% WZ: end our j-loopend% WZ: record results to filesave orifice.mat femdlmwrite

28、( orifice.dat ,output, ,)；三三 TJ7t -r t 三FT4Tt.二二二rmH:. :三5J2-: 一，二ff-sTF ，二三 MTrrgF.irnm组比2s E产打MlJll.W-74H二 457图4,速度场箭头矢量图及等压线（=0.05）该m文件脚本可能是本书中最长的，因此保留了部分注释。值得注意的第一点是几何结构是由脚本通过调整阻塞率生成的。通过建立练习中=0.4的几何对象。大多数的程序语句由 COMSOL Multiphysics GUI生成，选择“另存为 m 文件”进行保存。自行添加的主要是控制程序以及结果保存部分的代码。注意这里采用了自适应求解器，可以由

29、GUI来实现（将自适应求解器参数面板上的对 话框打勾即可）。对整个脚本的求解在一个常规的linux工作站下大约需要10CPU 分的时间。表2中对自由度的估计通过fem.sol.u长度来给出。表2孔板模型的主要输出结果。平均出口 x方向速度,积分粘性应力，自由度数目随孔板阻塞率的变化出口速度粘性剪切力p自由度0.051P -0.7852260,463905110.11-1.033661,551942350.151-1.210563,248100340 .0.21r -1.421165,60599151 10.251-1.554868,7291032600.31-1.824172.78810092

30、00.35 :1，-1.996378.043 |965790.4 11-2.320484.864 11064500.451-2.633893.8171030000.5 ；1 1-3.1446105.78 11041100.55 1-3.6672122.21 11038800.61-4.383145.46 ,1053400.651 1-5.2545179.9998070.7 I1-6.933233.9 1992310.751-9.2868325.69 ,92935表2由几何连续性研究的结果组成。明显的，附加压损随着阻塞因子的增加 而快速增加。粘性应力从 =0.05到=0.75也存在快速的增加。(

31、6)P231.0 0.61551.4726.955P0尽管（6）式在从=0.05到=0.75的范围内较好的拟合了所得数据，但如果在 一个较小的范围 氏0,0.5,对1 e的劳伦特级数可给出对数据更好的拟合。一 20.395652 0.3832P0(1)3图5和图6绘制了根据方程（5）和（7）得到的曲线以及计算所得数据点。由（7） 式得出的附加压损为 p=371,由（6）则为214,模型名&出309。（7）式的关键特点在 于如果考虑到式中分子中系数项几乎相等，该结构与量纲分析的结果基本一致。图5 Ap随e的变化（孔板厚度0.05）,曲线为&在0,0.5之间的三次曲线 拟合图6由1 e逆幕的劳伦特

32、级数进行拟合051 L.E图7 Re=1500, e =0.4。注意背流面粘滞流区域。流体的惯性导致流线跨过 该粘滞区域。在粘性流中，流体“转过拐角”并没有分离练习6.2Dugdale的孔板是锋利的。我们的孔板厚度则为整个流道厚度的 5%。试着 将孔板变得更薄：4%, 3%, 2%0对附加压损会有怎样的影响呢？根据文献12, 当缝隙宽度较小时，孔板厚度对曳力有一定的影响。如果缝隙是扁平的，则Dugdle 的量纲分析是正确的，如方程(3),但是如果颗粒是有限弧度的，则 Bungay和 Brenner的结果O(a-5/2)更为合理。练习6.3(a)改变孔板的上边界条件为对称边界条件。这个模型是二维

33、板状物的粘性绕流问题。尝试几何连续性。(b)改变m文件，不令上一次的几何结构作为下一次计算的初始条件。计算 快一些了吗？(c)这个例子并不是真正存在的物理现象，试着增加流函数-涡fangcheng用于考虑浮力对流的影响。平板问题由Kim13进行了研究，在缝隙宽度较小时采用了摄动法，将展开 的级数项加和得到曳力的奇异特性。练习6.4这些例子没有精确的采用几何连续性，这是因为由于Re=0几何形状改变成为线性问题，使得由于惯性力造成的非线性项随之消失。假设Re=1500.在fem.consts中给出常数，设置rho0=1500。图7中，注意背流面粘滞流区域。流 体惯性使得流体转过拐角，未造成分离。当

34、Re增大时，该流型结构与线性流有显著不同，求解较为困难。有两个有效的策略：对于Re的参数连续性或者几何结构角度的连续性。后者方法涉及到使用asseminit命令将旧的解内插到新的网fem.sol= asseminit(fem,init ,fem0)将 fem0.sol 使用内插方法转换到 fem 中的 网格上。由于新的几何对象和网格与旧的解和旧的几何对象并不兼容，因此这步工作是必要的。修改orfice.m文件以包含该步工作。3.瞬态几何连续性：膜干燥在之前章节中，除了考虑非线性情况(练习6.4)外，几何连续性不需要使用之前在不同几何对象上的解来作为新解的初始条件。几何连续性在一定范围的参数中使

35、用。在本章中，瞬态问题导致计算域将随时间发生变化，因此，上一时 间的解对于当前时间的解是非常重要的。这里所举例子为膜干燥问题，该模型在Mallego等人14实验工作的基础上做了一定简化。Exraporat&n图8膜干燥的两个前缘：顶部的蒸发前缘和内部的紧缩前缘包含颗粒(初始体积分数为0)的液膜在顶部表面以恒定速率蒸发。 如果颗 粒通过膜的扩散较小，则颗粒将在顶部表面聚集，堆积体积分数为m.随着时间 的增加，堆积层厚度也随着增加。整个膜的厚度则随时间线性下降。 计量垂直距离和时间随着初始膜厚度和特征干燥时间，膜表面可描述为y 1 To简单的质量平衡给出了聚集前缘以速度a移动，一匚。颗粒前缘位置可

36、得：m 0h； 1 Jf。颗粒传递如图8所示，在时间1/后没有进一步得聚集， pm 0这时达到了稳定的膜厚度。颗粒体积分数反映了颗粒前缘的变化， 可以很容易的 由阶梯函数表达：0 ( m 0)H(y hP)其中，H(y E)为海维斯特阶梯函数。图9朗谬尔等温吸附关系，k=10,A=1这也反映了膜的表面活性性能。当溶剂从膜中蒸发时，非挥发性的活性物质 保留下来，且既存在于溶液中也可以粘附于颗粒表面。假定初始容积表面活性物 质浓度为CS0,溶剂中以及颗粒表面活性物质浓度则与该初始浓度成比例。在颗粒填充动力学过程中，表面浓度由于聚集颗粒的吸附而不断变化。我们注意到溶液 中物质浓度的变化：初始条件：C

37、S 1C边界条件y 1:0,非渗透性表面无表面活性物质通过。y边界条件y 1 r：(1)T Pes(1)Cs一),y物质表面没有表面活性物质通过一一非挥发性物质停留在膜表面。H.F这里，Pes为表面活性物质的普朗特数，Pes 空E, H0是初始膜厚度，E Ds为上表面速度，Ds为表面扩散系数。为颗粒表面浓度，是心的函数。联系与 或称为吸附等温式。假设乳胶颗粒表面吸附遵循朗谬尔类型的吸附式，则吸附恒温式可写为其中的单位是moles/m2。因此以浓度的形式，我们可以获得无量纲恒温吸附式:CskCsACS0CsCs(8)其中，和A为朗谬尔等温吸附参数，k和A为朗谬尔等温吸附比例参数。图9以方程(8)

38、表示了朗谬尔等温吸附式。表面动态吸附式基础的，化学和石油化工中的常见情况。通过去污剂实现的 增强式油料循环已经应用了二十多年。吸附-脱附是液体色谱分离法的关键也是 另外一个例证。在紧缩前缘发生的解吸附过程是本模型的一个特点。Trogus等15 研究了油回收过程，提出了吸附/脱附速率的动力学模型。Ramirez等16对动态 吸附过程建立了双方程(浓度和表面活性物质担载)，1-D空时模型。他们的传递 模型仍然是针对均相多孔介质的，而我们的模型中，考虑到了紧密堆积和疏松堆 积层，表面物质的传递方程根据文献18给出：(Cs(1)fPesCs )T) y y其中，方程左手第二项代表了表面物质的聚集，第二

39、部分因子是由于吸附相的聚集。右手项代表了扩散。典型的参数值为m=0.64, 0=0.4图10坐标变换：一个前缘分析两个前缘问题，特别对于有效点源移动的情况特别复杂， 应尽可能通过 坐标变换至少消除一个移动前缘。我们也实验了使用非线性坐标变换在固定计算 域上消除两个移动前缘。令人惊奇的是这是可能的，但由于变化不是单调的，因 此大大增加了偏微分方程(9)的复杂性。较好的解决方案仍然是附着于一个前缘(内部的紧缩前缘)，将顶部前缘变换至固定域。变换较为简单：上；t1 t结果得出了对紧凑前缘的新的定义：p卜根据链式求导法则，微分方程可在新座标下得出:11 . ? . 、y 1 t t (1 t)(10)

40、(11)(12)可得变换后的表面物传递偏微分方程:d-Cs Cs1Pes(1)2(1Cs1)(Cs) m (1)p) 0(13)dCZ1这里，(p)是狄拉克6函数，出现在海维西特阶梯函数用于考虑颗粒体积分额。COMSOL Multiphysics提供了一个较容易的方法来使用海维西特函数以 及狄拉克6函数。两个m文件，flc1hs.m用于考虑海维西特阶梯函数，flc1dd.m 用于考虑狄拉克 6函数。表达式 Y=flc1hs(X,SCALE)通过在区间-SCALEX0)。解在0 T 1范围内较敏感。边界条件和初始条件分别为：初始条件：Cs(0) 1边界条件： &0,且有：0Cs (11PeS(1

41、)Cs)(1 )1表3包含了 COMSOL Multiphysics中对模拟的设置。应该注意这里使用的 COMSOL Multiphysics模拟策略提供了简单的方法来处理移动边界问题，且不需 要输出到m文件中。该策略避免了涉及修改 m文件。这里我们通过定义移动边 界为变量来处理移动前缘，然后利用内置函数来表示海维西特函数和狄拉克6函数。表3膜状干燥模型模型导航栏选择一维空间选择 COMSOL Multiphysics| Diffusion| Convection andDiffusion，瞬态分析完成Draw菜单设定物体直线0到1,完成Options 菜单 | Constants定义如下常数

42、：Theta_m 表达式：0.64Theta_0 表达式：0.40Alpha 表达式：theta_m/(theta_m-theta_0)Pe表达式：1K表达式：2A表达式：10Scale 表达式：0.001Options 菜单 |Scalar expressions定义表达式如下：Theta 表式：theta_0+(theta_m-theta_0)*flc1hs(x-frontpos, scale)Frontpos 表达式：(1-alpha*t)/(1-t) y 表式：fldc1hs(x-frontpos, scale) Ads_part 表达式：k*u/(A+u)Dtherm 表达式：k*A

43、/(A+u)A2)Physics 菜单 |Boundary settings选择边界1,设定逋量边界条件N0=-u*x*(1-theta+theta*dtherm)(1-t)选择边界2,设定逋量边界条件N0=(1-theta)*u+theta*ads_part)(1-t)-u* x*(1-theta+theta* dtherm)(1-t)Physics 菜单 |Subdomain settings选择子域1设定 涵=(1-theta+theta*dtherm)(1-t)A2设定 D=(1-theta)/Pe设定 R=-theta_m*(ads_part_u)y(1-t)设定 U=x(1-the

44、ta+theat*dtherm)(1-t)A2Init选项卡，设定u(t0)=1.0Mesh选择子域1,输入最大基元尺寸：0.001Solve菜单求解管理器，选择时间依赖求解器 设定输出时间0:0.001:0.374Advanced选项卡： weak形式点击工具栏上求解按钮(=)30图11不同干燥阶段表面活性物质的过量/消耗比率，当k, A和Pes值固定 时，它是距离基底距离的函数。图11为表面物质的分布。由 COMSOL Multiphysics模型得到的解给出 Cs 作为己的函数。但是在图11中，我们绘制了浓度百分比作为 己的函数，可表达 为：kC 。kC(m (1m)Cs)(1) (0

45、(10)Cs)|to)1O0(14)A CsA Cs(o -(1o)Cs) It oA Cs现在模型已被求解，但问题是解到底正确吗？这个问题可以从两方面来回 答。首先，如果解是正确的，图11中四条曲线代表了表面物质在四个不同时间 段的分布线。能够看到在溶剂蒸发阶段，表面物质分布在整个膜上并不连续， 而 是在颗粒前缘位置存在不连续。这个不连续是由于颗粒体积分数在颗粒前缘位置 的突变造成的。当颗粒体积分数最终达到均匀(=0.374, &=0),表面分布曲线变得连续。在整个干燥过程中，膜的上部具有较高的颗粒体积分额。当干燥进行过程中，表面物质浓度由于水的蒸发而增加， 结果，颗粒表面将发生更多的吸附作

46、 用。在颗粒体积分额以及颗粒表面吸附作用都增加的情况下，膜的上部达到饱和，因此，在p =0.374寸超过了空气水界面。证明模型精度的第二种方法是获得总的表面质量。表面物质是非挥发性的， 所以是守恒的。表达式(1 位 kCZ/(A 0)(1 )给出了表面物质在任意时间的总质量。因此，为了获得总质量，点击后处理标签，选择子域积分，输入(1 )u ku/(A u)(1)。在初始时刻你口末了时刻p =0.374寸计算积分COMSOL Multiphysics给出了表达式的数值，表明仅有大约 0.5%的误差。在较 宽范围内改变参数，误差均能够保持在 1%以下，这表明COMSOL Multiphysics

47、 在求解移动边界问题时的有效性。对于颗粒和表面物质的分布问题可参见我们之 前的出版物1718。4.结论本章中，我们研究了如何使用 COMSOL Multipysics模拟几何对象改变或者 由于前缘瞬态变化的情况，并介绍了这两方面的一些基础工作。特别的，在本章的前面部分运用了算子分割技术来处理瞬态几何对象改变的问题，几何图形在第一部分的时间步中发生，PDE方程则在第二部分中求解。这个较困难的方法要 求应用MATLAB自主编程，这里没有使用该方法，而是对内部移动边界使用了 平滑化界面的方法。根据质量守恒检验结果来看，这个新的方法是自洽的。止匕外，在许多情况下，可以通过在新的 ALE模式下定义移动边

48、界条件来避 免使用交错时间步方法来处理移动边界。两阶段方法要求在每一个时间步上重新 求解。因此，为了实现这个步骤，旧的解必须在新的网格上进行内插，这种情况 下模型的自由度可能不同，可以使用 asseminit()函数2合COMSOL Multiphysics 提供的代码处理这个问题。利用 MATLAB程序调用COMSOL Multiphysics函数 的方法是很常用的方法。在本章情况中，坐标变换和平滑界面方法也避免了使用 交错时间步方法。参考文献1 K. J. Ruschak, A method for incorporating free boundaries with surface te

49、nsion in finite element fluid-flow simulators, Int. J. Num. Methods Eng. 15 (1980) 639.2 R. T. Goodwin and W. R. Schowalter, Arbitrarily oriented capillary-viscous planar jets in the presence of gravity, Phys. Fluids 7(5) (1995) 954-963.3 H. Ockendon and J. R. Ockendon, Viscous Flow (Cambridge Unive

50、rsity Press, Cambridge, 1995).4 G M. Homsy et al., Multimedia Fluid Mechanics (CD-ROM) (Cambridge Unverisity Press, Cambridge, 2000).5 D. S. Dugdale, Viscous flow through a sharp-edged orifice, Int. J. Eng. Sci. 8 (1997) 725-729.6 W. N. Bond, Proc. Roy. Soc. 34(1922) 139.7 W. B. Zimmerman, The drag on sedimenting discs in broadside motion in tubes, I