固化界面失穩及型態發展之非線性解析及可視化研究(2/3)

行政院國家科學委員會專題研究計畫 期中進度報告

固化界面失穩及型態發展之非線性解析及可視化研究(2/3)

計畫類別： 個別型計畫 計畫編號： NSC92-2214-E-002-014- 執行期間： 92 年 08 月 01 日至 93 年 07 月 31 日 執行單位： 國立臺灣大學化學工程學系暨研究所 計畫主持人： 藍崇文 報告類型： 精簡報告 報告附件： 出席國際會議研究心得報告及發表論文 處理方式： 本計畫可公開查詢 中 華 民 國 93 年 5 月 28 日

國科會期中報告書

合金固化界面失穩及界面發展之可視化研究

執行單位：國立台灣大學化工系

計劃編號：NSC92-2214-E-002-014

執行時間：92 年 8 月 1 號至 93 年 7 月 31 號

計劃主持人：藍崇文 教授

目錄

中文摘要

1

英文摘要

2

第一章 緒論

3

第二章 實驗及模擬方法

6

第三章 結果與討論

9

3-1 平整界面

9

3-2 λ

c↔

λ

c

/2 細胞狀界面的轉換

1 2

3-3 Vλ

2

~常數的淺細胞界面

1 5

3-4 Multiplets 界面

1 8

第四章 結論

2 2

參考文獻 2 4

中文摘要

在本篇報告中，我們研究了一個 SCN/Acetone 合金在薄膜凝固過程 中，位於 Mullins-Sekerka 的不穩定理論所預測的臨界點附近的淺細胞狀 界面行為。我們進行了定量相場模擬以及極長時間尺度的實驗，兩者都是 從未被報導過的，我們的模擬結果在平整界面以及λc↔λc/2 區域與古典理 論呈現良好的一致性，在此λc是臨界波長，在拉速較高的區域，我們引入 了反包覆通量使我們的相場模擬具有定量性，而在此 Vλ2 ~常數的關係也 清楚顯現。最後，當波長接近λc/4，我們觀察到了 Multiplet 界面。就實驗 而言，我們進行了一個定速實驗來觀察界面的演化情形，因為溶質的慢慢 累積，以上模擬所描述的行為我們也在實驗上觀察到。另外一些有趣的非 線性動態行為，例如界面分裂、消除、溶質球吞沒、以及非對秤細胞狀界 面生長，也同時在模擬以及實驗上被詳細的紀錄下來。

Abstract

The morphological transitions of shallow cells near the onset of Mullins-Sekerka instability during the thin-film directional solidification of a succinonitrile/acetone alloy have been studied by both phase field simulation and long-time scale experiments. The simulation results show good agreement with classic theories at the planar and λc↔λc/2 cellular bifurcation regions; λc

is the critical wavelength. At a higher pulling speed V, the anti-trapping current is introduced for keeping local equilibrium, and the region for

Vλ2~constant is found. Finally, multiplets appear when λ is close to λc/4. In

the experiments, due to slow solute accumulation in the molten zone at a constant pulling speed, these morphological transitions were also revealed. Some interesting nonlinear phenomena, such as tip splitting, coarsening, pinch-off, and asymmetric cells, are also confirmed by both simulation and experiments.

第一章 緒論

在薄膜凝固程序中，界面在不穩定臨界點附近的行為，一直是合金凝 固領域中最重要的課題之一。在固定拉速或溫度梯度之溫場之下，溫度梯 度或拉速在超過某個臨界值之後，一個平整界面會產生不穩定的情形以至 於產生波紋狀的界面，古典的穩定性分析—MS 理論即完整描述這個行為 [1]，這種界面不穩定的現象來自於組成過冷克服了溫度梯度和固液之間 表面張力的穩定因素，進而發生界面崩潰的現象。從 1980 年代至今，許 許多多的實驗都觀察到：隨著組成過冷的大小，見面會從平整到淺細胞狀

（Shallow cells）到深細胞狀（Deep cells）最後到樹枝狀（Dendrites）的

界面，在這些形態改變的動態中，可以發現很多有趣的非線性行為，非常 值得探討。然而，因為在臨界點附近的拉速非常低，遠端濃度會隨時間改 變，且達到穩定型態的時間極長，故雖然在過去二十年，有許多理論[1-2]、 實驗[3-5]、以及模擬[6-8]的結果，而一個極長時間、空間尺度的研究，還 未被報導過。 在一個給定的溫度梯度下，隨著拉速的改變，我們可由 MS 理論知 道臨界拉速(Vc)和臨界波長(λc) ，然而，因為在臨界值附近的變化行為過 於劇烈，我們必須精準控這這些控制參數，不只是要非常長的熔區(等級 約為 DL/V，DL 是液相的擴散係數) ，試樣的純度也要極高，且拉速的控

制要極準確。根據 1993 年發表的 WL 理論[2]，界面位置大約要幾個 DL/kV 2 特性時間才能達到穩定，在此 k 是合金的分凝係數。因此，要探討這區域 的界面行為，長時間和長熔區的實驗還未被報導過。 因此，模擬便成了研究這個領域的一個重要手段，MIT 的 Brown 團 隊在 1985~1993 年間利用傳統的 Front-tracking 法來研究，並得到了很多 珍貴的結果，他們指出淺細胞狀的波長變化是由λc→λc/2→λc/4 [8]，且隨 著他們模擬能力的進步，他們也作了一些動態的結果[7]。然而，雖然實 驗上從未觀察過在λc/2 以下的波長，但是 Eshelman 等人[3]在實驗上觀察 淺細胞狀的波長變化時，他得出 Vλ2_{大約是一個常數的結論，且在由平面} 變成淺細胞狀界面的過程中，其行為為次臨界(subcritical)，與 Brown 等人 的模擬有很大的不同，其原因在於 Front-tracking 因為需要很長的時間在 作界面疊代，所以有很大的限制。除此之外，Front-tracking 法無法模擬界 面出現和消失的行為，意即它不能完全模擬實驗所看到的現象，故，想要 利用此方法模擬大規模區域的型態變化，仍然是不可能的任務。 相場模擬(Phase-field simulation)的出現[9-15]為這樣的問題帶來一絲 曙光，然而，相場模擬是最近十年才興起的模擬方法，想要相 Front-tracking 一樣可以定量化仍然相當困難，因為它有許多可調整的參數，會影響到模 擬的結果，其中界面厚度是一個關鍵值，為了要維持其熱力學守衡性，界 面厚度值要很小（~0.1µm），要模擬大規模時間和長度尺度（~1000µm） 的問題，幾乎不可能，因為界面厚度的限制，以往應用在薄膜凝固的相場

模擬的時間尺度大約都在µs 左右，遠較實驗的觀察尺度小，這些原因使

得與實驗一對一的相場模擬非常困難，尤其是因為界面過厚所導致的非平

衡現象—溶質包覆(Solute trapping)更是造成無法定量模擬的一大困難，最

近，Karma [16] 提出了反包覆通量(ATC)來消除溶質包覆效應，Lan 和 Shih

把這個方法加入好用的 WBM model 中，來模擬樹枝狀晶體生長，並取得 不錯的成果。 在本篇報告中，我們利用 WBM model 並在高速下加入 ATC 校正項， 來與一個極長時間尺度的時樣相比較，模擬結果與古典 MS 理論相比，取 得了突破性的進展，而我們也進行了一個利用超長時間尺度來完成觀察的 實驗，因為少部分的溶質累積，平整界面會慢慢崩潰成細胞狀界面，我們 的實驗和模擬結果也在型態上相符，雖然因為相場模擬的界面厚度限制， 我們還不能完全模擬真實系統，但這些結果還是讓相場模擬注入了一記強 心針，也在定性上印證了我們的實驗結果。

第二章 實驗及模擬方法

我們的相場模擬是利用由非平衡熱力學所推導出來的 WBM model， 而 ATC 項則在高速時使用，為了把主導方程式無因次化，我們所有的變 數都由一些特定物理量重新定量，且模擬區域跟著拉速 V 移動，其中濃 度以參雜濃度 c0來無因次化；時間則是

δ

2/DL；速度則是 DL/

δ

。如此，主 導方程式可整理如下：

(

2

)

* * 2 * * * * M S x y y x M y V t φ β β φ φ ηη φ ηη φ η ε φ φ ₋ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ − ∇ ⋅ ∇ = ∂ ∂ − ∂ ∂ ∗ ₍₁₎

{

[

(

1

)(

)

]

}

. * * * * * 0 * * * * * * a j + ∇ − − + ∇ ⋅ ∇ = ∂ ∂ − ∂ ∂ _φ A B S S c c c c D y c V t c (2) 我們模擬及實驗所考慮的系統是透明的 SCN/Acetone 合金，在 WBM

model 中所使用的物理常數都以純的 SCN 或 Acetone 來設定，除了 Acetone 的凝固潛熱除外，為了要配合 WBM model 中所使用的理想熔液假設，使 整個相圖符合真實情況，這個值要設成 9.651x107 J/m3。其他參數的推導 和設定，可參考其他文獻[17,18]。 唯一與傳統 WBM model 不同的地方在於高速時我們會引入反包覆通 量 ATC：ja，在移動座標下，他可以定義成： φ φ φ φ φ ∇ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ − ∂ ∂ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − + − = * * * ) ( ) 1 ( 1 * 2 ) 1 ( y V t h k k c k a * a j . (4) 在我們的模擬中，參數a=1/ 2 ，以使得模擬所得出來的分凝係數小於

0.14，有關於這個參數的討論，可參考其他的文獻[17]。 我們模擬所考慮的濃度是 0.01mol%，而溫度梯度是 100K/cm，而所 考慮的計算區域是 200µm(W) x 300µm(H)，如圖一所示： 圖一、相場模擬所使用的計算區域 我們並且忽略潛熱在凝固過程中所釋放的效應(FTA)，即溫度梯度 G 為一定值，在模擬上，此系統可以看成一個連續的熔液供料和合金出料， 在兩側，我們使用對稱的邊界條件(∂φ/∂x* =0,∂c*/∂x* =0)；在上邊界，則 使用固定濃度條件(c* =1)。為了要求解此方程式，我們使用適應性有限體 積法[14]來有效率的解決我們的問題。 實驗上，我們建造了一個大長度尺度的薄膜合金凝固系統，如同 Lee 等人的一樣[8]，其示意圖如圖二所示：

圖二、大長度尺度的薄膜凝固實驗示意圖 凝固的材料被封閉在 70µm 厚的兩片 Pyrex 玻璃之間，玻璃四周都以 Teflon 墊片塞住，銅製的冷卻和加熱系統可以製造差異大的高低溫，詳細 的溫場則用 K type 的熱電偶溫度計來測量。SCN/acetone 合金熔液則由外 注入此兩片玻璃之間，為了要製備極高純度的 SCN，我們使用真空蒸餾 法以及超過 60 次的精煉(Zone refining)過程純化。此外，為了要使生長方 向與熱流方向相等，我們使用快速凝固及冷卻的方法來減輕此效應，並且 靜置兩天以上，以確保裡面的容易可以均勻混合，然而因為參雜的控制很 難精準控制，我們使用 0.151mol%的試樣來進行實驗，此濃度是以崩潰速 度及 MS 理論反求而得。

第三章 結果與討論

3-1 平整界面 在 MS 理論中，在固定溫度梯度之下，當拉速超過某個臨界值，一個 平整介面會變得不穩定且開始呈現波紋狀，此波紋的波長我們稱之為λ。 在我們所呈現的模擬中(c0=0.01mol% and G=100K/cm)，利用 MS 理論所求 得的 Vc是 72.6µm/s，而臨界波長是 68.4µm。 圖三、模擬所求得的平整合金凝固的固相及液相濃度

因此，在測試不穩定的波長之前，我們應該先測試一個平整介面，看 看他的熱力學平衡性質如何，因此，我們首先模擬一個低於臨界拉速的凝 固問題(50 µm/s)，且給予適當的初始條件(液相濃度 cL*=1；固相濃度 cS*=0.1；界面溫度 Ti=331.2108K)。因為此拉速較低，非平衡熱力學校應 較不顯著，故我們在這個部分並沒有引入 ATC 來校正，當凝固開始進行， 如同基本冶金學所預測的，界面開始後退，溶質也隨之噴出，在此過程中， 液相濃度和固相濃度都被記錄起來以驗證平衡性質。如圖三所示，非平衡 效應很小以至於得出來的分凝係數和真實系統的值相當接近，除此之外， 濃度在液相的界面前沿的分佈是一指數衰減函數，與理論解析解相比有很 好的一致性。 為了要解析介面的動態行為，Warren 和 Langer[2]以濃度邊界層的特 性厚度為一 quasi steady-state 為假設，推倒出來一組界面位置和濃度邊界 層的特性長度偶合的 ODE，來改善傳統的 MS 理論，Losert 等人[19]也在 實驗證明了這個理論，在我們的模擬中，我們所求得的介面位置和特性厚 度也和 WL 理論呈現良好的一致性，如圖四(a)及(b)所示。在此圖中，z0 是界面在移動座標下的界面位置；lc是特性擴散厚度，很明顯的，在圖四 (b)中，界面位置需要 2DL/kV2 (10.2 秒)來達到穩態，換句話說，對於一般 濃度較高的實驗，其拉速約為 1µm/s，穩定的界面位置至少需要四五十個 小時才能發展完成。

圖四、相場模擬所求得的特性擴散厚度(a)及界面位置(b)之動態 為了避免溶質累積所造成實驗觀察上的不準確，越低濃度以及越大熔 區的系統是很重要的，以相場模式的發展而言，我們完成相場模擬在極低 濃度的模擬，算是踏出很重要的一步，以上的結果也告訴我們，當溶質包 覆效應很小時，模擬求得的結果可與理論定量化，這也讓我們對之後的結 果充滿信心。

3-2 λc↔λc/2 細胞狀界面的轉換 在相場模擬中，當拉速增加至 90µm/s，界面在 10DL/kV2 (15.68s)的時 間後仍然是維持平整，然而，當拉速再增至 100µm/s 時，在 5DL/kV 2 (6.3s) 的時間後此平整界面開始產生不穩定，如圖五所示： 圖五、模擬所求得的由平整界面到淺細胞狀界面演化圖 如圖五所示，在平整界面的某處會出現一個突起物，而界面在此突起 物的附近會開始產生不穩定，而不穩定就像一個波慢慢散佈到其他地方，

最後又再花了 5DL/kV 2_{才達到穩定狀態，而穩定狀態的波長是 33.3µm，根} 據 Brown 的理論，它位於λc/2 解的分支上(λc=68.4 µm)，這一系列的發展 和實驗上所觀察到的現象極為類似[20]。 我們逐漸增加拉速到 160µm/s 而其後再降低到 75µm/s，以檢驗其型 態的變化，我們延續 Brown 的研究，把所求得的細胞狀深度和拉速作圖， 其結果如圖六所示，模擬所求得的型態也同時畫於此圖上： 圖六、臨界點附近的波長分歧圖 在圖六中，我們可議看出拉速對細胞深度的關係，令人振奮的，我們

在低速時得到了與 MS 理論預測相符的λc波長，與 Brown 的 Front-tracking 模擬類似，在淺細胞界面由 85 (λc/2 mode) 到 80µm/s (λc mode)時，深度 會突然上升，而平整界面又會在拉速降至 75µm/s 得到，這與由 MS 理論 所預測的臨界拉速值 72.6 µm/s 相當接近，另外此圖也清楚的顯示了的由 平面到淺細胞狀的界面轉換是次臨界的，這個次臨界性質也由理論[21]及 實驗[4]所描述的類似。 圖七、溶質球的吞沒以及深度隨時間的震盪變化現象

除此之外，在我們的模擬中，深度的變化是隨時間穩定波動的，如同 Tsiveriotis 和 Brown[7] 的 結 果 類 似 ， 然 而 ， 因 為 他 們 的 結 果 是 利 用 Front-tracking 來進行模擬，所以無法觀察到由 Kurowski 等人[5]所觀察到 的溶質球吞沒的行為(pinch-off)，如他們在實驗上所觀察到的，在底部的 溶質會隨著深度的增加而被吞沒，如圖七所示。如前所述，到拉速為 160µm/s 為止，所求得的分凝係數仍在可接受的範圍內，故 ATC 並沒有 使用在這個區域。 圖八、實驗上所觀察到的λc/2 波長

另外，令人相當感興趣地，我們在我們的實驗中看到了λc/2 的波長， 到目前為止有報導過的文獻中，這是第一次看到，也映證我們的實驗系統 相當成功，如圖八所示。然而，因為這個分枝所存在的區域極小，所以他 只存在了一個小時就轉換到其他的波長，我們認為，在這一小時之中，因 熔區有限所造成的濃度變化不會超過 0.5%，可見此區域極難觀察。因此， 就算我們擁有極長的熔區，要把拉速調整在這個變動範圍仍然相當不容 易，我們使用一個在臨界拉速附近的固定拉速系統提供了一個很好的方法 來觀察這部分的行為，λc的波長也有可能利用稍微降速來加以求得。 3-3 Vλ2~常數的淺細胞界面 由 Tsiveriotis 和 Brown 所完成的分歧圖[7]顯示，λc/4 的分枝會出現在 λc/2 之後，與λc到λc/2 之間的行為類似，然而，λc/2 到λc/4 之間的波長還 是有在實驗上被觀察到[3]，並且 Vλ2_{大約是一個常數，我們為了要研究這} 個區域的行為，我們在模擬上把拉速更加提升。 在我們把拉速由 160 提升到 180 µm/s 之後，我們發現波長不再是 λc/2，而變成較小的波長，這時整個計算區域發生大規模的波長分裂、消 除，以及溶質球吞沒的非線性現象，這些現象的時間尺度仍然非常長，詳 細的非線性動態行為可以在我們出版的另一篇文章中看到[9]，然而，在 這個區域，一直到拉速為 750µm/s，因為拉速過快，我們觀察到很明顯的 溶質包覆現象，為了消除這個效應，我們使用 ATC 來消除這個非平衡效

應，波長與拉速之間的關係圖可參考圖九，許多得出的界面型態也畫於此 圖中，且我們的結果顯現出相當明顯的 Vλ2 ~常數關係，其偏差是因為我 們取的 domain 大小的緣故，如之前 Bi 和 Sekerka 探討的一樣[12]。 圖九、模擬所觀察到的 Vλ2 ~常數關係以及界面型態 在我們的實驗中，我們也觀察到了界於λc/2(255µm)和λc/4 (128µm)之 間的波長，如圖十所示，我們展示了平均波長為 206、194、184、171、 以及 160 的界面型態，與之前所觀察到的λc/2 分枝不同，它的維持時間相

長，約為十五個小時，因此，在一般的小長度尺度的實驗，這個關係是很 有可能可以被觀察到的[3-5,20]。簡單來說，我們的模擬和實驗支持淺細 胞狀在λc/2 分枝之後的行為不是λc/4，且 Vλ 2 ~常數關係似乎很符合此區域 界面轉換的行為。 圖十、實驗上所觀察到的較小波長的淺細胞狀界面 3-4 Multiplets 界面 Jamgotchian 等人[22]首次在即將淺細胞轉換為深細胞界面的附近發 現一排 Doublet 陣列，他們指出這是細胞壯界面的一個特殊分枝，且存在 區域非常狹小，Kopczynski 等人[23]更利用 Front-tracking 模擬得出其他分 枝，例如 Triplet 和 Quadruplet，在他們的分析中，這些都是解的分枝，並

不會同時出現，且這是由淺細胞到深細胞界面的一個過渡解。 圖十一、模擬和實驗所看到的 Multiplets 界面 在我們的相場模擬中，我們在拉速提升到 1000µm/s 時也看到了 Multiplets 界面，但是，我們看到的 Multiplets 會有很多分枝同時出現，如 圖十一所示，實驗所觀察到的行為也非常類似，所拍攝到的 Doublets、 Triplets，以及 Quadruplets 也在圖十一中與模擬結果比較，我們可以發現， 這些 Multiplet 群也在實驗中是同時出現的。 另外，由小波長的淺細胞界面到 Multiplets 之間的轉換是很有趣的，

並且從來沒有被實驗及模擬報導過，是故，在圖十二我們描述了拉速由 750 到 1000µm/s 的界面行為動態。

圖十二、拉速由 750 到 1000µm/s 的界面型態演化

在圖十二中，波長(~λc/4)在拉速上升之下幾乎不改變，且 Multiplets

175µm，甚至隨時間波長還有變大的趨勢，這代表深細胞區域已經接近 了，換句話說，Vλ2 ~常數的關係在此區域已經不適用。 圖十三、模擬和實驗所觀察到的非對秤的淺細胞界面生長 此外，我們在模擬和實驗也同時觀察到非對稱的淺細胞界面生長，此 現象首先由 Rappel [24]所報導，如圖十三所示，描述兩個細胞之間的時間 周期性競爭生長，它會出現在兩群 Multiplets 之間，且此搖擺頻率是拉速 的函數。

第四章 結論

在這份報告中我們描述了 SCN/Acetone 合金在臨界質附近的界面轉 換行為，包括定量的相場模擬和長時間尺度的實驗都算是一大突破。我們 的相場模擬在平整界面時與 MS 及 WL 理論相比，展現了良好的一致性， 而長時間且固定拉速的實驗也讓我們完成了從未被報導過的完整界面轉 換現象，我們把模擬所求得的波長和拉速與古典的 MS 理論比較，如圖十 四： 圖十四、模擬所求得的波長和拉速之關係與 MS 理論之比較

在上圖中，我們可以很清楚的看到四段明顯的型態變化，隨著拉速的 上升，界面會由平整界面到λc↔λc/2 細胞狀界面的轉換，再到 Vλ2~常數 的淺細胞界面，最後到 Multiplets 界面。除此之外，許多有趣的非線性動 態行為也清楚的呈現出來。 希望在不久的將來，透過主導方程式及模式的修正，我們可以達成使 用相場模式來模擬真實的問題，來與實驗作更精確的比較。

參考文獻

1. W.W. Mullins and R.F. Sekerka, J. Appl. Phys., 35, 444, 1964. 2. J.A. Warren and J.S. Langer, Phys. Rev. E, 47, 4072, 1993.

3. M.A. Eshelman, V. Seetharaman and R. Trivedi, Acta Metall., 36, 1165, 1988 4. M.A. Eshelman and R. Trivedi, Acta Metall., 35, 2443, 1987

5. P. Kurowski, S. de Cheveigne, G. Faivre, and C. Guthmann, J. Phys. (Paris), 50, 3007, 1989

6. L.H. Ungar and R.A. Brown, Phys. Rev. B, 29, 1367, 1984 7. K. Tsiveriotis and R.A. Brown, Phys. Rev. B, 48, 13495, 1993.

8. J.T.C. Lee, K. Tsiveriotis, and R.A. Brown, J. Cryst. Growth, 121, 536, 1992 9. C.W. Lan, C.J. Shih, and W.T. Hsu, J. Cryst. Growth, 264, 379, 2004

10. N.A. Ahmad, A.A. Wheeler, W.J. Boettinger, G.B. McFadden, Phys. Rev. E, 58, 3436 (1998).

11. W.J. Boettinger and J.A. Warren, J. Cryst. Growth, 200, 583, 1999 12. Z. Bi and R.F. Sekerka, J. Cryst. Growth, 237, 138, 2002

13. N. Provatas, N. Goldenfield, and J. Dantzig, J. Comput. Phys., 148, 265, 1999 14. C.W. Lan, C.C. Liu, and C.M. Hsu, J. Comp. Phys., 178, 464, 2002

15. M.J. Aziz and T. Kaplan, Acta. Metall., 36, 2335, 1988. 16. A. Karma, Phys. Rev. Letts, 87, 115701, 2001

17. C.W. Lan and C.J. Shih, Phys. Rev. E, 69, 031601, 2004

18. J.A. Warren and W.J. Boettinger, Acta Metall. Mater. 43, 689, 1995

19. W. Losert, B.Q. Shi, and H.Z. Cummins, Proc. Natl. Acad. Sci. USA, 95, 431, 1998 20. R. Trivedi and K. Somboonsuk, Acta Metall., 33, 1061, 1985

21. D.J. Wollkind and L.A. Segel, Philos. Trans. Roy. Soc. London 268, 351, 1970 22. H. Jamgotchian, R. Trivedi, and B. Billia, Phys. Rev. E., 47, 4313, 1993

23. P. Kopczynski, W.J. Rappel, and A. Karma, Phys. Rev. E, 55, R1282, 1997 24. W.J. Rappel, Phys. Rev. E, 48, 4118, 1993