二维稳态晶体生长控制方程的数值解

D0I:10.13374/j.issn1001053x.2005.05.043 第27卷第5期 北京科技大学学报 VoL.27 No.5 2005年10月 Journal of University of Science and Technology Beijing 0ct.2005 二维稳态晶体生长控制方程的数值解 廖福成”陶娟”刘贺平2) 1)北京科技大学应用科学学院，北京1000832)北京科技大学信息工程学院，北京100083 摘要分析了在均匀流场的作用下，金属凝固过程中晶体生长浓度的二维稳态方程的边 值问题.运用有限差分法将徽分方程数值离散化为线性代数方程组.用初等变换法将该代数 方程组分解为多个方程组进行处理，提高了计算效率，模拟结果揭示了在均匀流场作用下， 沿枝晶生长的方向，晶体生长的浓度呈现振荡衰减的本质特征 关键词晶体生长：偏微分方程；边值问题：数值解 分类号0175 1问题的提出 性方程组，再求解该线性方程组即得到原问题的 数值解.计算结果表明，在固液界面前沿晶体生 金属凝固过程的一个重要的问题是形核后 长浓度呈现振荡衰减的本质特征. 的晶体生长.关于晶体生长的理论和实验研究一 直非常活跃.文献[5]在电化学沉积实验中发 2差分格式 现在扰动电流的作用下，枝晶尖端附近的浓度场 是振荡的.文献[6-9]从数学的角度求解晶体生长 对于边界条件(3)、(4)作如下的处理.由于C 稳定态的控制方程，从理论上证实了固液界面前 对x的周期性，这里仅讨论一个周期即0≤x≤2l内 沿浓度呈现指数振荡衰减性.晶体生长控制方程 的解，再将解延拓到整个x轴上即得所有的解， 大多是复杂的偏微分方程回.在一定条件下，一 注意到： 类二维稳态晶体生长浓度控制方程及其边界条 C(0,z=C(21,zg(z) (5) 件为： 由于limC(x,z)=0,所以对任意e>0,存在L,使z>L, 98+8-0 (1) 时，对x∈[0,2]一致地有0≤Cx,z<e.为作计算设： Cx,L)=0(0≤x≤20 C(x,0)=f(x) (6) (2) C(x+21,z)=C(x,z) (3) 把它作为一个边界条件，利用上述性质可用gz) limC(x,z)=0 (4) 在无穷远处的性质和实际允许误差界来选择L. 下面求问题(1)，(2)，(5)，(6)的数值解.将矩形 其中D为扩散系数，速度分量y,及周期2l均为 区域-{(x,z:0≤x≤2l,0sz≤L}分成N×M个小矩 常数，x)是关于x的连续周期函数，周期为2l. 形，长、宽分别为：△x=h,△z=h,h,和h亦称空间 文献[7,10]已经给出以上问题的Fourier级数 步，这里h和h的选择应按实际允许误差界来确 形式的解析解.但对于应用来说，有时数值解更 为方便，本文用差分法求该问题的数值解，从数 定，使+与该实际允许误差界至少同阶.记 值解的角度来说明晶体生长的振荡衰减性.用差 1=x+h1,x-=x-h(=1,2,…,N-1),z州1=z+h,z-= -h0-1,2,…,M-1),则有 分法解偏微分方程，常用的方法有古典显格式、 古典隐格式以及Crank-Nicholson格式i.m.针对 C.6k,zC12》-2Cg3HC.-22+0m (7) 二维问题提出了一个截断误差阶为O(+)的显 C.)-Cm)-2C)+Ch) (8) 式差分格式，将偏微分方程离散化，得到一个线 C.)-CC) 2h1 (9) 收稿日期：200410-27修回日期：200503-18 C.k,zC62l-C2+O份) 基金项目：国家重点基础规划项目No.G2000067206-1) 2h2 (10) 作者简介：廖福成(1957一)，男，教授 将式(7)(10)中的C(x简记为C,代入式(1)，将

第2 , 卷 第 5 期 2 005 年 1 0 月 北 京 科 技 大 学 学 报 OJ u r . a l o f U n iv esr iyt o f s e ieD c e a n d 及 c血n o l o gy B e ij 恤g V b L2 7 N O 一 5 O Ct . 20 5 二维稳态晶体生长控制方程的数值解 廖福 成 ” 陶 娟 ” 刘 贺平 ” l )北京 科技大 学应 用科 学学 院 , 北京 10 0 83 2 )北京 科技 大学 信息 工程 学院 , 北 京 10 0 83 摘 要 分 析 了在均匀 流 场的 作用 下 , 金 属凝 固 过程 中 晶体 生 长浓 度 的二维 稳态 方 程 的边 值 问题 . 运用 有 限差分法 将微 分 方程数 值离 散化 为线 性代 数方程 组 . 用 初等 变换 法将 该代 数 方程 组分解 为 多个 方程 组进 行处 理 , 提 高 了计 算效 率 . 模拟 结 果揭 示 了在均 匀流 场作 用 下 , 沿 枝 晶生长 的方 向 , 晶体 生长 的浓 度呈 现振 荡衰 减 的本质 特 征 关键 词 晶体 生长 ; 偏微 分方程 ; 边 值 问题 ; 数 值解 分类号 0 1 7 5 1 问题 的提 出 金 属 凝 固 过程 的一 个 重要 的 问题 是 形 核 后 的晶体 生长 . 关 于 晶体 生长 的理 论和 实验 研 究一 直 非常 活 跃 ’冈 . 文 献 【5] 在 电化 学 沉积 实验 中发 现在 扰 动 电流 的作用 下 , 枝 晶尖 端 附近 的浓度 场 是振 荡 的 . 文献 [ 6 一 9] 从数 学的角 度 求解 晶体 生长 稳定 态 的控制 方程 , 从理 论上 证 实 了固液 界 面前 沿浓度 呈 现指 数振 荡 衰减 性 . 晶 体生 长控 制方 程 大 多是 复杂 的偏 微 分 方程 `10J . 在 一 定条 件 下 , 一 类 二 维 稳 态 晶 体 生长 浓 度 控 制 方程 及 其 边 界 条 件 为 : 性方 程 组 , 再求解该 线 性方 程组 即得 到原 问题 的 数值 解 . 计 算 结 果表 明 , 在 固液界 面 前沿 晶体 生 长浓 度 呈现 振 荡 衰减 的本质 特征 . , ( 刁 Z C . 刁 ZC 、 . 刁C . 刁C 。 代石厂百了+) 、 、丁产认飞百= u 侧方 , 0) 气尹’(x ) O (x + 2劫卜C ( x 力 l m C ( x 习二 O 其 中 D 为扩 散系 数 , 速 度 分 量 vx , 铸及 周期 21 均 为 常 数 沪以)是 关 于x 的连 续周 期 函数 , 周 期 为21 . 文 献 【7 , 10] 己 经给 出 以上 问题 的 F o iur er 级 数 形 式 的解 析 解 . 但对 于 应用 来 说 , 有 时数 值 解 更 为方 便 . 本 文 用 差分 法 求该 问题 的数 值解 , 从数 值 解 的角度 来 说 明晶体 生 长的 振荡 衰减 性 . 用 差 分 法 解偏 微 分 方程 , 常用 的方 法有 古 典显 格 式 、 古典 隐格 式 以及 C r a n 卜N i e h o l s o n 格 式 `, , 一 , ,」 . 针对 二 维 问题 提 出 了一 个截 断 误 差阶 为 以斌+ 从) 的显 式差 分格 式 , 将 偏 微 分方 程 离散 化 , 得 到 一个 线 收稿 日期 : 2 0 -4 1仁27 修 回 日期 : 2 0 5刁3 一 18 基金 项 目 : 国家重 点基础 规划 项 目N( .o 2G o 0() 6 7 2 06 一 l) 作者 简介 : 廖福 成 ( 19 57 一) , 男 , 教授 2 差 分格 式 对 于边 界 条件 (3 ) 、 (4 )作 如下 的处理 . 由于 C 对 x 的周 期性 , 这 里 仅讨 论 一 个周 期 即 O` x ` 21 内 的解 , 再 将解 延 拓 到 整 个 x 轴 上 即 得所 有 的解 . 注 意 到 : (C 0声 ) = (C 2今)气抓幻 ( 5 ) 由于 h m (C 琳力二0 , 所 以对 任意夕O , 存 在 L , , 使 艺> 乙 : 时 , 对 x 任 【0, 2月一 致地 有 O` C x( 习 e< . 为 作计 算 设 : C (才声 1 ) = 0 ( 0` x ` 2乃 ( 6 ) 把它 作 为 一个 边 界 条件 , 利 用 上述 性 质 可用以宕) 在无 穷 远 处 的性质 和 实 际 允许 误 差 界来 选 择L , . 下面 求 问题 ( l ) , ( 2 ) , ( 5 ) , (6 ) 的数值 解 . 将矩 形 区 域月七 { x( , z) :0 ` x ` 2人0` z ` L : } 分 成 Nx M个 小 矩 形 , 长 、 宽分 别 为 : xA = h , , 山= 爪 , h : 和爪亦称 空 间 步 . 这 里 h l和 h Z的选 择应 按 实际允 许 误差 界来 确 定 , 使材+h 圣与 该 实 际 允 许 误差 界 至 少 同阶 . 记 x +i , = +x, h l , x `一 , = x `一 h : (=1 1 , 2 , … , N 一 l ) , ’+z : = 与+ 人 2 , 为 一 : = jz 一 h Z价1 , 2 , … ,材` l) , 则 有 。 x( 声卜巫业裸严巫城 区h `, `, , 。 x( , 卜巫扯嘿笋鱼咧 + 口`” ,, “ , (2)l()34 cx hx( 二匹陈罕回 十 、 。 hx( = 恤斋哑 十 、 (9 ) ( 10 ) 将 式 (7) 一 ( 10) 中的 C示 j沟 )简记 为 q , 代 入式 ( l) , 将 DOI: 10. 13374 /j . issn1001 -053x. 2005. 05. 043

Vol.27 No.5 廖福成等：二维稳态晶体生长控制方程的数值解 ·561· O()和O(h)去掉，即知近似地有： 户M-2时， D(CR-2Cu+C-uCm-2Cu+Cu-i aC3m-2+bCiM-+cC2M-+dC2M-+eCiM-2-0, 4C0, bxx-2+dx--ex-CXx-a0; 2h2 j=M-1时， aC3M-1+bCiM-1+cCiv+dCiM-2+eCiM-1-0, 0}0 (11) bxw-+dx-u-+exaaxsw-n=cf (2). 会或，层-六6,0%，是元 D D (3)=W-2. -2D2D j=1,aCx-1+bCx-3.+cCw-12+dCx-20+eCx-21-0, 所一房e,则式1)变为： bXwex-3+CX-3X-2c-2x aCm/+bC-u/+cCim+dCu-1+eCu=0 (12) -d(N-2): 下面处理边界条件.一般地考虑C(x,0)(x), j2,aCx-12+bCx-32+cCx-23+dCx-+eCx-22=0, Cx,L)=(x),C0,z=g(a),C(24,z)=g(a).在此， bxow-xm+d-xM-x-3x-+CX-3xM-+ (x)=0,g(广g(z).离散化边界条件：C.(x, QXON-2XM-120; C.M=x),Cw=g(,Cw=g().为了简化记号， 可使用f()和5（团分别代替fc)和5(xi=0,1,, j户M-2时， M:用g)和g20)分别代替g(z)和g(ai=0,1,,M. aCx-IM-2+bCx-1M-2+cCx-2M-1+dCw-2M-+eCx-i-1=0, bxcN-sXM-mw-2+dkW-3KM-lnM-ytexw-3Xw-w-2+ 3线性方程组的建立及矩阵表示 CX(N-2KM-1+aX-2XM--O; 注意这里要求的是C(i=1,2,N-11,2,, j户M-1时， M-1).若令 aCx-M-+bCx-3M-1+cCx-2w+dCN-2M-2+eCx-1M-1=0, Ci-XUM-IXi-1 (13) bxo-wdxo-xex-2-x-1x- 则要求的即为向量x(=1,2,…,N-1)M-1)》.这是 -cf(N-2); (W-1M-1)个量.为了求出它们，用一个 (4)=N-1. (N-1)(M-1)x(N-1)M-1)的方阵来存储系数.下 1时， acm+bCx-21+cCN-12+dCv-10+eCN-1=0, 面分步讨论： (1)-1. bxw-Xw-exw-2xX-2XM-2 户l时，aC2i+bCo+cCa+dCo+eCi=0,即 -dW-1)-ag(1); ex:+cxa+axco-m=-df (1)-bgi(1); 2,aCx2+bCx-22+cCw-13+dCx-1+eCN-12=0, 户2时，aC22+bCaa+cCis+dC.a+eCz=0,即 bxw-3xw-12+dxo-2-ex-2x- dx+ex:+cx;+ax-m2=-bgi(2); CX(-2XM-1m=-ag2(2); j=3时，aC2a+bC+cCi4+dC2+eCia=0,即 dx2+ex;+cxs+axu-i=-bg(3); =M-2时， aCxM-2+bCx-1M-:+cCx-1M-1+dCx-1M-+eCv-1M--O, j=M-2时， bxw-3xw-w-2+dxw-2xw--exw-2m--2+ aCiM-2+bCoM-1+cCM-+dCiM-3+eCM-2=0. cxw-IxM--ag2(M-2); dxw-3+exn-2+cxx-+axw--2=-bg(M-2); j=M-1时， =M-1时， aCxM-1+bCx-1M-1+cCx-Ix+dCy-1M-2+eCx-1M-1=O, aCiM-1+bCox-1+cCi+dCim-2+eCiM-1-0, bxw-2xu-i+dxo-xw-I-+exow-x--(N-1)- dxw-2+exx-+axx-=-bg(M-1)-cf (1). ag(M-1)片 (2)-2. 将以上结果写为矩阵的形式，得到： 户l时，aC,+bC.+cC22+dC2o+eC2=0,即 AX=B (14) 「Eal bx+exw-m+cxo-im+axx-m=-dfi(2); bI E. =2时，aCa+bCz+cC,+dC,+eC2a=0,即 其中，A= bxz+dxw-+ex-+Cx-axO; ..al bI EJo-1xw-D)-w-1xM-1)

V b L 2 7 N 0 . 5 脚福 成 等 : 二 维稳 态 晶体生 长控 制方 程 的数 值解 . 5 6 1 . 以 h圣)和 以 h孟)去掉 , 即 知近 似地 有 : 血篇 过鱼杆` ` 赞 兰丛习 + ) = 0 , 。p {争式 ~ } C+, I J 一刹 。 一 砂!令六! ( D vz 、 ~ 医 一 厄蔺)叭 一 , C +J t + ( 11 ) 2 畴十 认 令 景 .卜命=al ZD ZD 一石不一荡{一 “ , D 斌 苦升二 b . 乙 n l 一 刹 二=0 静命 = 。 , 一 D 认 ; , 厂 一 下污一二口 , 刀三 乙月2 则式 ( 11) 变 为 : a 瓜 , , + bC ` , +c 心 , 十动以, 一 , +e q = 0 ( 12 ) 下 面处 理边 界条 件 一般 地考 虑 C ( x , 0) 书 (x) , C x( , 乙) = 关x() , (C 0 , z) = gl (z) , (C l2, )z = 今 ()z . 在 此 , 关x( 卜0 , gl (z ) 二 9 2 (z ) . 离 散 化 边 界 条 件 : C , 。丫 x(,) , 口M = 关(x,) , 0C 。 = gl (z,) , 蛛 . , 二 乡(z, ) . 为了简 化记 号 , 可 使用厂i() 和关( i) 分 别 代 替厂x(,) 和关 x(,) i( 二 0 , 1 , … , 劝 ; 用gl 仍和承功分 别代替gl (z, )和今 ()jz 价 o , 1 , … , 动 . 3 线性 方 程 组的 建 立 及 矩 阵 表示 注 意 这 里 要 求 的 是 q (=1 1 , 2 , … 万一 1产 l , 2 , … , 材` 1) . 若 令 q =x( , . x ,一 ,脚 ( 13 ) 则要 求 的 即为 向量 xk (=k l , 2 , … , (N 一 1) (泪件 1) . 这 是 (刀一 l) (材` l) 个 量 . 为 了 求 出 它 们 , 用 一 个 N( 一 I XM` 1) 、 N( 一 1)(力卜 l) 的方 阵来存 储 系数 . 下 面 分步讨论 : ( l )户 1 . =j l 时 , a 2C . : + b OC , L+ e C t ; + 动C t . 。 + e C : , 1= 0 , 即 ex :+ cxz + 。工伽 一 ,卜 : = 一 晰 ( l ) 一 坛 1 ( l ) ; =j 2 时 , a Q , 2+ b CO ; + c 口 , 3 十J C I . t +e 口 2 = 0 , 即 么 ,+ xe Z+ cxs + 口义浏 - , 。 = 一 坛 : ( 2) ; =j 3 时 , a 二 J + bOC , + e 已 , + dC : ; + e C : . 3 = 0 , 即 改 + exs + ’cx + 口义脚 一 : 卜 3= 一 bg l (3 ) ; 少= 材` 2时 , a 二* 2+ bOC , 一 2+ c C : * 一 l十 d C t * 3+ e C :踌 2 = 0 . 即击 加 - 3+ xe 材 一 2+ cxu 一 ,+ 以附 一 , ~ 2 =) 一 坛 1 (初件 2) ; 产材` 1 时 , a Q 茫 : + b OC * l + c C ; 对+ J C I茫 2 + e C t沁 1= 0 , 即击桥 2+ 己尤拟 一 1+ axz `A-, : )= 一 坛 , (材` l )一抓 ( l ) . ( 2 ) 扮2 . =j l 时 , a 3C . : + b C I . ; + e Q ; + 动C全 , 。 + e Q , ; = 0 , 即 bx :+ ex( , t卜 t+ cx , 1 、 + 以 2泌 , 卜 l = 一亦 (2 ) ; =j 2 时 , a C 二+ b口召 + e 已 , 3+ 成几 1+ e Q 津= 0 , 即 加 2十么伽 ~ 、 卜 ;+ 已工《材 一 1 、 + 伪 , t卜 , + 口笼 2《 , : 卜 2 “ 0 ; 少= 泪件 2 时 , a G 、 2 + b C I * 2+ c G、 : + 动C Z * 3+ e Q 、 2= 0 , 即 旅 翻 - 2+ 么山 - 1 ~ 3)十以困 一 , 、 , 2汁cx Z( , , )十“ 戈2淤 : 洲* 2)= 0 : 少= 泪件 1 时 , a 已、 t+ b C ,踌 ,+ e 已万叼已价 2仰 Q 阶一0 , 即石义“ - : + dx( * l ) * 2》+ exz 、 :汁以 」淤 , )= 叭 ( 2 ) . ( 3) ￡二万一 2 . 产 l 时 , a G 一 1 , 1+ bG es 3 . : + e 蛛 一 2 , 2 + 故几 ` 2 。 + e G 一 2 . : = 0 , 即 西为刀一 ;胜 .卜 : + ex 。 一 3 、 1卜 1+ cx 。 一 3 , 1 、 + 口为柑一 2 , , ) + : = 一晰 (万一 2 ) : =j 2 时 , a G es l , 2+ b G 一 3; + e G 一: 3 + dG 一 2 , + e G 一; 2 = 0 , 即 阮卿一 ; 、 1卜 2+ 击价一 3 、 : 卜 : + 已两刀一 3。 一 t 卜 2+ 以切一 3 x , 1。 + “ 浑伊 - Zx桥 1况 = 0 ; 产初卜 2 时 , a G - : * 2+ bG 一 , * 2+ e G 一 2对一 侧G - 2 * , + e 岛 se Z * 2= 0 , 即 云为卜 ;版 t曰 , 2十) 去侧 一 3脸 .〕 * 3)+ ` 义留一 3x , : 洲* 2,+ ` x 训 一 Zx * t )+ ax( 、 一 Z x桥 : 州桥 2 )= 0 : 少= 初卜 1 时 , a G - , 、 1+ b C解 3踌 1+ c G 一 2 , + J G 一 2 、 2+ e G 一 2茫 : = 0 , 即 奴卜。 、 , ,十么卜。 , , 州衬 ` Zj+ 载 入-r 2任 :汁“ 从刀一 , 、 , ,= 一抓 (刃一 2 ) ; ( 4 ) 扮 N 一 1 . ’=j 1时 , a q : + bG 一 2 , , + e G 一 , , 2+ J G 一 1 . 。 + e G 一 , , 1 “ 0 , 即 加侧 一 3。 一 1卜 ,+ 己从~ 一 Zx , : 卜 : + 以卿 一 2 、 : 卜 2= 一亦 (N 一 l) 一 agz ( l) ; =j 2时 , a q Z+ bG 一 2; + e C 即 - : . 3十d G 一 : . : + e瓜 一 , 2= 0 , 即 西为刀一 3 * ,卜 2+ 六恻一 2 , * l卜 : 十以阴 一 Z x桥 ,卜 2+ cx 淤 , x -M 1* 3= 一 ag2 ( 2) : 了= 习件2 时 , a G 月 es Z+ b G 一 、 2+ c G 一 : 、 1+ J 〔认 一 , * 3+ e G - : * 2 = 0 , 即西为刀一 3脸 . ~ 2 )+ 么叨 一 2征 : 。 * 3)+ 以侧一 2、 , : 川 -M 2,+ c 尤例 - , 、 , =) 一 口9 2 (为卜2 ) ; 产初卜 1 时 , a G 洲 - , + bG 一 2井 1+ e G 一 1 , + 口瓜 一 : * 2+ e G 一 : * : = 0 , 即 石为柑一 2 、 ,汁么侧 一 Z x , : ~ 2,+ xe 。 一 L、 -M :厂一 抓 (N 一 l ) 一 。 9 2 (M` 1) ; 将 以上 结 果 写 为矩 阵 的形 式 , 得 到 : 月无卜丑 ( 14 ) E 口 I bI E ’ 二 bI E ] , 一 : 隆 : 〕 、 伽 一 1。 一 1)

·562。 北京科技大学学报 2005年第5期 X B, 「Eal B X B bI E al B X- ,B= B [AB]= bI E al Xw-2 Bw-2 Xx-1w-1XM-01 Bx-1N-1XM-1 bI E alBx-2 e c bI EB-1] d e c 对[AB]作初等变换，记E=E,B=B,从[AB] 这里E= 5 的第1行（这里的行指的是分块矩阵的行）开始， dec 把第1行左乘以(-bE)加到第2行，记E= d el-1rw-1) X-IXM-IHI E-ab(E)',B=B2-b(E)Bi,依此类推，把第k行 X-1X-12 左乘以(-bE)加到第k+1行，记E=E-ab(E-)', X= k=1,2,…,(N-1), B=B-b(E-)B-(-2,3,…,N-2),则[AB]化为： Xu-ixa-1(M-2) 「E,al B XRM-1 E al -0-g [AB]- E al 年 -bg(2) B= Ex-2 al Bv-2 -bg(M-2) Ex-B-1 -d5(I)-bg(M-1)小wm (E=E -(k)] 0 其中， Bi=B E=E-ab(E-1) (=2,3,…,N-1)(15) B ,k=2,3,…,(N-2), Bi=B-b(E)Bi 0 这样就得到与方程组(14)同解的方程组： cf (k)kM- (E X+aX;=Bi -dN-1)-ag(1) 7 EX:+aX;=B: -ag(2) (16) BN- Ew-XN-2+aXN-1=B-2 -ag(M-2) EN-X-B- -d5(N-1)-ag(M-1)w-i 最后，从式(16)的最后一组方程出发，解出Xw-, 解方程(14)即得所有(k=1,2,…,(N-1)× 采用回代方法即可求出X-X-,,X,从而得到 (M-1).再用式(13)求出C(i=1,2,,(W-1), 线性方程组(14)的解.再用式(13)求出Cw(i=1,2, 广=1,2，，(M-1),加上边界条件即得所有的C ,N-1,=1,2,M-1),加上边界条件即得所有 (i=0,1,,N,J=0,1,,0. 的Cu(i0,1,…N,j=0,1,…0. 计算时仅需用一个(N-1(M-1)×(M-I)阶矩 4线性方程组(14)的求解 阵压缩地存储系数矩阵[E,,环-J『即可，这样就 如果A可逆，则方程组(14)的解为X=AB.但 大大地节省了存储空间，用循环语句求解方程组 是这一结果并不便于计算，因为通常A是一个阶 (14).这样，当z很大时也能计算，且能大大提高 数很高的矩阵，这是因为为了使得到的数值解 计算速度.以上求解线性方程组的算法如下： C更好地近似代替原问题的解C(x,z,N一般取得 (I)为矩阵E,B1,B2,,Bw赋初值； 较大：limC(xz-0又要求对很大的z进行计算，从 (2)E+E 而边界条件(6)中的L,需尽可能地取大，从而M应 for i=2 to (N-1) 取得很大，这样A将会占用大量计算机内存，再 E-E1-ab(E-),B,-B,-b(E))B, 对A求逆又会占用大量内存，用X=AB求解方程 (3)Xx-(E-)Bw-1, 组会使我们无法求解、但注意到A是一个大型带 for i=(N-2)step-1 to 1 型稀疏矩阵，利用A的这一特点，就可以克服上 X(E)-(B-aX): 述困难.对方程组(14)的增广矩阵为： (4)方程组的解为X=[X灯，，…，X还J

5 6 2 . 北 京 科 技 大 学 学 报 2 0 0 5 年 第 5 期 嚎卜 卜 「 瓶琳 尤 1 「 B门 { 凡 … { B , { =X } _ : … , =B 「 ’ } {圣 一 」 }肠一 { 【恙 一 , ] 扭一 以-M : ) 、 , L肠 一 ; ] 扭一 : x , L) X , …气 “ { 示… “ :.f1t … … “ 气 “ L “ e 」…《M ~ 1) x 《刊一 一) }聋 * 一 { } ~ tk一 I XM 一 1旧 l A[ } B ] = E 口 I lb E al lb E al 云了 E bI 这里 尤 , = 卜 1 , 2 , … , 伽一 l) 对 A[ }B] 作初 等 变换 . 记乙=E , 斌=B , , 从 回 }B] 的第 1 行 (这 里 的行 指 的是 分块 矩 阵 的行 ) 开始 , 把 第 1 行 左 乘 以卜 bE 内加 到 第 2 行 , 记 及= E 一 ab 伍J 一 , 』卜及一 b (E l犷 ’ B ; , 依 此类 推 , 把 第 k 行 左 乘 以(一 石及 , )加 到第 粉 1行 , 记及燕 E 一 ab *(E 一 l犷 ’ , 及扭厂 b (kE 一 : ) 一 1斌 一 :卜2 , 3 , … , 万一 2) , 则冈 la ]化 为 : 国阮卜嘛 乱aI 、 . 乱 lx( ` 一 1 权 日 l洲材- 2 ) 弋桥 , ) 」、 : ) 、 : 一 哪 ( l ) 一 坛 : ( l ) 一 坛 . ( 2 ) A[ } B 卜一 E , aI 及 aI 瓦 aI 一 坛 : (材` 2 ) 一联 ( l万坛 。 (材` l )」, ,卜 : 一 晰 (k) } 其 中 B 户 卜2, 3 , … , (N 一 2 ) , }巴=E 1竺`绍 1 }召{ 二甲 一口 西《马 - tB 三二B 厂 b咬E 几一 t : ) 一 (=k 2 , 3 , 一 N 一 l ) ( 15) ) 一 ’ )B ; 一 t 0 一 叭 (k) 这样 就 得到 与 方程 组 ( 14) 同解 的 方程 组 : 一 岁(万一 l ) 一 agz ( l) 一口g 泛( 2 ) 抚 _ . = 一 口9 2 (M ` 2 ) 一 联 N( 一 1)一。 肠 (材` l) ( 16 ) : = 刀二 一 2 (材 一 1) x l 解 方 程 ( 14 ) 即 得 所 有 x ` (k = l , 2 , … , (N 一 l ) x (材` 1)) . 再 用 式 ( 13) 求 出 ilC , ( i = l , 2 , … , N( 一 l ) , j = l , 2 , … , (扣卜 l) ) , 加 上 边 界 条件 即得 所 有 的 q (’=l 0 , 1 , … , N, j = 0 , 1 , … , 功 . 4 线 性方 程 组 ( 1 4) 的 求解 如果A 可逆 , 则方程 组 ( 14 ) 的解 为 方华刁 一 BI . 但 是这 一结 果并 不便 于计 算 , 因为通 常 通是 一 个阶 数 很 高 的矩 阵 . 这 是 因为 为 了使得 到 的数 值 解 q 更 好地 近似 代 替原 问题 的解 以尤 .z) , N 一般 取 得 较 大 ; h m C。 力= 0 又 要 求对 很 大 的z 进行 计 算 , 从 而 边 界条 件 (6 )中的L : 需尽 可 能地 取 大 , 从 而对应 取 得很 大 . 这 样A 将 会 占用 大 量计 算机 内存 , 再 对 A 求逆 又会 占用 大 量 内存 , 用方片月 一 ’ B 求 解 方程 组会 使我 们无 法 求解 、 但 注意 到 A 是一个 大 型带 型 稀疏 矩 阵 . 利 用月 的这 一 特 点 , 就 可 以克 服上 述 困难 . 对方 程 组 ( 14) 的增广 矩 阵 为 : 最 后 , 从式 ( 16 ) 的最后 一 组 方程 出发 , 解 出恙 - : , 采用 回代 方 法 即可 求 出恙 ~ 2 ,恙 一 3 , … 犬 , 从 而 得 到 线性 方 程 组 ( 14) 的 解 . 再 用式 ( 13) 求 出 C,J (=1 1 , 2, … , N 一 1 =,j 1 , 2 , … , 材` l) , 加 上 边 界条 件 即得所 有 的 q (卜0 , 1 , …入 =j 0 , 1 , … , 劝 . 计 算 时仅 需 用 一 个 伽- 1) (泪件 1) x (M ` 1) 阶矩 阵压缩 地 存储 系数 矩 阵 [耳 , … , 虱 一 」 T 即可 , 这 样就 大 大地 节 省 了存储 空间 , 用循 环 语句 求解 方程 组 ( 14 ) . 这样 , 当 z 很 大 时也 能计算 , 且能 大 大提 高 计 算速 度 . 以上求 解 线性 方 程 组 的算 法 如下 : ( l ) 为矩 阵 E , B , , 及 , … , 几 一 : 赋初 值 ; ( 2 ) E : ~ E fo r 扮Z ot (N 一 1) 忿 一 E l 一 ab (赚 : ) 一 ’ , B `一 B : 一 b田犷 , )B 卜 : , ( 3风 一 1= (乙 - , ) 一 ,肠 一 : , fo r 卜N( 一 2 ) s t eP 一 1 t o l 戈= (E,) 一 1 (B `一 口琴 十 J ; (4 )方程组 的解 为方阵男 , 盯 , … , 蛛 一 : J T

VoL27 No.5 廖福成等：二维稳态晶体生长控制方程的数值解 ·563· 5 仿真实例 1.0 在式(1)，(2)，(5)，(6)中取D=1,,=1,=1,N=20, 0.5 L=10及L=100,相应地取M=30及仁300，Cx,0)= 0 cos(x),C(0,z)=C(21,z)=exp(-z/2)cos(3z).Mat- 0.5 lab语言，编出上节算法的程序.得到的模拟结果 如图1和图2.从图中可以看出，在固液界面前沿 50 晶体生长浓度呈现振荡衰减的.这与已有的实验 00 2 结果是一致的.图3是在L=100时，沿z方向晶 图2晶体生长浓度C的分布（化，=100） 体生长的浓度C分布剖面图， Fig.2 Distribution of mass concentration for crystal growth (L,=100) 0.6 1.0 0.4 0.5、 0.2 0↓ 0 0.5 ·00 0.2 -1 0.4 10 0102030405060708090100 00 2 图3沿z方向晶体生长浓度分布的剖面图化=100) 图1晶体生长浓度C的分布(L=0) Fig.3 Section diagram of mass concentration distribution for Fig.1 Distribution of mass concentration for crystal growth(L,=10) crystal growth along the direction (L,=100 参考文献 制方程及解析解.中国有色金属报，2002,12（专集1少57 [)廖福成，郑连存.稳态品体生长控制方程的解析解.见：全 [1]Nash G E,Glicksman M E.Capillarity-limited steady state de- 国信息与计算科学学术研讨会论文集，西安：陕西人民出 ndritic growth(I):theoretical development.Acta Metall,1974 版社，2002.31 22:1283 [⑧)廖福成，祖翠娥，郑连存，等，控制熔体浓度三维稳定态方 [2]Xu J J.Interfacial wave theory of solidification:dentritic pattern 程的精确解.北京科技大学学报，2004,26(153 formation and selection of tip velocity.Phys Rev A,1991,15(43): [9]王自东，周永利，常国威等，控制单向合金凝固界面形态 930 [3]Xu JJ.Generalized needle solution,interfacial instabilities and 非线性力学方程.中国科学(E辑)，1999,29(1)：】 [10]王凤英，陈明文，孙仁济，等.三维稳态晶体生长的物理本 pattern formation.Phys Rev E,1996,53(5):5051 质.北京科技大学学报，2003,25(3)：230 [4]Xu J J.Interfacial Wave Theory of Pattern Formation.Berlin [11】徐长发，李红.偏微分方程数值解法，武汉：华中理工大学 Springer-Verlag,1998 出版社，2000 [5]Wang M,Zhong S,Yin X B.Nanostructured copper filaments in [I2]Mathews JH,Fink K D著.数值方法(MATLAB版)陈渝等 electrochemical deposition.Phys Rev Lett,2001,86:3827 译.北京：电子工业出版社，2002 [)]孙仁济，王自东，陈明文，等.铝合金晶体生长稳态时的控 Numerical solution of governing equations for two-dimension steady state crystal growth LIAO Fucheng,TAO Juan,LIU Heping 1)Applied Science School,University of Science and Technology Beijing,Beijing 100083,China 2)Information Engineering School,University of Science and Technology Beijing,Beijing 100083,China ABSTRCT A boundary value problem of governing equations for the concentration of crystal growth is solved in the two-dimension steady state considering the effect of uniform convection field.The differential equation is nu- merically discretized into a system of linear algebraic equations by using the finite difference method.In order to improve computational efficiency,the system of linear algebraic equations is decomposed to several sub-systems. The result of numerical simulation shows that the concentration of crystal growth in steady state presents oscillating attenuation along the direction of dendrite growth in the action of uniform convection field. KEY WORDS crystal growth;partial differential equation;boundary value problem;numerical solution

Vb L2 7 N 0 . 5 廖福 成 等 : 二 维稳 态 晶体生 长控 制方 程 的数值 解 一 5 6 3 - 5 仿真 实例 在式 ( l ) , (2 ) , ( 5 ) , (6 ) 中取 =D l , 护1 , 护1 , =N 20 , L l = 10 及 L t = 100 , 相 应 地 取人卜 3 0 及几卜3 0 0 , 侧沐 , 0 ) = e o s (x ) , (C 0习= (C 2切 = e x P卜z2/ ) e o s ( 32) . 采 用 M at - lab 语 言 , 编 出上节 算法 的程序 . 得 到 的模拟 结 果 如 图 1 和 图 2 . 从 图中 可 以看 出 , 在 固液 界 面前 沿 晶体 生长 浓度 呈现 振荡 衰减 的 . 这 与 已 有 的实验 结果 是 一致 的 `5] . 图 3 是在乙二 10 0时 , 沿 z 方 向晶 体 生 长 的浓度 C 分 布剖 面 图 . 0 . 5 Q 0 刁 . 5 一 、a u U 图 2 晶体 生长 浓度 C 的 分布 (L . = 1 0 ) iF g . 2 D is t ir b u it o . o f m a s , e o . e e n t r a it o n fo r e yr s at l g r o 叭沈h ( L , = 10 0 ) 0 . 6 0 . 4 0 . 2 n ù`, 欲 ǎ的z81 . à口0 刁 . 4 0岭1 0 2 0 一3 0 4 0 5 一0 6 0 7 0一8 0 9 0 1 0 0 图 l 晶体 生长 浓度 C 的分布 ( L . = 0) F ig . I D 七t r i b u it o . o f m朋 5 c o n e e n t ar iot n fo r e 叮s扭 1 g or 袱h (L . = 10) 参 考 文 献 [ l 」N a s h G E , G li e ks m an M E . C iaP l iar ty 一 lim i te d s t e a d y s at e d e - 口 d ir it e gID w ht ( 1) : ht e o re it c a l d e ve l o四 e nt . A e 加 M e at .l 19 7 4 , 2 2 : 12 8 3 [ 2 ] Xu J J . nI te arf e i a l w va e ht e o ry o f s o lid i if c a t i o n : d e n itr t i e P a t e rn of rm at i o n an d s e l e e it on o f it P v e loc i .ty P b y s R e v A , 19 9 1 , 15 ( 4 3 ) : 9 3 0 [ 3 ] X u J J . eG n e r a lize d ne e id e s o lut ion , in te arf e ial in s atb ili t i e s an d P a t e m of rm a t lon . P b ys R e v E , 19 9 6 , 5 3 ( 5) : 5 0 5 1 4[ ] Xu J J . iht e r af c lal 叭a/ v e T七e o ry o f P at em F o arm it o n . B e r li n S P n o g e r ` Ve r l a g , 19 9 8 [ 5 ] w 山l g M , z b o n g s , Y in x B . N an o ’strU e姗 d e 叩衅 if l am e n t s i n e l e e加hc em i c a l de P o s iit on . P b y s R ve L e t , 2 0 0 1 , 86 : 3 82 7 6[] 孙仁济 , 王 自东 , 陈明文 , 等 . 铝合 金晶体 生 长 稳态 时的控 Z 图 3 沿z 方 向晶 体 生长 浓度 分布 的剖面 图 (L I = 10 ) F i咨3 S e e 6 o n d i a g ar m o f m a s s e o . e e n t r a 找o n 由 , itr b u it o n fo r e yr s t a l g门 w ht a l o n g th e d i溉it o n (L . , 10 ) 制方 程及 解析 解 . 中国 有色金 属 报 , 2 0 02 , 12( 专集 l) : 57 7[ 〕 廖福 成 , 郑连存 . 稳 态晶 体生长 控制 方程 的解析 解 . 见 : 全 国信 息与 计算科 学学 术研讨会 论文集 , 西 安 : 陕 西人 民出 版社 , 2 0 0 2 . 3 1 18 ] 廖福 成 , 祖 翠娥 , 郑连存 , 等 . 控制 熔体 浓度三 维稳 定态 方 程的 精确解 . 北京科 技大 学学报 , 2 0 0 4 , 2 6( l ) : 5 3 19 1 王 自东 , 周 永利 , 常 国威等 . 控制 单 向合金凝 固 界面形 态 非线 性力 学方程 . 中国 科 学( E 辑 ) , 1 99 9 , 2 9 ( I ) : ] 【10 ] 王 凤 英 , 陈 明文 , 孙仁 济 , 等 . 三维 稳态 晶体 生长 的物理 本 质 . 北京科 技大 学学报 , 2 0 0 3 , 2 5(3 ) : 2 3 0 【川 徐 长发 , 李 红 . 偏 微分 方程 数值解 法 . 武汉 :华 中理工 大学 出版社 , 2 0 0 0 [ 12 ] M a ht e w s J H , r ink K D 著 . 数值方 法 ( M A T L A B 版 )陈渝等 译 . 北京 : 电子工 业 出版 社 , 2 0 02 N u m ier e a l s o l u t i o n o f g o v em i n g e q u a t i o n s ofr tw o 一 d im e n s i o n s t e a d y s at e c ry s at l g or w th L IA O uF e h e叮 , , , TA O uJ a n , , , IL U H op i n犷 ) l ) A P P li e d S e ie cn e S e h o o l , U n i v e rs ity o f s c i e n e e an d 介 e ho 0 1 o gy B e ij i n g , B e ij in g l 0 00 8 3 , C h i n a 2 ) ih of n 刀 a t lon E n ig ne n n g S e ho o L U n i V ers ity o f s e i en e e an d 介 e h n 0 1 o gy B e ij 运g , B e ij ign 1 00 0 83 , C ih an A B S T R C T A b o un d a ry v a l u e P or b l e m o f g o v em in g e q u iat o n s of r ht e e o n e e n tr at i o n o f e ry s at l 脚袱h 1 5 s o vl e d in ht e tw o 一 d im e n s i o n s te a d y s at e e o ns ide ir n g ht e e fe e t o f u in of rm e o n v e e t i o n if e l d . T h e d i fe re n t i a l e qu iat o n i s nu - me ir e al ly d i s e re itz e d i n t o a s y s t e m o f line ar a lg e bar i e e q u iat o n s b y u s i n g ht e if n i t e d i fe o n e e m e ht o d . I n o 记e r ot im Por v e e o m P u at t i o an l e if e i e n e又 ht e s y s t e m o f li n e ar a lg e b r a i C e q u iat o n s 1 5 de c o m P o s e d t o s e v er a l s u b 一 s y s te m s · hT e 把s u lt o f n um e ir e a l s imu lat i o n s h ow s ht at ht e e o n e en tr a t i o n o f e ry s at l g or wt h i n s t e a d y s t ate P re s e n t s o s e il l a t in g at e n u a t i o n al o gn hte d ire e t i o n o f d e n d ir t e g ID wt h i n ht e a e t i o n o f un i of rm e o n v e e t i on if e ld · K E Y W O R D S c ry s at l g or wt h ; P a rt i a l d i fe 化n t i a l e qu at i o n : b o un d a ry v a l u e Por b l e m : nu m e ir e a l s o l u ti o n