（工程力学专业论文）波形板上薄膜流动的动力学分析.pdf

摘要 本文数值研究了波形板上薄膜流动的动力学特性。研究中考虑两种流动模型， 种为固定边界流动，不可压缩流体在固定波形底板和水平运动平板之间流动；另外 一种为自由表面流动，不可压缩流体在固定波形板上流动。选取两种波形底板模型： 一种为正弦板基面；另一种为三角板基面，其轮廓曲线分别为周期变化的正弦和三 角形曲线。流动中采用两种不同流体介质：牛顿流体和非牛顿流体。 运用有限元方法，从斯托克斯方程，连续性方程和流函数方程出发，通过有限 元离散、等参元欧射、数值积分和牛顿迭代，得到有限元方程的精确数值解，利用 后处理软件对得到的数值结果进行了可视化处理。本文分析了薄膜流动在不同情况 下的流线结构和底板上壁面剪切应力值，并定性研究了涡旋在流动过程中的生成及 演化情况。 在薄膜流动过程中，考虑波形板的轮廓形状参数变化、流动中雷诺数的大小、 薄膜厚度的变化以及非牛顿流体的剪切交稀系数对薄膜流动的流变特性影响。当波 形板的波动度超出临界值，波形板壁面上的流动出现流线分离，流动中有涡生成： 增大波动度，流动中出现第二个涡；月p = 0 ，流线在流场中对称分布；雷诺数增大， 流线驻点在水平方向出现位置偏移，涡的范围增大；对于非牛顿流体，随着剪切变 稀指数的减小，涡的作用范围减小。 本课题数值计算应用f o r t r a i l 程序w a v y p la 【e 和w a v ) r f i l m ；数据后处理部分使用 m a t l a b 编译器，对模拟计算得郅的数值结果进行图像可视化处理 关键词：薄膜流动；波_ 形板；涡生成；剪切变稀系数；波动度 a b s t r a c t u n d e rc r e e p i n gn o wc o n d n i o n s ，t h es t e a d ys t a t eo f t h i n6 l mn o wo v e rs t r o n 百yu n d u l a t e db o t t o m p r o f n e si sn 啪e r i c a l l ya n a l y z e di n 山i sp 印e ta st h eb o u n d a r yo f 恤d o m a i ni sd 旅r e m ，也e r e 对et w o k i n d so f m i nf i l mn o w ，o n eo f t h e mi sc r e 印i n gc o u e t t en o w ，a n dm eo 吐1 e ro n ei s 抒e e - s l l r f a c ef l o wf o r 疗e e - s u r f a c en o wm o d e l ，n u i dn d wo v e rf i x e d 哪d u l a t e dp r o 丘1 e ；f b rc o u e n ef l o wm o d e l ，n u i dn o w b e 似- e e nt 、op i a t e s ，o n eo f t h e mi sam o v i n gp l a l l e ，t h eo m e ro n ew i t hau n d 词a t e dp r o f i i e s i n u s o i d a l p m f i l ea 1 1 d 们a r i g l l l a rp r o n l e 甜ec h o s e n 站t h eb o 咖p r o f i l e ，a n dt h eu 1 1 d u l 砒e dp r o f i l ev 州e s p 甜o d i c a l l y n e 姒o n i a nf i u i da 1 1 dn o n 州e 谢o n i a nn u j da r ec o n s i d e r e dj ne a c hn o wm o d e l 1 1 1 es t r e a m l i n e so ft h ef l o wa r eo b t a i n e df b ma ne x a c ta 1 1 a l y t i c 8 l 仃e a t m e n to ft h es t o k e s e q u a t i o n sb ym a k i n gu s eo f 丘n i t ee l e m 曲tm e t h o d s a n dc o m p l e xf u n c t i o nt h e o r y b a s e d0 nm e c a l c u l a t i o n s ，t h ef b r n l a t i o na n de v o l u t i o no fv o r t i c e si n 1 en o wa r ev i s u a l i z e d f o rn 州0 n i a l ln u i d ， l ei n n u e n c eo f l l l ep l a t eu n d u l a t i o n s ，t h er e ”0 1 d s 加m b e r 锄dt 1 1 em i c k n e s so f t h en l ma r es h o w nb y p a 蛐g t e rs m d i e s ；蠡竹n o n n e 耐o n i a n 丑u i d ，出es h e 静t h i n n i n gc o e 蜀c i e n t 。nm e 盎h nn o wi s 曲o w nb y p 釉e t e r s t u d i e s f o 曲。也n e w t o n i a nn o wa n dn o n _ n e 、v t o n i a nn 0 w ，t h ev o n i c e sa r eg c n e r a t e dw h 锄也ew a v i n e s s o f 也eb o t t o mp r o 丘l ei sb e y o n dac 石c a ln u m b e lb yi n c r e a s i n gt l l ew a v i n e s s ，as e c o n dv o r t e xc a nb e g e n e r a t e d w i t h 也ei i l c r e 弱i n go fi 蛳o l d sn u m b e t l l es y m m e t f yo ft h es 心e 啪l i i l e sa t 忍= 0i s d e s t r o ”da n dt h ev o 而c e sa r es l o p e da n db e c o m e1 a r g e li na d d i t i o i l ，f 研1 en o n j n e w c o n i a nn o w ，t h e s 仃e a l l l l i n e sa r es y m m e t r i c a la n dt h ev o n i c e sb e c o m es m a l i e rw i t ht b ed e c r e a s i n go f s h e a 卜t l l i n n i n gr a t e n t h e s u a lf 0 m a i lc o d ew a v y p l a c e 鲫dw a v y f i l ma r c 印p l i e di nn u m e r i c a is i m u l a t i o n ，a n dm a t i a b m a t h e m 8 t i c a jl i b r a r yi su s e di np o s t p r o c e s s 啦g k e y w o r d s f i l mf l o w ；w a v y p l a t e ；v o r t e xf b n n a t i o n ：s h e a 曲n i n gr a t e ；w a v i n e s s - i i 独创性声明 本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取得的研 究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其 他人已经发表或撰写过的研究成果，也不包含为获得北京工业大学或其它教育 机构的学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何 贡献均已在论文中作了明确的说明并表示了谢意。 签名：蕴送日期：壬! 蔓! ! ：兰! ： 关于论文使用授权的说明 本人完全了解北京工业大学有关保留、使用学位论文的规定，即：学校有 权保留送交论文的复印件，允许论文被查阅和借阅：学校可以公布论文的全部 或部分内容，可以采用影印、缩印或其他复制手段保存论文。 ( 保密的论文在解密后应遵守此规定) 签名：换茁导师签名：窆l 丝i 整目期： 1 1 引言 第1 章绪论 研究薄膜流动在无孔介质基面上的涂抹或沉积的力学行为特性，对于许多工业 制造过程来说，具有十分重大的意义1 1 2 l 。薄膜的应用主要有以下几类：1 ) 装饰薄膜， 主要应用在手机l c d 保护屏、家用电器面板、数码产品的镜头保护屏以及各种其它 装饰性材料；2 ) 光学镀膜，比如说各种增透膜、防水膜、减反膜和多层膜等；3 ) 光电子镀膜，主要应用在i t o 透明导电玻璃、触摸屏用i t o 导电玻璃。许多工业产 品，光盘、磁盘、光学设备，特别是电子产品，如传感器、集成电路、显示屏面， 微电子设备等，通常需要涉及到多层的连续的薄膜涂抹。比如说，目前笔记本电脑 以及移动电话等热门电子产品都会因为高频电磁波之间的互相干扰，造成杂讯以致 影响通讯质量，而且人体长期暴露在强力电磁场下则易忠癌症。因此手机、p d a 、 携带式d v d 、掌上电脑都已经应用e m i 电磁波防干扰镀膜。手机按钮、按键类以及 l e n s 等，都采用装饰性金属镀膜，籍以提升产品质感。为了减少或增大光在元件表 面上反射，反射光学元件上将会涂抹一些光学增透膜或高反射膜，这类光学薄膜的 镀制是通过在元件表面沉积介质和金属材料如s i 0 2 ，t i 0 2 或铝，以达到预期的效果。 1 2 概述 薄膜涂抹的工艺流程大致可以分为以下几个步骤：1 提供涂覆用的溶液，这是 涂覆的准备阶段，根据涂覆预期达到的目的来调制涂覆溶液，选择的涂覆溶液是否 恰当直接影响着涂覆质量；2 将涂覆溶液涂到基底表面，这一过程有多种方法。比 如浸渍涂覆( d pc o a t i n g ) 、刮除涂覆( b l a d ec o a t i n g ) 、滑动涂覆( s l i d ec o a t i n 曲等；3 ，涂 覆溶液在基底表面自由分配，在涂覆过程结束之后，涂覆液体会在基底表面咀薄膜 流动的形式自由流动分配；4 形成均匀的薄膜，在某些涂覆方法中，最后还要通过 刮除等处理，确保最终形成的薄膜的厚度以及均匀程度符台预期的要求。其工艺流 程也可以简单概括为在器件表面涂覆一层或多层的均衡的液层，随后经过烘干，后 加工处理等来完成。薄膜通常用来保护物体器件或起装饰美化作用，同时也能做一 些其它的功用，比如说信息记录。在一般人眼中，薄膜涂膜通常与粉刷油漆和涂抹 颜料等联系在一起，由于日常的这种粉刷装饰应用，难免让人联想到薄膜涂膜研究 就是一种简易的，鞘细的，应用价值不高的工作。其实不然，在实际的工业应用中， 通常在薄膜涂膜过程中遇到以下的一些情况：薄膜厚度尺寸非常小而且尺寸精度要 求非常之高；薄膜流动的基面形状复杂多样；薄膜流动的流变特性在很大范围内变 化；由于某些涂抹的实际要求，流体特性不可以随意改变；几个离散的涂层需要同 化；由于某些涂抹的实际要求，流体特性不可以随意改变；几个离散的涂层需要同 时涂覆；高涂覆速度的工业生产要求；薄膜流动为自由边界等等因素的影响，实际 的薄膜流动是非线性的p 】。 现在，薄膜技术已经在工业生产上得到广泛的应用，而且，工业生产中涉及到 的薄膜技术一般都是多层涂膜技术，在实际的多层涂膜生产过程中，涂层厚度可能 会出现变化，这样，如果某一镀层的厚度涂抹得不均匀，将会导致整个薄膜的厚度 不均匀，甚至产生不稳定性，因此导致出现一些工业生产者们不希望看到的涂抹缺 陷。同时，在实际的生产过程中，由于一些机械设备、操作以及生产工艺等方面的 因素，薄膜涂抹层也将出现一些无法完全避免的缺陷，例如薄膜中出现气泡( b u d d l e ) 、 针孔( p i n h o l e ) 、波纹( w a v e s ) 、起脊( 脯) 等等诸如此类的各种缺陷 4 1 。随着对产品镀膜 质量以及外观完美度的要求的不断提升，如何控制或者减少缺陷的出现，就需要对 薄膜流动的动力学特性进行深入的研究。 在流体动力学中，关于薄膜流动这一领域的研究己经十分深入了。但是，绝大 多数的对此领域的研究都集中在对平面板( f l a tp l a t e ) 基面上的薄膜流动及其稳定性 的分析。s m l w a g o n 和l a r s o n 较早从事这方面流动的理论和实验的研究工作 5 - 7 1 。1 9 9 0 年，通过使用润滑近似理论，s t i l l w a g o n 和l a r s o n 最早计算出阶梯基面上薄膜流动的 一维润滑方程的解析解【”。1 9 9 2 年，p r i t c h a r d ，s c o t t 和t a v e n e r 对倾斜平面上的重力 作用流进行了分析【“。 对于辊筒之间的薄膜流动，模型见图1 1 ，g a s k e u 做了许多有实用价值的研究 工作，在试验中，采用光学器件和粒子成像技术，通过在流动中注射入颜料粒子， 从而观察流动的流线结构变化，研究两种情况，分别为两个辊筒相对正转和相对反 转的流动模型，研究中，考虑两个辊筒之间的距离、辊筒转速、流速以及毛细管数 对流动的影响【。在前面研究的基础之上，g a s k e l l 等( 1 9 9 8 ) 对反转锟筒 ( c o n 廿a - r o t a t h 培) 之间的定常牛顿流动进行深入研究，物理模型仍如图1 1 所示：辊 图1 1 辊筒薄膜流动示意图i 9 ) f i g 1 1r o l lc o a c i n g f 9 】 第1 覃绪论 筒的切线速度分别为u l 和u 2 ，半径分别为r 。和r 2 ，辊筒流动的薄膜厚度分别为马 和岛。在计算过程中分别考虑斯托克斯流的无量纲流动速度；两辊筒的剪切速度比 为弘u l 醍、毛细管数g 以及定常牛顿流动的薄膜厚度这三个参数的变化对薄膜流 动中的驻点一鞍点位置和流动模式的影响，整个研究结果通过实验和有限元数值计 算分析得出【1 0 1 。随后，g a s k e u 和s a v a g e 对三种实际工业应用中的斯托克斯流( m = o 的粘性流) 模型进行了详细研究，依据于一些流动假设条件，将流动模型转化为关 于流函数的双调谐边界值问题，同时给出了流动的流线结构图 1 1 】。 1 9 9 8 年，w i l s o n 和d u f 母采用润滑近似方法，研究了缓慢变化的基面上的单向 牛顿流，在研究过程中，考虑不同的横断截面【1 2 l 。同年，c h o 和k i m 对正弦波动壁 面上的流动进行了流动分析 1 3 。1 9 9 9 年，m a l 锄a t a r i s 和b o m o z o g l o u 通过使用迦辽 金法，对倾斜波动壁面上的定常牛顿流进行二维数值计算，研究中考虑两种不同形 状基面模型：“s ”型和“c ”型，研究流动对象为：层流、不可压缩、两维等温流。 其中考虑了高雷诺数和低雷诺数对流动的影响；考虑了壁面、自由表面的曲线形状 对流动的影响。同时对两个模型中流动区域中的流线结构、速度和压力分布进行了 分析1 1 4 j 。d e c r e ，f e m a i l d e z p a r e n t 和l 瑚e r s 在1 9 9 9 年再次对s t i l l w a g o n 和l a r s o n 【6 】 所研究的问题进行流动分析，为这种流动建立了格林函数格式 ” 。 1 9 9 9 年，国内学者李学文在长波解的情况下，对沿斜面流动的二层广义牛顿流 体的稳定性进行分析，推导出二层广义牛顿流体流动稳定性控制基本方程，研究了 广义牛顿流体在重力作用下沿斜面流动的稳定性，得到了在长波情况下，临界雷诺 数作为定常流动速度和平板斜率的函数。其中用到以下边界条件：1 ) 在自由表面处 切应力必须为零；2 ) 在界面处速度分量连续，切应力连续：3 ) 在固定的底板壁面 处流动无滑移【l 。对于双层流动，2 0 0 0 年，p a s c a l 对沿固定倾斜角度的具有中等雷 诺数的双层流动模型进行了深入研究，考虑重力流外层流体的厚度对重力流形状的 影响，以及一些其它因素如雷诺数、流体密度以及重力流流动的流变特性对重力流 形状的影响。h m t h o m p s o n 和n k a p l l r 对刚性辊筒反向转动之间的流动进行数 值计算，计算结果与试验结果十分吻合【1 8 】。2 0 0 4 年，j i a i l g ，h e l e n b r o o k 从全纳维 斯托克斯方程出发，通过数值计算，模拟分析倾斜表面上的无惯性双层薄膜流动的 非稳定性 1 9 】。 k a j l i a d a s i s 和h o m s y 使用润滑近似理论对具有沟槽形状基面上的薄膜流动及其 稳定性分析进行了详细的研究 2 0 2 l 】。对于在复杂形状基面上的重力作用薄膜流动， d e c r e 和b a r e t 做了很多有重要意义的研究工作 2 2 。 2 0 0 3 年，w i e r s c h e m 和a k s e l 对小波动度、低雷诺数下的倾斜波形板上的定常 牛顿流进行了线性稳定分析 2 3 。w i e r s c h e m ，s c h 0 1 l e 和a k s e l 在实验室中，利用荧 光粒子影像速度仪，对大波动度的正弦板上的重力驱动薄膜流动进行模拟试验，取 雷诺数为这1 0 。5 一数量级，研究在如此低的雷诺数下，波动度、薄膜厚度的改变对 流动特性的影响。 针对牛顿流动，传统的研究方法是采用润滑近似法，摄动方法以及数值模拟方 法。2 0 0 4 年，g a s k e l l ，j 疏a c k 和s e l l i e r 研究了无孔介质倾斜平板上的两维和三维重 力作用流的流动特性，运用润滑理论，通过多重网格方法，得到非线性控制方程的 精确数值解【2 “。 润滑近似法和摄动方法都有一个适用范围，仅仅在一定有效范围之内，才能得 出可以让人接受的结果。如果薄膜流动的基面极度不平，例如在波动度较大的波形 板上的流动，采用传统的理论方法将会失效。同样，如果波形板的波长超出一定范 围，因为收敛原因而会出现大量的问题，依据于无穷级数而得出的双调谐函数的经 典分解算法也将失效【2 6 】。 1 1 1 i e l e ，k n o b l o c h 对倾斜固体表面上的热薄膜流动特性进行了分析，并且考虑了 液气界面的热交换影响拉”。 2 0 0 4 年，s c h 0 1 l e 对波形板上的斯托克斯流进行了研究，考虑在两板块之间的定 常牛顿薄膜流动，两板块中，顶板为平面，底板为波形板。通过采用复合函数理论， 依据于正则函数，得到双调谐函数的解析解，将斯托克斯函数以及相应的边界条件 简化为单变量函数的一般积分和微分方程，从而将二维问题简化为一维问题，通过 傅里叶分析和代数计算得到方程的解，同时对波形板的波长以及两板块之间的平均 距离对流动特性的影响也进行了深入研究 2 8 】。 1 3 薄膜流动研究的一般方法 在流体力学中，针对固定波动板上薄膜流动的特性分析研究，通常会用到以下 几种不同的理论方法： 润滑理论：润滑理论属于粘性力起主要作用的小雷诺数流动。应用润滑理论， 可以将流动问题简化到一个仅仅关于局部膜厚度的方程，前提假设为在z 方向的典 型长度与在x 方向的典型长度的比值极小，速度域成为局部抛物线，对于周期变化 的波形板，前提假设就成为一个长波近似( 波长比深度大很多) ，仅仅适用于对于波动 度较小的情况。对于波动度大的波动板上的的液膜流动，润滑理论将不适用【2 9 1 。 摄动理论：摄动一词渊源于天体力学，但它同流体力学的发展有着密切的关系。 奇异摄动理论在流体力学中的早期应用，可以追溯到斯托克斯( 1 8 4 7 ) 关于非线性 深水波理论的研究。摄动理论由于不同的坐标系、不同的比例和摄动参数而相区别。 忽略级数展开中的高阶项，所以，摄动理论仅仅对于小摄动参数的情况才适用。即 波动板轮廓曲线的振幅远远小于波长，或者说，两板间的液膜厚度远远小于波动板 - 4 - 轮廓曲线的曲率半径。后一种情况要求一个正的局部倾斜角度作为其附加的几何约 束条件。当波长、底面轮廓曲面的振幅以及薄膜厚度在同一数量级上时，渐近展开 的级数不再收敛，所以说，大波动度波形板上的流动特性分析不可以采用摄动理论， 超出了摄动理论的适用范围。 数值方法：数值方法不受上面所提到的几何条件的约束，在计算流体动力学中， 所求解的都是一些非线性偏微分( 积分) 方程组。数值模拟的基本方法就是通过离 散方程来逼近描述流体运动基本规律的数学方程及其边界条件。从而通过采用一定 的方法来求解这些离散方程。离散的方法一般有：有限差分法、有限体积法、有限 元法和谱方法。由于在数值求解中将引入误差，如离散误差，求解方法中的误差( 迭 代误差、舍入误差等) ，即使所选择的数学方程是所研究物理问题的准确方程，所得 到的数值解也仅仅是按一定的精度来逼近准确物理解的近似解p 0 3 1 】。在数值计算中 经常遇到大型代数方程组的求解，而且计算问题越来越复杂，相应的代数方程组的 形式也越来越复杂，阶数越来越高，求解难度也越来越大。最常用的求解代数方程 组的方法是迭代法，但是迭代法一般来讲计算收敛比较慢，近十年发展起来的多重 网格法是一种加速收敛的比较好的方法，而且已经在计算流体动力学( c f d ) 中得 到比较广泛的应用。q 1 4 本文研究的目的和内容 在实际的工业生产中，薄膜流动的基面一般来说不仅仅有平面，而且还有一些 表面轮廓形状不规则变化的表面，本文将表面轮廓曲线周期变化的涂抹基面统称之 为波形板，关于波形板上薄膜流动的系统、全面的力学特性分析，以及应用有限元 理论方法离散求解流动有限元方程及其边界条件，这仍然是全新的课题。本文对这 方面的内容进行了比较深入的研究。 本文选取两种基面模型，即两种不同表面形状的波形板，波形板保持固定，基 面无限延伸，这样，可以将空间问题简化为平面问题，其轮廓曲线分别为周期变化 的正弦和三角形曲线，为了在以后的章节中方便介绍，这样的波形板分别称之为正 弦板和三角板。考虑牛顿流体和非牛顿流体在不同形状固定底板( 即正弦板和三角 板) 上的薄膜流动，这里，凡是应力和应变速度之间的关系满足牛顿内摩擦定律的 流体称为牛顿流体( n e 叭o n i a nf l u i d ) ，凡是应力和应变速度之间的关系不满足牛顿 内摩擦定律的流体称为非牛顿流体( n o n n e 叭o n i a l lf l u i d ) 。广义来讲，牛顿流体是 非牛顿流体的一种特例，对于非牛顿流体，主要考虑剪切变稀系数变化对薄膜流动 特性的影响，剪切变稀表示流动的剪切速率增大，流体的粘度变小，剪切变稀系数 是非牛顿流体的一个特性参数，剪切变稀指数等于1 的流体为牛顿流体。同时，对 于薄膜流动的上表面，本文也考虑两种边界，其中，一种薄膜流动的上表面为以匀 一5 一 速水平移动的平板，与传统意义上的牛顿库埃特流相类似，不同之处仅仅在于流动 的基面不是平板，而是形状周期变化的波形板，流动中不考虑重力的影响，这种流 动称之为固定板和运动板之间薄膜流动；另外一种薄膜流动的上表面则为自由表面， 液体在自由表面上与空气接触，与液体的粘度相比较，空气的粘度可以在研究中忽 略不计，底板与水平方向有一定的倾斜角度，流动过程中考虑重力的影响，这种流 动称之为固定板上自由表面薄膜流动。在薄膜流动的过程中，着重研究波形板的轮 廓形状变化、流动中雷诺数的大小、薄膜厚度的变化以及非牛顿流体的剪切变稀系 数对薄膜流动的流动特性影响。在研究波形板轮廓形状变化对薄膜流动的影响时， 引入波动度这样一个概念，对于正弦板，波动度为振幅与波长的比值；对于三角板， 波动度为三角形的高与底边的比值，波动度即为一个描述底板轮廓形状变化大小的 参数。在自由表面薄膜流动的研究中，采用的非牛顿流体介质为b c l o o o c s ，其物理 参数：密度为0 9 6 9 c m 3 ，运动粘度为1 2 1 3m m 2 s ，表面张力为2 0 ，4m n m 。为了 与实验数据进行对比，采用另一种流体介质：b c 5 0 0 0 c s 。b c 5 0 0 0 c s 的物理参数：密 度为0 9 7 2 c m j ，运动粘度为5 7 8 0n h n 2 s ，表面张力为2 1 4m n m 。 在研究中，运用有限元数值计算理论方法，从n s 方程、连续性方程和流函数 方程出发，通过有限元离散、等参元映射、数值积分和牛顿迭代，得到有限元方程 的精确数值解，利用后处理软件对得到的数值结果进行了可视化处理。分析了薄膜 流动在不同情况下的流线结构变化、底板边界上流动的流线分离情况以及漩涡在流 动过程中的生成以及演化情况，同时对底板上避免剪切应力值做了数值计算，根据 普朗特流线分离判据，验证固定底板边界上的流动流线分离现象的产生和变化结论。 2 1 引言 第2 章薄膜流动的数值计算理论、方法 在计算流体力学( c f d ) 中，求解流体力学微分控制方程的主要方法概括起来 有以下几种：有限差分法，有限元法，有限体积法，边界元法，有限分析法，谱方 法，和格子l a t t i c e b o l t z r n a n n 方法等。 有限差分方法是计算流体力学中最重要的数值方法之，它理论上系统成熟， 应用广泛、有效。采用有限差分方法数值求解问题，首先要求在求解域内，以差分 网格或差分节点将连续的求解域转化为有限的离散点集，所得的数值解是定解域内 网格节点上的离散数值解。 有限体积法从描写求解问题的守恒型控制方程出发，对控制方程在控制范围内 积分，运用高斯公式，按照一定的方法离散，因而导出的方程具有守恒型，而且对 区域的适应性比有限差分法好，是目前应用较普遍的一个数值方法，特别是数值传 热学领域。 边界元法的特点是可以降低求解问题的维数，只须求解计算区域边界上的未知 量，然后求得整个流场内的解，先求基本解，用格林公式导出区域边界上的积分方 程，然后利用有限元法的分段插值思想和过程，对边界积分方程数值求解。此方法 优点除维数降低一位外，精度一般比有限元高，但是其缺点也非常明显( 不适用非线 性问题，即对于雷诺数不等于o 的情况不适用) ，同时对一般问题难以找到基本解【l 】。 有限单元法是流体力学在数值模拟中广泛采用的种方法，究其原因有：迦辽 金( g a l e r k i l l ) 法的加权平均是一种最优近似法；能精确地处理不规则几何边界，因 为有限单元法对求解区域的单元划分没有什么特殊的限制，对流动中的未知自由表 面或不同介质之间的交界面也比较容易处理，这对模拟那些具有复杂流道的流动是 非常合适的 3 0 3 1 ，35 1 。 基于本课题的研究性质，在薄膜流动过程中，流体的惯性力不可以被忽略。流 体具有一定的粘度，因此，流动控制微分方程为非线性：流动底板的形状为周期变 化的波形板，当物理模型为自由表面流动时，自由表面形状不是平面，其表面的形 状是动态变化的，位置不可以假设为固定。因此，薄膜流动的研究过程中，采用有 限元数值计算方法，有限元方法是一种求解偏微分方程的有效数值方法，这一方法 首先需要在这些控制方程所在的空间对其进行离散化。所谓离散化就是将一个大的 计算区域划分成一些结构简单但形状任意的局部小区域( 即有限单元) ，单元的形 状一般为三角形或四边形，单元之间仅在边界上的节点处连接。单元离散化后，最 初的偏微分方程将变成具有某种形式的矩阵方程，它们把单元中某些特定点( 节点) 7 一 ，。一一，：。一型型鎏譬耋墼! ：二，一! ，。，。一，一 上的输入( 已知量) 与同一点上的输出( 待求解量) 联系起来。为求解这个大型计 算区域上的控制方程，可以将区域上的矩阵方程按节点叠加起来，得到某个总体矩 阵方程从而求解，也可以采用逐个单元法求解，以避免产生( 巨大的) 总体矩 阵3 9 0 。 2 2 流动方程 2 2 1 微分控制方程 在稳态条件下，不可压缩流的纳维一斯托克斯方程： p u v u = 一v p + p v2 u + p g 以及流体连续性方程： v u = o ( 2 一1 ) ( 2 2 ) 在上面的两个方程中，d 卢分别表示流体的密度和粘度，p 和【厂分别表示流体压力 和流动速度，其中g 为重力加速度。对以上两个微分方程进行无量纲处理，取d 为 长度比例因子，压力和粘性力比例因子为“，经过转换处理后，得到无量纲形式 的纳维斯托克斯方程和无量纲形式的连续性方程： r e h v 雎= v - 盯+ s 辔 v h = 0 ( 2 3 ) ( 2 4 ) 在无量纲方程( 2 - 3 ) 和( 2 - 4 ) 中，雷诺数船= 烈，d 4 堤惯性力与粘性力的比值，斯托克 斯数盯= 咫_ d 锄u 为重力与粘性力的比值，o = 哆n 【v 时( v h ) 7 】是应力张量，单元向 量蜃表示重力的作用方向。其中，j 为单位张量，p 和h = v ) 分别是无量纲压力和 速度。 2 2 2 有限元离散 将控制方程的空间进行离散化，即将流场区域离散为 有限数量的单元，每个单元上分布有一定数目的节点，节 点上的速度和压力为待求量，单元之间仅在节点上建立联 系。插值函数项的次数由单元上的节点数目来决定，因此 在有限元计算过程中必须慎重选定合适的单元类型。本文 中采用的单元类型是六节点三角形单元，每个节点上有三 - 8 图2 - 1 单元类型 f i g2 - 1 伊3e l e m e n t 第2 章薄膜流动的数值计算理论、方法 个自由度，即( “，v 。，脚，“、v 、p 分别表示x 方向的速度、y 方向的速度和匪力， i 为第f 个节点在每个三角形单元中，三角形三条边上的中间节点的压力自由度可以 通过其相邻角点的压力值插值计算得到，因此，中间节点上只有速度自由度。整个 单元具有压力和速度自由度的节点减少到3 个，只有速度自由度的3 个节点位于三 角形边的中点。本文中采用的单元类型详见图2 1 ，方框表示的是具有压力和速度自 由度的节点，圆点表示的则是只有速度自由度的节点。 通过引入形函数对连续函数进行离散化，得到速度和压力函数的有限元格式为 p b ，y ) = 芦，v ，b ，_ y ) ，一1 在这里，稚= ( “b 哟表示为单元上第五节点上的速度， 压力。 k 惭为插值函数，它们的数值满足为 ( 2 - 5 ) ( 2 6 ) p j 则为单元上第三节点上的 在节点址 在其它节点上 在节点工上 在其它节点上 ( 2 7 ) ( 2 - 8 ) 将( 2 5 ) 和( 2 6 ) 分别代入到方程( 2 3 ) 和( 2 4 ) 中，利用变分原理，通过加权余值法 近似得到有限元形式的控制微分方程，即关于误差值的方程。在数值计算中采用的 加权余量方法为迦辽金( g a l e r k i n ) 法，即把对速度和压力的插值函数作为权函数， 对于单个单元来讲，单元中有2 k + l 个自由度，动量方程和连续方程分别用坼，i 圻 作为加权函数，即有2 k + l 个非线性代数方程。经过迦辽金法加权处理以后，动量误 差和连续性误差分别为 硝= 【 e 。0 v 扣。v 口一m j 嗜地= o ( 2 9 ) 磁= k f ，v 髓d q = o ( 2 - 1 0 ) 在动量误差方程( 2 9 ) 中，代入等式胍v 盯= v ( 肌0 一砜历这样，上面的动量 误差方程转换为 置盏= 【陋e 。如v k v 一帆盯) + v 。盯一。s 嘻脚 ：t 陋州。( 川k + 眦一m 孵妞一“确西 2 1 ” 、jybi 一 。h i i 、jy0 h 1 0 l 0 ，1l r，、l = = 妒 程方元臣n口p有 j21 【 从误差方程( 2 1 1 ) 中可以看到，动量误差积分式由两部分积分组成，一部分为面积积 分，另一部分为流场边界线积分。对于固定板和运动板之间薄膜流动来说，固定边 界满足相对无滑移条件，同时，由固定壁面上的边界条件和运动壁面上的边界条件 可以知道，固定底板壁面上的流体速度为零，接近运动平板的流动速度和运动平板 的移动速度一致，而且整个流场的区域形状在流动过程中没有发生变化，所以，边 界上节点的( 蜥，v l ，p ，) 为确定值，因此与这些边界上的节点相对应的方程将会在数值 计算过程中被抛弃，对于动量误差方程的求解来讲，可以不需要求解关于线积分的 那部分。而对于固定底板上自由表面薄膜流动来讲，液体与空气的交界面( 即自由 表面) 上的节点的自由表面的法线方向的速度已知( 空气的粘度忽略不计) ，但是 水平方向上的速度为未知量，所以需要考虑式中关于流场边界的线积分。 对于将要研究的固定底板和运动平板之间的薄膜流动，完整的动量误差方程和 连续性误差方程整理为 砖，= h 帆m 警+ 叭脚，警卜芦，警 、 。 7 、( 2 - 1 2 ) 材。盟盟盟盟忆，盟盟k ：o 砂砂1 缸砂1 缸缸j 月乱= h m 蛳，警川v ，m ，警卜芦，警 、 、 、( 2 - 1 3 ) ，盟掣帆盟盟+ 2 v 。盟盟k ：o 四劣鲫积卵鲫i 耻酬虬警警肛。 ( 2 - 1 4 ) 对于固定板上自由表面薄膜流动，由于自由表面节点上的自由度未知，在积分 方程中需要考虑边界条件对积分方程的作用，方程( 2 1 1 ) 中的线积分部分必须在计算 中加以考虑，对于固定板上自由表面薄膜流动的动量误差方程中的线积分部分，将 在下面的章节2 4 中进一步讨论给出。 2 3 有限元方程的求解 2 3 1 等参元映射 在直角坐标系中，对于一个已经划分完毕的有限元网格来说，网格中所有三角 形单元的尺寸和形状一致的时候，选取的插值函数就会比较简单，但是在有些实际 的流动问题中，要求划分的网格单元形状一致，是不切实际的。因此，在整体直角 - 1 0 坐标系中，插值函数的表达式就比较复杂。为了简化计算，求解过程中，引入等参 元映射方法，即选用一种自然坐标，通过坐标变换处理，将所有的三角形单元映射 为规则的三角形单元，如图2 2 。 在等参元映射处理过程中，引入局部面积坐标( 自然坐标) 上l ，三2 ，上3 ，在三角单 元中，三角形内任意一点与三个顶点连线，单元被划分为三个子三角形，这样，三 角单元中的点就可以通过这三个面积坐标来确定，局部面积坐标工1 ，e 2 ，厶分别 定义为 三l = 4 1 埘，三2 = 爿2 0 ，3 = 3 ，4( 2 1 5 1 图2 2 三角形单元的映射 f i g 2 21 1 1 ei s 叩a r m e m cm 婶p i n g 其中，爿是三角形单元的总面积，4 l ，一2 ，彳3 分别为三角形单元三个顶点相对应的 予三角形的面积，三角形单元的总面积：爿= 爿1 十4 2 十爿3 ，即面积坐标上1 ，工2 ，三3 中只有两个是独立的，且具有与插值函数相等的性质，见图2 。3 。 图2 - 3 局部面积坐标系 f i g 2 3s u a 1 1 9 l e su s e dt od e f i n em ea r e ac o o r d i n a t el o c a la r e ac 0 0 r d i n a t e 选用如图2 - 2 所示的自然坐标和相同的节点标号，局部面积坐标上1 ，三2 ，厶 的坐标表达式为 厶= ；( 1 蟛) ，厶= ；6 一掌+ 娲l 厶= ；6 一毒一向) ( 2 - 1 6 ) ，。一，。：一。，一一。竖! ! 型童鎏! 圭! ：些0 ，一! 。，一。，。： 通过坐标转换后，单元上六个节点上的插值函数在局部坐标系( 自然坐标) 中 的表达式为 - 2 厶( 2 厶一1 1 。5 e ：( 2 厶一1 l 虬2 厶( 2 ，一1 ) r 2 1 7 、 4 = 4 上l 三2 ，j = 4 2 3 ，d = 4 l 三3 、。 在一个单元中，所有节点上的插值函数的代数和等于单位1 。即 ， m 皓，彳) = 倍，7 ) = 1 ( 2 一1 8 ) 女= l，= l 等参元的含义就是在所有单元中利用相同的插值函数来定义单元的几何形状和位移 场，比如 x = k m 偕，7 ) ( 2 - 1 9 ) 这里的巩是第个局部节点的直角坐标。 在直角坐标中插值函数是关于x ，y 变量的函数，而自然坐标中的插值函数是关 于f 和_ 1 7 变量的函数，引入等参元后，积分计算也应从直角坐标平面转换到自然坐 标平面上，在以后的计算中必须对两种系统进行统一处理。这里，先对f 和目进行 求导， a 鸳 a a 打 之， 卦 在上面的等式中，l ，为雅可比矩阵，t ，的系数由下面的方程来决定 因此， j = 缸 钙 叙 孙 劫 芘 咖 _ _ 一 卸 嚣= 参等，骞= 扣等a 善f ：f 。a 善7a 7f ： a 玎 通过采用等参元映射，整个流场离散成 g ，芦l 劬2 蒌 l o ，g 代，1 ) ，】，i 鹾加j 其中，式中卢为未知变量向量卢7 = lv l p 勺 - 1 2 ( 2 2 0 ) ( 2 乞1 ) ( 2 2 2 ) ( 2 2 3 ) 、_llj_-，_tlj a一缸a一砂 w i i h i i 八砂一鸳分一卸缸一西鱼却 第2 章薄膜流动的数值计算理论、方法 2 3 2 数值积分 前一节已经推导出有限元控制方程的离散形式，对于离散形式中的积分项，需 要采用一种数值积分方法对得到的基于单元的积分式进行求解，这里采用一种有效 的数值积分方法g a l l s s 积分( 高斯积分) ，在有限元方法的数值计算中，g a u s s 积 分是一种经常采用的数值积分形式。这种数值积分方法就是利用函数在若干点的函 数值，近似的计算该函数在某个区间的定积分 4 ”。求积公式中的插值节点称为g a u s s 点，取的g a u s s 点越多，得到的计算精度就越高，单元刚度矩阵会出现过分刚性， 当然，对计算量也有较大的影响，取的g a u s s 点较少，可能会引起不稳定，所以在 g a u s s 积分，需要多方面考虑来选取合适的g 锄s s 点数目。通过应用g a u s s 积分方法， 上一节中的积分式可以转换为 【，g 刃1 ，l d 影叼= 厂慨靠h ( 2 2 4 ) 1 l = l 式中是对应于第f 个高斯点的权数，h 表示高斯点的数目。对于将要研究两种薄膜 流动，在标准的规则三角形单元中，对于前面章节中推导出来的动量误差积分方程 的求解中，多方面综合考虑来选去合适的高斯积分点的数目。 在正弦板上薄膜流动模型中，无论是正弦板和运动板之间薄膜流动( 固定边界 薄膜流动) ，还是正弦板上自由表面薄膜流动，在有限元数值计算过程中，本文都采 用相同的起始网格划分，即采用相同的网格数目和相同的网格分布，具体的网格分 布详见图2 4 。在对正弦板上流动区域进行划分的有限元网格中，x 方向最多有1 2 5 图2 - 4 有限元网格示意图 f i g 2 4t h eg r i du s e di nl h ef e m 1 3 个节点，y 方向上最多有6 1 个节点，整个网格中一共划分有1 8 4 8 三角形单元，每个 单元上有6 个节点，每个节点上具有3 个自由度：“、v 、p ，分别表示x 方向的速度、 y 方向的速度和压力。这样，网格中总共具有3 8 2 9 个节点和8 1 4 9 个自由度，对于三 角板上的薄膜流动，包括固定边界薄膜流动和自由表面薄膜流动，网格划分有所变 化，具体的两种波形板上薄膜流动的网格基本划分信息如表2 ，1 所示。 表2 1 网格信息 t a b l e2 - lt h ei n f o r n l a d o no f m ef i d 正弦板薄膜流动三角板薄膜流动 x 一节点数( m a x )1 2 51 6 9 y - 节点数( m “) 6 l1 0 1 单元总数1 8 4 84 9 2 2 单元节点数66 节点自由度 33 节点总数 3 8 2 91 0 0 3 9 总自由度 8 1 4 92 2 4 6 8 针对本文研究的固定边界薄膜流动和自由表面薄膜流动，在标准的规则三角形 单元中，经过多方面综合考虑，对于前面章节中推导出来的动量误差积分方程的求 解，在关于区域面积积分部分中选取四个g a u s s 积分点，在线积分部分中选取3 个 g a u s s 积分点，具体的面积g a u s s 积分点和线积分g a u s s 积分点在自然坐标系中的位 置以及对应的权数分别见表2 2 和表2 3 。当在高斯积分公式中采用4 个g a u s s 积分 点时，其计算代数精确度可以达到9 。 表2 - 2 面积积分中高斯点坐标以及其对应的权数 t a b l e2 2g a u s sp o j n t s 柚dc o r r e s p o n d m gw e i 曲tf o rd o m a i ni n t e 舭l s 善野 w 110 j 1 6 21 ，2 历2 j 1 6 31 2 一i 2以1 6 4 00 9 3 1 6 1 4 表2 3 线积分中高斯点坐标以及其对应的权数 t a b l e2 - 3 g a u s sp o i m sa dc o n 它s p o n d i n gw e i 曲t sf o r b o u n d a r yi n t e g r a l s 工】 w 1