数理解析研究所講究録
1129 巻 2000 年 29-37

29

Numerical Analysis of$\mathrm{E}\mathrm{g}\mathfrak{U}\mathrm{C}\mathrm{h}\mathrm{i}_{-}\mathrm{o}\mathrm{k}\mathrm{i}- \mathrm{M}\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{u}\mathrm{m}\mathrm{u}\mathrm{r}\mathrm{a}$Model(for Phase Separation)
花田孝郎 (HANADA Takao, Chiba Institute of Technology)

中村正彰 (NAKAMURA MasaAki, Nihon University)
島近義 (SHIMA Chikayoshi, Nihon University)
Abstract
In thermodynamics, phase$\mathrm{J}.\mathrm{W}$.Cahn and$\mathrm{J}.\mathrm{E}$ations in binary alloys are interesting phenomena.
.Hilliard introduced the free energy and derived the famous Cahn$\mathrm{s}\mathrm{e}\mathrm{p}_{C7}x$

Hilliard equation to analyze the spinodal decomposition. The relative minimizers forthe free energy are very interesting in mathematics, especially whose instability gives

difficult problems to makers of several products using alloys.
T. Eguchi, K. Oki, and S. Matsumura introduced the degree of order in binary
alloys adding to the concentration of components to investigate the kinetics of phaseseparations. Using this model, we shall show that the local concentration begins todiverge by small perturbations in the degree of order though the distribution at thebeginning is homogeneous.

1Introduction

We introduce the Eguchi-Oki-Matsumura equation describing a phase separation for a substitutional binary alloy $A_{1-m1m}B+$ consisting of A and atoms filled in a vessel. In thecontinuum theory, the local concentration is first introduced to be conserved as$\mathrm{B}$$u$$\frac{1}{|\Omega|}\int_{\Omega}u(t, X)dX=m$,(1)

where is a bounded domain in$n=1,2,3$ , with the smooth boundary. Next thelocal degree of order is introduced to describe the thermodynamic potential of this system$\Omega$$\mathbb{R}^{n},$$v$$\partial\Omega$

30
as$F(u, v)= \int_{\Omega}(f(u, v)+\frac{1}{2}H|\nabla u|^{2}+\frac{1}{2}K|\nabla v|^{2})dx$.

(2)Here, $f(u, v)$ is the density of the bulk free energy,$f(u, v)= \frac{1}{2}au^{2}-\frac{1}{2}bv^{2}+\frac{1}{4}b_{1}v^{4}+\frac{1}{2}gu^{2}v^{2}$(3),are positive constants, and is a physical parameter depending on the temperature such as to be positive only below the critical one. And $K,$ $H$ are the surface energy
per unit area considered to be positive constants. Then the equations of equilibriurn state

where$a,$$b_{1},$$b$$g$

are given as follows,$\frac{\partial f(u,v)}{\partial u}=\mu,$ $\frac{\partial f(u,v)}{\partial v}=0$

is a chemical potential.
Hence we obtain the kinetic equations forwhere(4),$\mu$

and ,

$u$

$v$

$L^{-1} \frac{\partial u}{\partial t}=-\nabla/^{2}(H\nabla^{2}u-\frac{\partial f(u,v)}{\partial u})$(5,in$\frac{\partial v}{\partial t}=K\nabla^{2}v-\frac{\partial f(u,v)}{\partial v}$$\Omega$(6,

is the coefficient of diffusion speed of the material to one of the degree of order.
These equations are analyzed with the boundary conditions

where$L$,(7)$=0$ ,(8)$0$$=$$\nu\cdot\nabla u(t)$$\nu\cdot\nabla^{3}u(t)$$\nu\cdot\nabla v(\theta)$

$=$

$0$

on

$\partial\Omega$(9),such that the equation (1) is satisfied, and the initial conditions2$u(0)$$=$$u_{o}$$v(0)$$=$$v_{o}$

(10)

,

in

$\Omega$(11).

Mathematical Results

We can show the well posedness of the problem based on the$\mathrm{E}\mathrm{g}\mathrm{u}\mathrm{c}\mathrm{h}\mathrm{i}-\mathrm{O}\mathrm{k}\mathrm{i}- \mathrm{M}\mathrm{a}\mathrm{t}_{\mathrm{S}}\mathrm{u}\mathrm{m}\mathrm{u}\mathrm{r}\mathrm{a}$model.

31

Theorem 1 For any $T>0$ and any$(u_{o}, v)\mathit{0}\in L_{2}\cross(L_{2}\cap L_{4})_{f}$there exist a solution$(u, v)\in C_{w}(0,T;L2)sati_{S}Mng$$u\in L^{\infty}(0,\tau;L_{2}1\cap L_{2}(0, T;H^{1})$,$v\in L^{\infty}(0,\tau;L_{2}\cap L_{4})\cap L_{2}(0,\tau;V)\cap L_{6}(0, T;L6)$(12)$\frac{1}{L}\frac{d}{dt}(u, \phi)=-H(\nabla^{2}u, \nabla^{2}\phi)+(\frac{\partial f(u,v)}{\partial u}, \nabla^{2}\phi)$$\forall\phi\in H^{2}(\Omega),$ $\nu\cdot\nabla u=0$

on$\partial\Omega$$\frac{d}{dt}(v, \psi)=-K(\nabla v, \nabla\psi)+(\frac{\partial f(u,v)}{\partial v}, \psi)$$(u(\mathrm{o}), v(\mathrm{O}))=(u_{o},v)\mathit{0}$$\forall\psi_{\in H^{1}}(\Omega)$

.

(14)

The solutions homogeneous in space

2.1

If the distributions$u,$ $v$

are constant in space, the equations (5-6) is reduced to$\dot{u}=0$,

(15)
(16)

$\dot{v}=(gu^{2}+v-2b)v$ ,where$\beta=0$(13)$\dot{u},\dot{v}$denote $du/dt.dv/dt$ . Since we have $u(t)=m$ , let $\beta=b-gm^{2}$ . Then we have, if,$v(t)= \frac{v_{o}}{(2tv_{O}^{2}+1)^{\frac{1}{2}}}$,(17)or otherwise$v(i)=v_{o}( \frac{\beta}{?J_{o^{2}}-(v^{2}-\circ\beta)\mathrm{e}_{\lrcorner}\mathrm{x}\mathrm{P}^{-2}\beta t})^{\frac{1}{2}}$.

(18)

In more realistic situations, the quantity increases as the alloy gradually cooled down.
So the growth, or changing of in time must be considered strictly. However, in this paper,
we separate the heat convection in the alloy from the model of phase separations.
First, if the temperature of the binary alloy is still high enough for, then we have
$b$

$b$

$\beta\leq 0$$\lim_{tarrow\infty}v(t)=0$.

32$\succ$Figure 1:$\mathrm{E}\mathrm{v}o$lution of homogeneous solutions, $\beta<0,$ $\beta=0,$ $\beta>0$ , respectively

Next, if the temperature is below the critical value such that $b>0$ , there exist othersituations such that$\lim_{tarrow\infty}|_{U(t)|=}\sqrt{\beta}$,

if $v(t),\neq 0$ at some . In this case, $v(t)=0$ still satisfy the equations (15-16), but it is$t$unstable$(\mathrm{F}\mathrm{i}\mathrm{g}.2^{\backslash })$.

Figure 2: Bifurcation of uniform solutions $v(t)$ inAs a stationary state of$u$

and$|v|$, let$(m,0\rangle (\beta\leq 0)$,$(\beta>0)$.

$(\overline{u}(\beta),\overline{v}(\beta))=\{$$(\gamma\gamma l, \sqrt{\beta})$Then, we have$tarrow\infty \mathrm{i}\mathrm{i}\mathrm{m}v(t)=\pm\overline{v}(\beta)$

in accordance with the sign of$v(\mathrm{O})$.

, or ,
$0$

$b$

33

3Problems of One-dimension in space
In this section, the problems of and depending onstudied. Then the equations (5-6) is reduced to$u$$v$$\mathrm{t}\mathrm{i}_{1}\mathrm{n}\mathrm{e}t$

and only one direction

$x$

are,

(19)

$\dot{v}=Kv’’-(gu^{2}+v^{2}-b)v$ ,(20)$L^{-1}\dot{u}=-(Hu^{\prime/}-(a+gv^{2})u)^{\prime/}$, etc. denote $du/dx,$, etc.,.
For the problem (19-20), there are not only homogeneous solutions described in the previous section, but also nonhomogeneous ones. In order to examine whether the homogeneous
solutions are stable or not, and look for other solutions, numerical simulations are used.
where3.1$d^{2}u/dx^{2}$$u_{\mathrm{z}}’\backslash u’f$$\mathrm{r}\mathrm{e}\mathrm{s}_{\mathrm{P}}\mathrm{e}\mathrm{C}\mathrm{t}\mathrm{i}_{\mathrm{V}\mathrm{e}\mathrm{I}}\mathrm{y}$Discretized Schemes

The equations (19-20) in the domain $(0, \infty)\cross(0,1)$ is discretized with forward differences$\triangle x=1/n$ , then the equations forin time and central differences in space. Let
approximations . of $u(t, x_{k})$ andof $v(t,X_{k})$ are as following:$x_{k}=k\triangle x,$$U_{k}$$V_{k}$$L^{-1_{\frac{\overline{U}_{k}-U_{k}}{\triangle l}=-}}H \frac{U_{k-2}-4c\Gamma-1+k6U_{k}-4U_{k+}1+Uk+2}{\triangle x^{4}}$$+ \frac{(a+gV_{k-}1^{2})Uk-1-2(a+gV_{k}^{2})U_{k}+(a+gVk+1)2U_{k}+1}{\triangle x^{2}}..$

,

where
and
are the approximations at the next time step atconditions $u’=0,$
are discretized as$\overline{U_{k}}$

(21)
(22)

$\frac{\overline{V}_{k}-V_{k}}{\triangle t}=I\sigma\frac{V_{k-1}-2Vk+Vk+1}{\triangle x^{2}}-(gU^{2}k.+V_{k}^{2}-b)V_{k}$$(0\leq k\leq n)$,$\overline{V_{k}}$$t+\triangle t$. The boundary$u^{f}\prime\prime=0$$U_{-2}=U_{2},$$U_{-1}=U_{1}$,$U_{n-1}=U_{n+1},$$U_{n-2}=U_{n+2}$

,

(23)

and $v’=0$ is as$V_{-1}=V_{1},$$V_{n-1}=V_{n+1}$.

(24)

In numerical simulations taking constants as$L=1/1024,$ $H=K=1/1000,$ $m=0.25,$ $a=0.25,g=8,$ $b_{1}=1$ ,the dependence of solutions on , and the stability of the homogeneous states are studied.
$b$

34

3.2

Instability of the homogeneous solutions

In this’subsection, we show that the homogeneous solutions are not stable if $\beta>0$ . For that
purpose, initial conditions are perturbed slightly from the homogeneous solutions.

3.2.1

Zero solution $v=0$

Taking initial conditions as$u(0, x)=$$m$,$v(0,x)=\epsilon\cos(\pi x)$

for small$|\epsilon|$.

The solutions, especially $u(t)$ , are firstly going away from the initial state for all values of, then we may conclude
if. Then, $u(t)$ is converging to the constant functionthat the homogeneous zero solution is stable. However, if $\beta>0$ , then the perturbation for
$m$

$b$

$\beta\leq 0$$\mathrm{v}$$\mathrm{u}$Figure 3: Behavior of $u(t, x),$ $v(t, X)$ (evolving from behind)$v$

is expanding to the values$\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\pm\sqrt{b}$, and it is followed by the separation of the valuesof (which is observed as the phase separation phenomena). It is shown by mathematicalanalysis in the previous subsection 2.1. In Figure 3, the simulated results are shown in
$u$

35

the case of $b=0.99$ (, therefore $\beta=0.49$ ), $\epsilon=0.001,$$0\leq t\leq 128$ .

3.2.2

Instability of$\triangle t=$1/256 and $\triangle x=1/64$ , while$v=\sqrt{\beta}$

Taking a value of such that $\beta>0$ , there exists the solution $(u, v)=(m.\sqrt{\beta})$ . Therefore,
it is important to simulate solutions starting from the initial states near that homogeneoussolution as$b$,$u(0,x)=$$m$$v(0, x)=$$\sqrt{\beta}+\epsilon\cos(\pi x)$.

For small such thatis positive, and $v(t, x)$ is kept to be positive also. In
this case, $v(t)$ approaches to a function oscillating between a number nearand another$\epsilon$$\epsilon^{2}<\beta,$$v(\mathrm{O}, x)$$\sqrt{b}$positive near zero. Then, oscillations of $u(t)$ appears according to ones of $v(t)$ . In Figure 4,
the simulated results are shown in the same parameters, while $0\leq t\leq 4096$ . However, thenumber of oscillations of $v(t)$ and $u(t)$ varies from one at first, to three secondly, and to one
after a long time. Also other simulations result in the same solution having one oscillation,other than having many oscillations by syrrumetries.

4Conclusion

In numerical simulations, we have observed that some homogeneous solutions are not stablein one-dimensional problems. Then, it is shown that the phenomena of phase separation inbinary alloy are caused by perturbation only in the order parameter.
Since it is thought that the stationary problems in the.M. model have many solutions
stable or not, the structure of all solutions is left to be made clear.
$\mathrm{E}.\mathrm{O}$

References[1] Alikakos, A., Bates, P.W., and Fbsco, G.: Slow motion for the Cahn Hilliard equationin one space dimension, Journal of Differential Equations 90, 1991, pp.81-135.
[2] Alt, H.W. and Pawlow, I.: Dynamics of non-isothermal phase transition, Internationalser. Numer. Math. 95(1990), pp.1-26.

36$\mathrm{u}$$\mathrm{v}$Figure 4: Behavior of $u(t.$ x), $v(t,$X)

37[3] Carr. J. and Pego, L.: Metastable$pab\partial ern\mathit{8}$in solutions of $u_{t}=\epsilon^{2}u_{xx}-f(u)$ , Communi-cations Pure and Applied Mathematics, 42(1989), pp.523-576.
[4] Cahn,$\mathrm{J}.\mathrm{W}.$: On$\mathit{8}pinodal$

decomposition, Acta Metal, 9(1961) pp.795-801.

. and Hilliard,energy of nonuniform
energy, J. Chem. Phys. 28(1961) pp.259-267.

[5] Cahn,$\mathrm{J}.\mathrm{E}.:Free$$\mathrm{J}.\mathrm{W}$[6] Eguchi T., Oki K. and Matsumura$\mathrm{S}.:KineticS$$\mathit{8}y_{\mathit{8}te}mS$

I. Interfacial freeof ordering with phase separation, Math.

Res. Soc. Symp. Proc. 21(1984) pp.589-594.

. and French,of the Cahn-Hilliard equation forseparation, IMA J. Appl. Math. 38(1987) pp.97-128.

[7] Elliott,$pha\mathit{8}e$$\mathrm{D}.\mathrm{A}.:Numerical\mathit{8}tudies$$\mathrm{C}.\mathrm{M}$[8] Furihata Daisuke and Mori Masatake: General demeans of a$di_{\mathit{8}Cr}ete$[9] Fusco, G. and Hale,$7\dot{T}vation$of finite difference schemes byvariation, Trans. of the Japan SIAM 8(1998) pp.317-340.

motion manifold, dormant instability and singular perturbation, Journal of Dynamical and Differential Equations, 1 (1989) pp.57-94.
$\mathrm{J}.\mathrm{K}.:Sl_{\mathit{0}}w$. and Matano, H.: On the structure of equilibrium phase transitions with thegrad\’ient theory of fiuids, Quarterly Appl. Math. 46(1988) pp.301-317.

[10] Gurtin,$\mathrm{M}.\mathrm{E}$[11] Hanada Takao, Nakamura MasaAki and Shima Chikayoshi: On $Eguchi- Oki-Ma\theta Sumura$equations, The fourth Japan-Chinajoint seminar on numerical mathematics (to appear)[12] Qiang Do and Nicolaides,:Numerical analysis of a continuum model of phase transition, SIAM J. Numer. Anal. 28(1991) pp.1310-1322.
$\mathrm{R}.\mathrm{A}.$