 Applied Physics 应用物理, 2012, 2, 35-40 http://dx.doi.org/10.12677/app.2012.22006 Published Online April 2012 (http://www.hanspub.org/journal/app) Generalized Difference Methods for a Fluid Mixture Model Suli Lin, Quanxiang Wang, Zhiyue Zhang Jiangsu Provincial Key Laboratory for Numerical Simulation of Large Scale Complex Systems, School of Mathematics Science, Nanjing Normal University, Nanjing Email: linsuli2008@163.com, zhangzhiyue@njnu.edu.cn Received: Jan. 18th, 2012; revised: Feb. 24th, 2012; accepted: Feb. 29th, 2012 Abstract: In this paper, we propose two numerical methods for a fluid mixture model. The model is usually used to describe the tissue deformations. It contains a nonlinear hyperbolic equation and an elliptic equation. The first numerical method is the generalized difference method based on linear element function space and piecewise constant function space. Numerical experiments show that our scheme is second-order accuracy in space. To eliminate the oscillation near the discontinuities, we design a generalized upwind difference method to solve the fluid model. Numerical results show that the two methods are effective for the considered fluid mixture model. Keywords: Fluid Mixture Model; Generalized Difference; Upwind Schemes 一类流体混合模型的广义差分法 林素丽，王全祥，张志跃 南京师范大学数学科学学院，大规模复杂系统数值模拟江苏省重点实验室，南京 Email: linsuli2008@163.com, zhangzhiyue@njnu.edu.cn 收稿日期：2012 年1月18 日；修回日期：2012年2月24日；录用日期：2012年2月29日 摘 要：本文针对一类流体混合模型设计了两种数值格式。该流体混合模型是关于可收缩间叶细胞组织变形的 模型,是由非线性的双曲型方程和椭圆型方程组成的混合方程组。第一种方法通过选取试探函数空间和检验函数 空间为一次元函数空间和分片常函数空间,针对光滑情形,得到的广义差分格式具有二阶精 度。为消除解在间断 处的数值震荡,我们设计求解该流体混合模型的广义迎风差分格式。数值结果表明两种数值方法对考虑的混合模 型是有效的。 关键词：流体混合模型；广义差分；迎风格式 1. 引言 过去几十年里，人们已经很好地把连续介质理论 和固体流体运动变形的原理结合起来，应用于生物学 组织和细胞的研究中。这些研究集中于组织的生长和 建模、力学、非反应混合物的传播以及组织工程学上。 1986 年，Dembo 等人[1]基于反应的内部敏锐流这一概 念，提出了关于可收缩生物聚合体网状系统的物理模 型。1996 年，Manoussaki 等人[2]提出了体外血管网络 构成的一种机械模型。近年，又有不少作者基于人类 软组织变形建立了其它模型[3,4]，然而数值计算方面的 研究却很少。Lubkin 等人[5]利用带有一步预处理–校 正的有限差分方法来求解输送细胞和纤维的体积分 数及新生物的细胞–纤维态分数。He 等人[6]在带有自 适应网格的圆柱几何体上应用了有限元方法对表面 张力梯度作用下石油液滴的变形进行了数值模拟。 Sta st na[7]利用显 式有限差分格 式计算了其中 与时间有 关的情形。 本文研究的是一维可收缩间叶细胞组织变形的 Copyright © 2012 Hanspub 35  一类流体混合模型的广义差分法 数学模型，该模型的建立参见相关文献[1,2,5]。另外 Jiang 等人[8]给出了下面三个假设： 1) 间叶组织细胞是由水成和细胞纤维状态组成； 2) 水成状态具有Stokes 流的特征； 3) 细胞纤维状态的压力可由相关时间尺度上的 永久变形来消散，而且可当做Stokes 流来处理。 在此基础上建立的一维数学模型为： 0, 0, v x L tx (1) 10, pv xx (2) 2ln1 v Mp xx 0, (3) 其中 表示细胞和纤维的体积分数， v表示细胞纤维 态的速度， 表示牵引系数， 是压力，p 表示收缩 系数， 表示膨胀系数， M 为细胞–纤维分数的粘 性系数，由此可知，01 。相关的边界条件如下： 0,, 0,0,, 0, 0,, 0. pp vt vLttLt xx tLt xx (4) 利用边界条件(4)，并从方程组(1)~(3)中消去 ， 原方程变为： p 0, 0, v x L tx (5) 2ln 1 vv M xxx x 1 (6) 本文的目的是针对上面简化的方程组进行数值 求解，来得到不同方法求解该流体混合模型的有效 性。该方程组是由一个类似对流的双曲方程和椭圆型 方程组成。我们将用广义差分方法求解该问题。此方法 于上世纪 80 年代提出[9]，因为广义差分方法简单，计算 精度高，目前已被广泛应用于计算流体力学中[9-13]。 本文结构安排如下：第一部分利用广义差分方法 对方程组(5)，(6)进行离散求解，试探函数空间采用一 次元空间，检验函数空间采用分片常数函数空间。给 出了两个不同模型的计算来验证方法的有效性。第二 部分采用广义迎风差分方法对方程组进行求解，同样 给出两个不同模型的数值实验来验证方法的有效性。 2. 广义差分方法 本节将采用一次广义差分方法对方程组(5)，(6) 进行求解。 2.1. 广义差分格式 首先对区间 0, L进行剖分 ，节点为： h T 01 1 0NN x xxx L 剖分 1 :,,1,2,, hiiii TIIxxi N ，单元 i I 的长度为： 1,1,2,, iii hxxiN 。再做对偶剖分 h T ，节点为： 012 12 0NN x xxx L ， 12 12 00 1212 ,,1,2,, ,,,, hi ii NNN TIxx iN IxxI xx 1, 其中 121 2 ii xxx 。现取试探函数空间 为相应 于剖分 的一次元空间，。与 h U I h T 0 hh Uuc i x 对 应的基函数为： 1 1 1 1, 1, 0, i ii i i ii i xx , , i x xx h xx x xxx h 其他. h 所以 h uU 可以唯一表示为 n h i uui i ，其中 ihi uux。因此在单元i I 上有： 11 1 hiii ii uuxxhuxx h i 1,,1,2,, hiii i uuuhxIi N . 选取检验函数空间为 0 0, , hhNh i Vvx xIxIvI 或在每个上为常数 ， 所以 的基函数可取为： h V 1,, 1, 2,,1. 0, , j j j xI xj xI N 对方程(5)两边分别乘以检验函数 j ，然后在区 间 0, L上积分，可得下面的离散格式： 11 11 d2 jjj jj Ixv v t 0 (7) Copyright © 2012 Hanspub 36  一类流体混合模型的广义差分法 对(7)式时间导数项采用向前差分离散，可得方程 (5)的全离散格式： 11 11 1 2jj jj nn nnnn jj tvv h , (8) 其中表示时间步长， h表示这里采用等距剖分。 t 类似地，可得(6)式的广义差分格式： 11 11 11 11 11 11 24 2 1 1ln 1 2 1 ln1. 2 nnnnn jjjjj nnn j jj n j nn nn jj jj nn nn jj jj MMh M vvv hh h 1 (9) 对于边界条件的处理，我们将通过引入虚网格点 来构造二阶精度的格式[14]，从而得到了求解流体混合 模型的广义差分格式。 因为 0 ，所以易知上式中 v的系数矩阵是 三对角阵，且是对角占优的，从而椭圆型方程的解 n j v 存在而且唯一。求出n j v的值代入(8)式即可求出 1n j 的 值，求出的 1n j 又可代入(9)式求出1n j v。如此循环就 可以把未知量 和在节点的值全部求出。下节，我 们将基于上面给出的离散格式给出不同模型的数值 结果，从而说明我们格式的有效性。 v 2.2. 数值例子 因为此处考虑模型的分析解并不知道，我们将用 半网格原理[15]来说明新构造格式的精度和收敛性。首 先定义以及 范数： L2 L 2 ,max 2 hh Lhtui ui ， 2 22 1 ,2 N hh i Lhtuiui h . 定义如下的误差比率： 2 2 hh hh uu ruu ， 这里 表示u 或v，表示空间步长为 时的数 值解。如果趋近于 2说明格式为一阶精度，如果趋 近于 4说明格式是二阶精度。 h uh r 模型一 根据 Jiang 的研究[8]，我们取 ，0.5M0 ， 并且假设考虑的区间为[0 。相应的边界条件为： , 1] 0,1,0,0,1, 0vt vttt xx . 数值实验一 对方程(6) 中的其余参数分别取为 1.0 ， 1.0 ，取初值 01212sin 2π x 。表 1给出的 是在 0.05T 时刻 和的误差，其中满足v2 th 。 从表中可以看出，随着 成倍的增加， N 两种误差范 数均相应减小，而所定义的范数的比率逐渐趋于 4， 另外，观察 的误差范数，也能得到类似的结论。这 表明广义差分方法解模型一是 可行的 ，且使 v 和 都 达到了二阶精度。 v 数值实验二 对于方程(6)取1.8 ，1.0 ，初值为： 0 0.1, 0.650.35, 0.5, 0.350.65. xx x 或 用一次广义差分方法求解的 ， 5T200N 时， 解的图像见图 1和图 2。从图 1可以看出，这种取值 下解 在中间位置向上增长，但 的解随着时间发展 在间断处出现了很大的震荡。另外，从图 2可以看到， 由于 在间断出现了较大的震荡，导致 v的解也出现 了一些震荡，这与实际情况不符合。其中图中的虚线 表示初始值，实线表示T1, 2, 3, 4 中间时刻解的情况， 点线表示 5T 时的解。这表明广义差分法不能消除 解在间断处的震荡。为了能够对两个相邻组织进行真 实数值模拟，我们将在第二节构造适合求解间断问题 的数值格式。 模型二 此模型中取 0.5 M e ，e 。 数值实验三 对方程(6)中其它的参数分别取为 2.0 ，1.0 ， Table 1. Error norms of θ and v using generalized difference method at T = 0.05 表1. T = 0.05时，广义差分方法得到的 θ和v的误差范数 θ v N L∞范数 L2范数 L∞范数 L2范数 406.3769e–004 2.9506e–004 1.8009e–004 1.1243e–004 801.6427e–004 7.6023e–005 4.3977e–005 2.8012e–005 1604.1326e–0051.9155e–005 1.0936e–005 7.0304e–006 3201.0340e–0054.8010e–006 2.7447e–006 1.7735e–006 Copyright © 2012 Hanspub 37  一类流体混合模型的广义差分法 Figure 1. The change of θ at T = 5 and N = 200 图1. T = 5，N = 200时，θ的变化 Figure 2. The change of v at T = 5 and N = 200 图2. T = 5，N = 200时，v的变化 初值为 01314sin 2π x 。边界条件为： 0,1, 0,0,1, 0vt vttt xx . 表2给出的是在 时刻0.05T 和v的误差，其 中满足 。从表中可以看出，随着 成倍的增加， 2 th N 两种误差范数均相应减小，而所定义的范数的比率 逐渐趋于 4，另外，观察 的误差范数，也能得到类 似的结论。这表明广义差分方法解模型一是可行的， 且使 v 和v都达到了二阶精度。 3. 广义迎风差分方法 为了求解初值间断问题，本节将采用一次广义迎 Table 2. Error norms of θ and v using g meth T = 0.05 风差分方法对方程(5)进行求解。 eneralized difference od at 表2. T = 0.05时,广义差分方法得到的 θ和v的误差范数 θ v N L∞范数 L2范数 L∞范数 L2范数 401. 6. 2.4 1.44930e–48390e–55259e–003447e–00 803.9073e–5 1.7835e–5 6.1802e–005 3.2485e–005 1609.8752e–6 4.5064e–6 1.5287e–005 7.9692e–006 3202.4820e–6 1.1304e–6 3.7802e–006 1.9694e–006 .1. 广义迎风差分格式 元函数空间 为对应 于对 3 对于第二节的剖分，设有限 h U 偶剖分 h T 的一个分片常函数空间。则对任意的 hh U 有下面的表达式： hi i n i ，其中 1, ,xI ()0,1,2, , 0, , j j j x jN xI , , hjj x tt 显然 。在方程(5)两边分别乘以任 意的 j h U ,并且在[0, L]上积分有： dd jj j hh jjh II I d0 j x tx x v xvx . 由格林公式可得 d j j h jhj I I vx v x , 其中 j I 表示对偶单元 j I 的边界。根据迎风格式的方 法，当 0 j时，有 v 1 d jj h jhj jjj II xv v x 当0 j v 时，有 1 d jj h j hjjj j II vx vv x . t 设是时间步长，利用向前差分公式可得全离散 的广义迎风格式： 1nnn v 11 1 11 0, 0 0, 0 n jjj j nnnn n jjjj j nnn n jjj j nnnn n jjjjj vv v th h vvv v th h (10) 结合离散格式(9)，就可以对考虑的流体混合模型 进行求解，其中对边界条件的处理按照与第二节同样 Copyright © 2012 Hanspub 38  一类流体混合模型的广义差分法 的方 3.2. 数值例子 然使用半网格原理来说明广义迎风 效性。 研究[8]，我们取 ， 式，数值算法与广义差分格式的类似。下面我们 将用广义差分方法的数值实验来说明推导出的广义 迎风格式的有效性和精度。 本节我们仍 格式的有 模型一 根据 Jiang 的0.5M0 ， 并且的区间为 。相 假设考虑 [0, 1]应的边界条件为： 0,1, 0,0,1, 0vt vttt xx 数值实验四 对方程(6) 中的其余参数分别取为 1.0 ， 1.0 ，取 初值 01215sin 2π x 。表 的3给出 是在 0.05T时刻 和v的误差，其中 t 表中可以 增加， 满足 从 ，随 成倍的 h 。 看出着N 两种误差范数 小，所义的范数的比率逐于2，另 外，观察 v的误差范数，也能得到类似的结论。这表 明广义迎风差分方法解模型一是可行 的，且使 均相应减 而渐定 趋 和v 都达到了一阶精度。 数值实验五 对于方程(6)取1.8 ，1.0 ，初值为： 用一次广义差分方法求解的 时解 的图像见图 3和图 4。 00.5,0.350. 0.1, 0.65 0.35, 65. xx x 或 5T，200N 从图 3可以看出，这种取值下解 在 置向 上增长，但 中间位 的解随着时间发 没有出现 数值 表3. T = 0.05时,广义差分方法得到的 θ和v的误差范数 θ v 展在间断处并 v震荡。另外，从图 4可以看到， 的解也没有出 现非物理的震荡，这与实际情况基本一致。其中图中 Table 3. Error norms of θ and v using generalized difference method at T = 0.05 N ∞2∞范数 2 L范数 L范数 L L范数 40 7.5719e–004 3.9250e–0043.6020e–004 1.7796e–004 80 3.4 1.42.4 1.47693e–009702e–008765e–004722e–00 160 1.8920e–004 9.8728e–005 1.6830e–004 8.7304e–005 320 9.4605e–005 4.9424e–005 8.9982e–005 4.6959e–005 Figure 3. The change of θ at T = 5 and N = 200 图3. T = 5，N = 200时，θ的变化 Figure 4. The change of v at T = 5 and N = 200 图4. T = 5，N = 200时，v的变化 的虚线表 时刻解 情况，点线表示 示初始值，实线表示 1, 2,3, 4T中间 的5T 时的解。这表明我们构造的 处的震荡 行数 广义迎风差分法能消除解在间断 ，适合对两 个邻近的细胞组织进 值模拟。 模型二 此模型中取 0.5 M e ,e 。 六 其它的数分别取为 数值实验 对方程(6) 中 参2.0 ， 1.0 ，初值为 01314sin(2π) x 。边 ： 界条件为 Copyright © 2012 Hanspub 39  一类流体混合模型的广义差分法 Copyright © 2012 Hanspub 40 Tifference 表4. 广义差分方法得到的 θ和 able 4. Error norms of θ and v using generalized d method at T = 0.05 T = 0.05时, v的误差范数 参考文献 (References) θ v N L∞范数 L2范数 L∞范数 L2范数 40 2.4098e–004 1.2476e–0042.7195e–004 1.4290e–004 80 1.1958e–004 6.57.5 3.2457e–007477e–008577e–005 160 5.8638e–005 3.0999e–005 2.4488e–005 1.1629e–005 320 2.8928e–005 1.5410e–005 8.7890e–006 4.0442e–006 0,1, 0,0,1, 0vt vttt xx . 表4给出的是取不同的空间步长和时间步 用广义迎风差分格式得到的数值实验六的误差结果。 从表 长时， 中我们可以看出，当 th 时，随着空间步长和 时间步长都几乎变为原来的一半，相应解的误差范数 都变为原来的12，这说明 度为 Oth 。 4. 结论 解的精 本文针对 察细胞组织 混合流体模型，设计了两种数值格式考 的变形，数值结果表明，标准的广义差分 格式 自然科学基金(11071123)以及江苏 金(BK2009397)的资助。 [1] M. Dembo, F. Harlow. Cell motion, constractile networks, and the physactive flow. Biophysical Journal, B. Vernon, et al. A mechanical omputers nite volume element s for problems with initial and boundary layers. ics of interpenetrating re 1986, 50(1): 109-121. [2] D. Manoussaki, S. R. Lubkin, R. model for the formation of vascular networks in vitro. Acta Biotheoretica, 1996, 44(3-4): 271-282. [3] G. Szekely, C. Brechbuhler, R. Hutter, et al. Modelling of soft tissue deformation for laparoscopic surgery. Simulation, Medical Image Analysis, 2000, 4(1): 57-66. [4] W. Mollemans, F. Schutyser, J. Cleynenbreugel and P. Suetens. Tetrahedral mass spring model for fast soft tissue deformation. International Symposium on Surgery Simulation and Soft Tissue Modeling, 2003, 1673: 145-154. [5] S. R. Lubkin, T. Jackson. Multiphase mechanics of capsule formation in tumors. Journal of Biomechanical Engineering, 2002, 124(2): 1-7. [6] X. He, M. Dembo. Numerical simulation of oil-droplet cleavage by surfactant. Journal of Biomechanical Engineering, 1996, 118 (2): 201-209. [7] M. Stastna. A moving boundary value problem in soft tissue mechanics. Journal of Canadian Applied Mathematics Quarterly, 2005, 13(2): 183-198. [8] Q. Jiang, Z. Li and S. R. Lubkin. Analysis and computation for a fluid mixture model. Communication in Computational Physics, 2009, 5: 620-634. [9] R. Li, Z. Chen and W. Wu. Generalized difference methods for differential equation: Numerical analysis of finite volume methods. New York: Marcel Dekker, 1999. [10] Z. Zhang. Error estimate of finite volume element method for the pollution in ground-water flow. Numerical Method for Partial Differential Equations, 2009, 25(2): 259-274. [11] 王平, 张志跃. 有限体积元数值方法在大气污染模式中的应 用[J]. 计算物理, 2009, 26(5): 656-664. [12] F. Gao, Y. Yuan. The characteristic finite volume element method for the nonlinear convection-dominated diffusion problem. C 能够很好的模拟单细胞变形且具有二阶精度，而 对于多细胞问题—间断问题，则会出现虚假的数值震 荡。广义迎风方法虽然对单细胞问题仅有一阶精度， 但它能够模拟多细胞问题，这是广义迎风方法的优 势。在实际应用中，我们可以根据具体情况选取适当 的方法来模拟组织细胞的变形。 5. 致谢 and Mathematics with Applications, 2008, 56(1): 71- 81. [13] E. Richard, L. Raytcho and Y. Lin. Fi approximations of nonlocal reactive flows in porous media. Numerical Method for Partial Differential Equations, 2000, 16(3): 285-311. [14] 张文生. 科学计算中的偏微 分方 程有限 差分法[M]. 北京: 高 等教育出版社, 2006. [15] E. P. Doolan, J. J. H. Miller and W. H. A. Schilders. Uniform numerical method Dublin: Boole, 1980. 本文感谢国家 省自然科学基 |