This post builds on the formulations I showed in my previous posts by
introducing their nonlinear versions .
In a typical nonlinear problem, the variational setting leads to the weak formulation
Find u ∈ V u\in V u ∈ V
F ( u , v ) = 0 \begin{equation}
F(u,v) = 0
\htmlId{eq:femnonlinear2}{}
\tag{1}\end{equation} F ( u , v ) = 0 ( 1 ) for all v ∈ V v\in V v ∈ V F F F u u u v v v
We linearize F F F
L i n [ F ( u , v ) ] u = u ˉ = F ( u ˉ , v ) + D u F ( u , v ) ⋅ Δ u ∣ u = u ˉ \begin{equation}
\Lin [F(u,v)]_{u=\bar{u}} = F(\bar{u}, v)
+ \varn{F(u,v)}{u}{\Var u}\evat_{u=\bar{u}}
\htmlId{eq:femnonlinear4}{}
\tag{2}\end{equation} Lin [ F ( u , v ) ] u = u ˉ = F ( u ˉ , v ) + D u F ( u , v ) ⋅ Δ u u = u ˉ ( 2 ) Equating (2) to zero yields a linear system in terms of
Δ u \Var u Δ u
a ( Δ u , v ) = b ( v ) \begin{equation}
\boxed{
a(\Var u, v) = b(v)
}
\htmlId{eq:femnonlinear6}{}
\tag{3}\end{equation} a ( Δ u , v ) = b ( v ) ( 3 ) where
a ( Δ u , v ) = D u F ( u , v ) ⋅ Δ u ∣ u = u ˉ b ( v ) = − F ( u ˉ , v ) . \begin{equation}
\boxed{
\begin{aligned}
a(\Var u, v) &= \varn{F(u,v)}{u}{\Var u}\evat_{u=\bar{u}} \\
b(v) &= -F(\bar{u}, v)
.
\end{aligned}
}
\htmlId{eq:femnonlinear5}{}
\tag{4}\end{equation} a ( Δ u , v ) b ( v ) = D u F ( u , v ) ⋅ Δ u u = u ˉ = − F ( u ˉ , v ) . ( 4 ) We can compute the components of the matrices and vectors according to (3)
A I J = a ( N J , N I ) = D u F ( u , N I ) ⋅ N J ∣ u = u ˉ b I = b ( N I ) = − F ( u ˉ , N I ) . \begin{equation}
\boxed{
\begin{alignedat}{3}
\Aelid{I\!J}{} &= a(N^J,N^I)
&&= \varn{F(u,N^I)}{u}{N^J}\evat_{u=\bar{u}} \\
b^{I} &= b(N^I)
&&= -F(\bar{u}, N^I).
\end{alignedat}
}
\htmlId{eq:femnonlinear9}{}
\tag{5}\end{equation} A I J b I = a ( N J , N I ) = b ( N I ) = D u F ( u , N I ) ⋅ N J u = u ˉ = − F ( u ˉ , N I ) . ( 5 ) Then the update vector Δ u = [ Δ u 1 , Δ u 2 , … , Δ u n n ] T \Var \Bu = [\Var u^1, \Var u^2, \dots, \Var
u^{\nnode}]\tra Δ u = [ Δ u 1 , Δ u 2 , … , Δ u n n ] T
A Δ u = b \begin{equation}
\BA \Var \Bu = \Bb
\tag{6}\end{equation} A Δ u = b ( 6 ) Letting Δ u \Var u Δ u
u ← u ˉ + Δ u \begin{equation}
\boxed{
\Bu \leftarrow \bar{\Bu} + \Var\Bu
}
\tag{7}\end{equation} u ← u ˉ + Δ u ( 7 ) Example: Nonlinear Poisson’s Equation
Consider the following nonlinear Poisson’s equation
− ∇ ⋅ ( g ( u ) ∇ u ) = f in Ω u = 0 on ∂ Ω \begin{equation}
\begin{alignedat}{4}
- \nabla \dtp (g(u)\nabla u) &= f \quad && \text{in} \quad && \Omega \\
u &= 0 \quad && \text{on} \quad && \del\Omega
\end{alignedat}
\htmlId{eq:femnonlinear8}{}
\tag{8}\end{equation} − ∇ ⋅ ( g ( u ) ∇ u ) u = f = 0 in on Ω ∂ Ω ( 8 ) The weak formulation reads
Find u ∈ V u\in V u ∈ V
− ∫ Ω ∇ ⋅ ( g ( u ) ∇ u ) v d v = ∫ Ω f v d v \begin{equation}
- \int_\Omega \nabla \dtp (g(u)\nabla u) v \dv= \int_\Omega f v \dv
\tag{9}\end{equation} − ∫ Ω ∇ ⋅ ( g ( u ) ∇ u ) v d v = ∫ Ω f v d v ( 9 ) for all v ∈ V v\in V v ∈ V V = H 0 1 ( Ω ) V=H^1_0(\Omega) V = H 0 1 ( Ω )
Applying integration by parts and divergence theorem on the left-hand side
∫ Ω ∇ ⋅ ( g ( u ) ∇ u ) v d v = ∫ Ω ∇ ⋅ ( g ( u ) ∇ ( u ) v ) d v − ∫ Ω g ( u ) ∇ u ⋅ ∇ v d v = ∫ ∂ Ω g ( u ) v ( ∇ u ⋅ n ) d a ⏟ v = 0 on ∂ Ω − ∫ Ω g ( u ) ∇ u ⋅ ∇ v d v \begin{equation}
\begin{aligned}
\int_\Omega \nabla \dtp (g(u)\nabla u) v \dv
&= \int_\Omega \nabla \dtp (g(u)\nabla (u) v) \dv
- \int_\Omega g(u)\nabla u\dtp\nabla v \dv \\
&= \underbrace{\int_{\del\Omega} g(u) v (\nabla u\dtp\Bn) \da}_{v = 0
\text{ on } \del\Omega}
- \int_\Omega g(u)\nabla u\dtp\nabla v \dv \\
\end{aligned}
\tag{10}\end{equation} ∫ Ω ∇ ⋅ ( g ( u ) ∇ u ) v d v = ∫ Ω ∇ ⋅ ( g ( u ) ∇ ( u ) v ) d v − ∫ Ω g ( u ) ∇ u ⋅ ∇ v d v = v = 0 on ∂ Ω ∫ ∂ Ω g ( u ) v ( ∇ u ⋅ n ) d a − ∫ Ω g ( u ) ∇ u ⋅ ∇ v d v ( 10 ) Thus we have the semilinear form
F ( u , v ) = ∫ Ω g ( u ) ∇ u ⋅ ∇ v d v − ∫ Ω f v d v = 0 \begin{equation}
F(u,v) = \int_{\Omega} g(u) \nabla u \dtp \nabla v \dv - \int_{\Omega} f \,
v \dv = 0
\tag{11}\end{equation} F ( u , v ) = ∫ Ω g ( u ) ∇ u ⋅ ∇ v d v − ∫ Ω f v d v = 0 ( 11 ) The linearized version of this problem is then with (4)
a ( Δ u , v ) = ∫ Ω ( d g d u ∣ u ˉ Δ u ∇ u ˉ + g ( u ˉ ) ∇ ( Δ u ) ) ⋅ ∇ v d v b ( v ) = ∫ Ω [ f v − g ( u ˉ ) ∇ u ˉ ⋅ ∇ v ] d v \begin{equation}
\begin{aligned}
a(\Var u,v) &= \int_{\Omega}
\rbr{\deriv{g}{u}\evat_{\bar{u}} \Var u\, \nabla \bar{u}
+ g(\bar{u})\nabla(\Var u)} \dtp \nabla v \dv \\
b(v) &= \int_{\Omega} [f \, v - g(\bar{u}) \nabla \bar{u} \dtp \nabla v] \dv\\
\end{aligned}
\tag{12}\end{equation} a ( Δ u , v ) b ( v ) = ∫ Ω ( d u d g u ˉ Δ u ∇ u ˉ + g ( u ˉ ) ∇ ( Δ u ) ) ⋅ ∇ v d v = ∫ Ω [ f v − g ( u ˉ ) ∇ u ˉ ⋅ ∇ v ] d v ( 12 ) and the matrix and vector components are with (5)
A I J = ∫ Ω ( d g d u ∣ u ˉ N J ∇ u ˉ + g ( u ˉ ) B J ) ⋅ B I d v b I = ∫ Ω [ f N I − g ( u ˉ ) ∇ u ˉ ⋅ B I ] d v \begin{equation}
\begin{aligned}
\Aelid{I\!J}{} &= \int_{\Omega}
\rbr{\deriv{g}{u}\evat_{\bar{u}} N^J \, \nabla \bar{u}
+ g(\bar{u})\BB^J} \dtp \BB^I \dv \\
b^{I} &= \int_{\Omega} [f \, N^I - g(\bar{u}) \nabla \bar{u} \dtp \BB^I ] \dv\\
\end{aligned}
\tag{13}\end{equation} A I J b I = ∫ Ω ( d u d g u ˉ N J ∇ u ˉ + g ( u ˉ ) B J ) ⋅ B I d v = ∫ Ω [ f N I − g ( u ˉ ) ∇ u ˉ ⋅ B I ] d v ( 13 ) where the previous solution and its gradient are computed as
u ˉ = ∑ I = 1 n n u ˉ I N I and ∇ u ˉ = ∑ I = 1 n n u ˉ I B I . \begin{equation}
\bar{u} = \suml{I=1}{\nnode} \bar{u}^I N^I
\eqand
\nabla \bar{u} = \suml{I=1}{\nnode} \bar{u}^I \BB^I
.
\tag{14}\end{equation} u ˉ = I = 1 ∑ n n u ˉ I N I and ∇ u ˉ = I = 1 ∑ n n u ˉ I B I . ( 14 ) Nonlinear Time-Dependent Problems
In the case of a nonlinear time-dependent problem, we have the following weak
form:
Find u ∈ V u \in V u ∈ V
m ( u ˙ , v ; t ) + F ( u , v ; t ) = 0 \begin{equation}
m(\dot{u}, v; t) + F(u,v; t) = 0
\htmlId{eq:nonlineartimedependentweak1}{}
\tag{15}\end{equation} m ( u ˙ , v ; t ) + F ( u , v ; t ) = 0 ( 15 ) for all v ∈ V v \in V v ∈ V t ∈ [ 0 , ∞ ) t \in [0,\infty) t ∈ [ 0 , ∞ ) F F F
Discretization yields the following nonlinear system of equations
M ( t ) u + f ( u ; t ) = 0 \begin{equation}
\BM(t)\Bu + \Bf(u; t) = \Bzero
\tag{16}\end{equation} M ( t ) u + f ( u ; t ) = 0 ( 16 ) where
M I J ( t ) = m ( N J , N I ; t ) f I ( u ; t ) = F ( u , N I ; t ) . \begin{equation}
\begin{aligned}
M^{I\!J}(t) &= m(N^J, N^I; t) \\
f^{I}(u;t) &= F(u, N^I; t).
\end{aligned}
\tag{17}\end{equation} M I J ( t ) f I ( u ; t ) = m ( N J , N I ; t ) = F ( u , N I ; t ) . ( 17 ) Explicit Euler Scheme
We discretize in time with the finite difference
u ˙ ≈ [ u n + 1 − u n ] / Δ t \dot{u} \approx [u_{n+1}-u_n]/{\Delta t} u ˙ ≈ [ u n + 1 − u n ] / Δ t
\begin{equation}
\boxed{
m(\dot{u}, v; t)
\approx \frac{1}{\Delta t} [m(u_{n+1}, v; t_{n+1}) - m(u_n, v; t_n)]
}
\label{eq:discretetimedependent1}
\end{equation}
We discretize the variational forms in time according to
\eqref{eq:discretetimedependent1}, and evaluate the remaining terms at t n t_n t n
1 Δ t [ m ( u n + 1 , v ; t n + 1 ) − m ( u n , v ; t n ) ] + F ( u n , v ; t n ) = 0 \begin{equation}
\frac{1}{\Delta t} [m(u_{n+1},v;t_{n+1}) - m(u_{n},v;t_{n})] + F(u_n, v; t_n) = 0
\tag{18}\end{equation} Δ t 1 [ m ( u n + 1 , v ; t n + 1 ) − m ( u n , v ; t n )] + F ( u n , v ; t n ) = 0 ( 18 ) The corresponding system of equations is
1 Δ t [ M n + 1 u n + 1 − M n u n ] + f n = 0 \begin{equation}
\frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n]
+ \Bf_n = \Bzero
\tag{19}\end{equation} Δ t 1 [ M n + 1 u n + 1 − M n u n ] + f n = 0 ( 19 ) where f n = f ( u n , t n ) \Bf_n = \Bf(u_n, t_n) f n = f ( u n , t n )
u n + 1 = M n + 1 − 1 [ M n u n − Δ t f n ] \begin{equation}
\boxed{
\Bu_{n+1} = \BM_{n+1}\inv [\BM_n\Bu_n - \Delta t \Bf_n]
}
\tag{20}\end{equation} u n + 1 = M n + 1 − 1 [ M n u n − Δ t f n ] ( 20 ) For a time-independent m m m
u n + 1 = u n − Δ t M − 1 f n \begin{equation}
\Bu_{n+1} = \Bu_n - \Delta t \BM\inv\Bf_n
\tag{21}\end{equation} u n + 1 = u n − Δ t M − 1 f n ( 21 ) Implicit Euler Scheme
For the implicit scheme, we evaluate the remaining terms at t n + 1 t_{n+1} t n + 1
G ( u n + 1 , v ) : = 1 Δ t [ m ( u n + 1 , v ; t n + 1 ) − m ( u n , v ; t n ) ] + F ( u n + 1 , v ; t n + 1 ) = 0 \begin{equation}
G(u_{n+1}, v)
:= \frac{1}{\Delta t} [m(u_{n+1},v;t_{n+1}) - m(u_{n},v;t_{n})] + F(u_{n+1}, v; t_{n+1}) = 0
\tag{22}\end{equation} G ( u n + 1 , v ) := Δ t 1 [ m ( u n + 1 , v ; t n + 1 ) − m ( u n , v ; t n )] + F ( u n + 1 , v ; t n + 1 ) = 0 ( 22 ) We will hereon replace u n + 1 u_{n+1} u n + 1 u u u G ( u , v ) G(u, v) G ( u , v )
L i n [ G ( u , v ) ] u = u ˉ = G ( u ˉ , v ) + D u G ⋅ Δ u ∣ u = u ˉ = 0 \begin{equation}
\Lin[G(u,v)]_{u=\bar{u}}
= G(\bar{u}, v) + \varn{G}{u}{\Var u}\evat_{u=\bar{u}} = 0
\tag{23}\end{equation} Lin [ G ( u , v ) ] u = u ˉ = G ( u ˉ , v ) + D u G ⋅ Δ u u = u ˉ = 0 ( 23 ) We thus have the following linear setting for the Newton update Δ u \Var u Δ u
a ( Δ u , v ) = b ( v ) \begin{equation}
a(\Var u, v) = b(v)
\tag{24}\end{equation} a ( Δ u , v ) = b ( v ) ( 24 ) where
a ( Δ u , v ) : = D u G ⋅ Δ u ∣ u = u ˉ = 1 Δ t m ( Δ u , v ; t n + 1 ) + D u F ( u , v ; t n + 1 ) ⋅ Δ u ∣ u = u ˉ b ( v ) : = − G ( u ˉ , v ) = − F ( u ˉ , v ; t n + 1 ) − 1 Δ t [ m ( u ˉ , v ; t n + 1 ) − m ( u n , v ; t n ) ] \begin{equation}
\begin{aligned}
a(\Var u, v)
&:= \varn{G}{u}{\Var u} \evat_{u=\bar{u}}
= \frac{1}{\Delta t} m(\Var u, v; t_{n+1})
+ \varn{F(u, v; t_{n+1})}{u}{\Var u} \evat_{u=\bar{u}} \\
b(v) &:= -G(\bar{u}, v)
= - F(\bar{u}, v; t_{n+1})
-\frac{1}{\Delta t} [m(\bar{u},v;t_{n+1}) - m(u_{n},v;t_{n})]
\end{aligned}
\tag{25}\end{equation} a ( Δ u , v ) b ( v ) := D u G ⋅ Δ u u = u ˉ = Δ t 1 m ( Δ u , v ; t n + 1 ) + D u F ( u , v ; t n + 1 ) ⋅ Δ u u = u ˉ := − G ( u ˉ , v ) = − F ( u ˉ , v ; t n + 1 ) − Δ t 1 [ m ( u ˉ , v ; t n + 1 ) − m ( u n , v ; t n )] ( 25 ) Discretization yields
( 1 Δ t M n + 1 + A ~ ) Δ u = b \begin{equation}
\rbr{\frac{1}{\Delta t} \BM_{n+1} + \tilde{\BA}}\Var \Bu
= \Bb
\tag{26}\end{equation} ( Δ t 1 M n + 1 + A ~ ) Δ u = b ( 26 ) where
A ~ I J : = D u F ( u , N I ; t n + 1 ) ⋅ N J ∣ u = u ˉ and b I : = b ( N I ) \begin{equation}
\tilde{A}^{I\!J} := \varn{F(u, N^I;t_{n+1})}{u}{N^J} \evat_{u=\bar{u}}
\eqand
b^I := b(N^I)
\tag{27}\end{equation} A ~ I J := D u F ( u , N I ; t n + 1 ) ⋅ N J u = u ˉ and b I := b ( N I ) ( 27 ) The Newton update is rendered
u ← u ˉ + Δ u and Δ u = [ 1 Δ t M n + 1 + A ~ ] − 1 b \begin{equation}
\boxed{
\Bu \leftarrow \bar{\Bu} + \Var\Bu
\eqwith
\Var \Bu
= [\frac{1}{\Delta t} \BM_{n+1} + \tilde{\BA}]\inv\Bb
}
\tag{28}\end{equation} u ← u ˉ + Δ u and Δ u = [ Δ t 1 M n + 1 + A ~ ] − 1 b ( 28 ) which is repeated until the solution for the next timestep u \Bu u
Nonlinear Coupled Problems
For a nonlinear coupled problem, the weak formulation is as follows
Find u ∈ V 1 u\in V_1 u ∈ V 1 y ∈ V 2 y\in V_2 y ∈ V 2
F ( u , y , v ) = 0 G ( u , y , w ) = 0 \begin{equation}
\begin{aligned}
F(u, y, v) &= 0 \\
G(u, y, w) &= 0 \\
\end{aligned}
\htmlId{eq:nonlinearcoupled1}{}
\tag{29}\end{equation} F ( u , y , v ) G ( u , y , w ) = 0 = 0 ( 29 ) for all v ∈ V 1 v\in V_1 v ∈ V 1 w ∈ V 2 w \in V_2 w ∈ V 2 F ( ⋅ , ⋅ , ⋅ ) F(\cdot,\cdot, \cdot) F ( ⋅ , ⋅ , ⋅ ) G ( ⋅ , ⋅ , ⋅ ) G(\cdot, \cdot, \cdot) G ( ⋅ , ⋅ , ⋅ ) u u u y y y v v v w w w
We linearize the semilinear forms about the nonlinear terms:
L i n [ F ( u , y , v ) ] u ˉ , y ˉ = F ( u ˉ , y ˉ , v ) + D u F ( u , y , v ) ⋅ Δ u ∣ u ˉ , y ˉ + D y F ( u , y , v ) ⋅ Δ y ∣ u ˉ , y ˉ L i n [ G ( u , y , w ) ] u ˉ , y ˉ = G ( u ˉ , y ˉ , w ) + D u G ( u , y , w ) ⋅ Δ u ∣ u ˉ , y ˉ + D y G ( u , y , w ) ⋅ Δ y ∣ u ˉ , y ˉ \begin{equation}
\begin{alignedat}{4}
\Lin[F(u, y, v)]_{\bar{u},\bar{y}}
&= F(\bar{u},\bar{y},v)
&&+ \varn{F(u, y, v)}{u}{\Var u} \evat_{\bar{u},\bar{y}}
&&+ \varn{F(u, y, v)}{y}{\Var y} \evat_{\bar{u},\bar{y}} \\
\Lin[G(u, y, w)]_{\bar{u},\bar{y}}
&= G(\bar{u},\bar{y},w)
&&+ \varn{G(u, y, w)}{u}{\Var u} \evat_{\bar{u},\bar{y}}
&&+ \varn{G(u, y, w)}{y}{\Var y} \evat_{\bar{u},\bar{y}}
\end{alignedat}
\htmlId{eq:nonlinearcoupled2}{}
\tag{30}\end{equation} Lin [ F ( u , y , v ) ] u ˉ , y ˉ Lin [ G ( u , y , w ) ] u ˉ , y ˉ = F ( u ˉ , y ˉ , v ) = G ( u ˉ , y ˉ , w ) + D u F ( u , y , v ) ⋅ Δ u u ˉ , y ˉ + D u G ( u , y , w ) ⋅ Δ u u ˉ , y ˉ + D y F ( u , y , v ) ⋅ Δ y u ˉ , y ˉ + D y G ( u , y , w ) ⋅ Δ y u ˉ , y ˉ ( 30 ) where the evaluations take place at u = u ˉ u=\bar{u} u = u ˉ y = y ˉ y=\bar{y} y = y ˉ
Equating the linearized residuals to zero, we obtain a linear system of the form
a ( u , v ) + b ( y , v ) = c ( v ) d ( u , w ) + e ( y , w ) = f ( w ) \begin{equation}
\begin{alignedat}{3}
a(u, v) &+ b(y, v) &&= c(v) \\
d(u, w) &+ e(y, w) &&= f(w) \\
\end{alignedat}
\htmlId{eq:coupledweakform1}{}
\tag{31}\end{equation} a ( u , v ) d ( u , w ) + b ( y , v ) + e ( y , w ) = c ( v ) = f ( w ) ( 31 ) with the bilinear forms a a a b b b d d d e e e c c c f f f
a ( Δ u , v ) : = D u F ( u , y , v ) ⋅ Δ u ∣ u ˉ , y ˉ b ( Δ y , v ) : = D y F ( u , y , v ) ⋅ Δ y ∣ u ˉ , y ˉ d ( Δ u , w ) : = D u G ( u , y , w ) ⋅ Δ u ∣ u ˉ , y ˉ e ( Δ y , w ) : = D y G ( u , y , w ) ⋅ Δ y ∣ u ˉ , y ˉ and c ( v ) : = − F ( u ˉ , y ˉ , v ) f ( w ) : = − G ( u ˉ , y ˉ , w ) \begin{equation}
\begin{gathered}
\begin{alignedat}{4}
a(\Var u, v) &:= \varn{F(u, y, v)}{u}{\Var u} \evat_{\bar{u},\bar{y}}
\quad &
b(\Var y, v) &:= \varn{F(u, y, v)}{y}{\Var y} \evat_{\bar{u},\bar{y}}
\\
d(\Var u, w) &:= \varn{G(u, y, w)}{u}{\Var u} \evat_{\bar{u},\bar{y}}
\quad &
e(\Var y, w) &:= \varn{G(u, y, w)}{y}{\Var y} \evat_{\bar{u},\bar{y}}
\end{alignedat}
\\
\text{and}
\\
\begin{aligned}
c(v) &:= -F(\bar{u},\bar{y},v) \\
f(w) &:= -G(\bar{u}, \bar{y}, w)
\end{aligned}
\end{gathered}
\tag{32}\end{equation} a ( Δ u , v ) d ( Δ u , w ) := D u F ( u , y , v ) ⋅ Δ u u ˉ , y ˉ := D u G ( u , y , w ) ⋅ Δ u u ˉ , y ˉ b ( Δ y , v ) e ( Δ y , w ) := D y F ( u , y , v ) ⋅ Δ y u ˉ , y ˉ := D y G ( u , y , w ) ⋅ Δ y u ˉ , y ˉ and c ( v ) f ( w ) := − F ( u ˉ , y ˉ , v ) := − G ( u ˉ , y ˉ , w ) ( 32 ) Discretizing as done in the previous section, we obtain the following linear system of
equations
[ A B D E ] [ Δ u Δ y ] = [ c f ] \begin{equation}
\begin{bmatrix}
\BA & \BB \\
\BD & \BE
\end{bmatrix}
\begin{bmatrix}
\Var \Bu \\ \Var \By
\end{bmatrix}
=
\begin{bmatrix}
\Bc \\ \Bf
\end{bmatrix}
\tag{33}\end{equation} [ A D B E ] [ Δ u Δ y ] = [ c f ] ( 33 ) whose solution yields the update values Δ u \Var \Bu Δ u Δ y \Var \By Δ y
u ← u ˉ + Δ u y ← y ˉ + Δ y . \begin{equation}
\begin{alignedat}{3}
\Bu &\leftarrow \bar{\Bu} &&+ \Var\Bu \\
\By &\leftarrow \bar{\By} &&+ \Var\By
.
\end{alignedat}
\tag{34}\end{equation} u y ← u ˉ ← y ˉ + Δ u + Δ y . ( 34 ) Example: Cahn-Hilliard Equation
The Cahn-Hilliard equation describes the process of phase separation, by which
the two components of a binary fluid spontaneously separate and form domains
pure in each component. The problem is nonlinear, coupled and time-dependent.
The IBVP reads
∂ c ∂ t = ∇ ⋅ ( M ∇ μ ) in Ω × I ∇ c ⋅ n = 0 on ∂ Ω × I ∇ μ ⋅ n = 0 on ∂ Ω × I c = c 0 in Ω , t = 0 μ = 0 in Ω , t = 0 \begin{equation}
\begin{alignedat}{4}
\partd{c}{t} &= \nabla\dtp(\BM\nabla \mu)
\qquad&& \text{in} \qquad&& \Omega\times I \\
\nabla c\dtp\Bn &= 0 && \text{on} && \del\Omega\times I\\
\nabla \mu\dtp\Bn &= 0 && \text{on} && \del\Omega\times I\\
c &= c_0 && \text{in} && \Omega, t = 0 \\
\mu &= 0 && \text{in} && \Omega, t = 0 \\
\end{alignedat}
\htmlId{eq:cahnhilliard1}{}
\tag{35}\end{equation} ∂ t ∂ c ∇ c ⋅ n ∇ μ ⋅ n c μ = ∇ ⋅ ( M ∇ μ ) = 0 = 0 = c 0 = 0 in on on in in Ω × I ∂ Ω × I ∂ Ω × I Ω , t = 0 Ω , t = 0 ( 35 ) where
μ = d f d c − ∇ ⋅ ( Λ ∇ c ) \begin{equation}
\mu = \deriv{f}{c} - \nabla\dtp(\BLambda\nabla c)
\htmlId{eq:cahnhilliard2}{}
\tag{36}\end{equation} μ = d c df − ∇ ⋅ ( Λ ∇ c ) ( 36 ) and t ∈ I = [ 0 , ∞ ) t\in I = [0,\infty) t ∈ I = [ 0 , ∞ )
c c c μ \mu μ f : c ↦ f ( c ) f: c \mapsto f(c) f : c ↦ f ( c ) M \BM M Λ \BLambda Λ
The fourth-order PDE governing the problem can be formulated as a coupled
system of two second-order PDEs with the variables c c c μ \mu μ (35)
and (36) .
The weak formulation then reads
Find c ∈ V 1 c \in V_1 c ∈ V 1 μ ∈ V 2 \mu\in V_2 μ ∈ V 2
∫ Ω ∂ c ∂ t v d x − ∫ Ω ∇ ⋅ ( M ∇ μ ) v d x = 0 ∫ Ω [ μ − d f d c ] w d x + ∫ Ω ∇ ⋅ ( Λ ∇ c ) w d x = 0 \begin{equation}
\begin{aligned}
\int_\Omega \partd{c}{t} v \dx
- \int_\Omega \nabla\dtp(\BM\nabla \mu) v \dx &=0 \\
\int_\Omega \sbr{\mu - \deriv{f}{c}} w \dx
+ \int_\Omega \nabla\dtp(\BLambda\nabla c) w\dx &= 0
\end{aligned}
\tag{37}\end{equation} ∫ Ω ∂ t ∂ c v d x − ∫ Ω ∇ ⋅ ( M ∇ μ ) v d x ∫ Ω [ μ − d c df ] w d x + ∫ Ω ∇ ⋅ ( Λ ∇ c ) w d x = 0 = 0 ( 37 ) for all v ∈ V 1 v \in V_1 v ∈ V 1 w ∈ V 2 w \in V_2 w ∈ V 2 t ∈ I t \in I t ∈ I
We discretize in time implicitly with ∂ c / ∂ t ≈ ( c n + 1 − c n ) / Δ t \del c/\del t \approx
(c_{n+1}-c_n)/\Var t ∂ c / ∂ t ≈ ( c n + 1 − c n ) /Δ t We also denote the values for the next timestep
c n + 1 c_{n+1} c n + 1 μ n + 1 \mu_{n+1} μ n + 1 c c c μ \mu μ
Using integration-by-parts, the divergence theorem, and the given boundary
conditions, we arrive at the following nonlinear forms
F ( c , μ , v ) = ∫ Ω 1 Δ t ( c − c n ) v d x + ∫ Ω ( M ∇ μ ) ⋅ ∇ v d x = 0 G ( c , μ , w ) = ∫ Ω [ μ − d f d c ] w d x − ∫ Ω ( Λ ∇ c ) ⋅ ∇ w d x = 0 \begin{equation}
\begin{alignedat}{3}
F(c,\mu,v) &= \int_\Omega \frac{1}{\Var t} (c-c_n) v \dx
+ \int_\Omega (\BM\nabla \mu)\dtp \nabla v \dx &&= 0 \\
G(c,\mu,w) &= \int_\Omega \sbr{\mu - \deriv{f}{c}} w \dx
- \int_\Omega (\BLambda\nabla c)\dtp \nabla w\dx &&= 0
\end{alignedat}
\tag{38}\end{equation} F ( c , μ , v ) G ( c , μ , w ) = ∫ Ω Δ t 1 ( c − c n ) v d x + ∫ Ω ( M ∇ μ ) ⋅ ∇ v d x = ∫ Ω [ μ − d c df ] w d x − ∫ Ω ( Λ ∇ c ) ⋅ ∇ w d x = 0 = 0 ( 38 ) which is a nonlinear coupled system of the form (29) .
We linearize the forms according to (30) and obtain
the following variations
D c F ⋅ Δ c = ∫ Ω 1 Δ t Δ c v d x D μ F ⋅ Δ μ = ∫ Ω ( M ∇ ( Δ μ ) ) ⋅ ∇ v d x D c G ⋅ Δ c = − ∫ Ω d 2 f d c 2 Δ c w d x − ∫ Ω ( Λ ∇ ( Δ c ) ) ⋅ ∇ w d x D μ G ⋅ Δ μ = ∫ Ω Δ μ w d x \begin{align*}
\varn{F}{c}{\Var c}
&= \int_\Omega \frac{1}{\Var t} \Var c\, v \dx \\
\varn{F}{\mu}{\Var \mu}
&= \int_\Omega (\BM\nabla (\Var\mu))\dtp \nabla v \dx \\
\varn{G}{c}{\Var c}
&= - \int_\Omega \dderiv{f}{c}\Var c \, w \dx
- \int_\Omega (\BLambda\nabla (\Var c))\dtp \nabla w\dx \\
\varn{G}{\mu}{\Var \mu}
&= \int_\Omega \Var\mu \, w \dx
\end{align*} D c F ⋅ Δ c D μ F ⋅ Δ μ D c G ⋅ Δ c D μ G ⋅ Δ μ = ∫ Ω Δ t 1 Δ c v d x = ∫ Ω ( M ∇ ( Δ μ )) ⋅ ∇ v d x = − ∫ Ω d c 2 d 2 f Δ c w d x − ∫ Ω ( Λ ∇ ( Δ c )) ⋅ ∇ w d x = ∫ Ω Δ μ w d x We substitute basis functions and obtain our system matrix
and vectors
P I J = ∫ Ω 1 Δ t N J N I d x Q I L = ∫ Ω ( M B L ) ⋅ B I d x r I = ∫ Ω 1 Δ t ( c ˉ − c n ) N I d x + ∫ Ω ( M ∇ μ ˉ ) ⋅ B I d x S K J = − ∫ Ω d 2 f d c 2 ∣ c = c ˉ N J N K d x − ∫ Ω ( Λ B J ) ⋅ B K d x T K L = ∫ Ω N L N K d x u K = ∫ Ω [ μ ˉ − d f d c ∣ c = c ˉ ] N K d x − ∫ Ω ( Λ ∇ c ˉ ) ⋅ B K d x \begin{align*}
P^{I\!J}
&= \int_\Omega \frac{1}{\Var t} N^JN^I \dx \\
Q^{IL}
&= \int_\Omega (\BM\BB^L)\dtp\BB^I \dx \\
r^{I}
&= \int_\Omega \frac{1}{\Var t}(\bar{c}-c_n)N^I \dx
+ \int_\Omega (\BM\nabla\bar{\mu})\dtp\BB^I\dx \\
S^{K\!J}
&= - \int_\Omega \dderiv{f}{c}\evat_{c=\bar{c}} N^J N^K \dx
- \int_\Omega (\BLambda \BB^J)\dtp \BB^K\dx \\
T^{K\!L}
&= \int_\Omega N^L N^K \dx \\
u^{K}
&= \int_\Omega \sbr{\bar{\mu} - \deriv{f}{c}\evat_{c=\bar{c}}} N^K \dx
- \int_\Omega (\BLambda\nabla \bar{c})\dtp \BB^K\dx
\end{align*} P I J Q I L r I S K J T K L u K = ∫ Ω Δ t 1 N J N I d x = ∫ Ω ( M B L ) ⋅ B I d x = ∫ Ω Δ t 1 ( c ˉ − c n ) N I d x + ∫ Ω ( M ∇ μ ˉ ) ⋅ B I d x = − ∫ Ω d c 2 d 2 f c = c ˉ N J N K d x − ∫ Ω ( Λ B J ) ⋅ B K d x = ∫ Ω N L N K d x = ∫ Ω [ μ ˉ − d c df c = c ˉ ] N K d x − ∫ Ω ( Λ ∇ c ˉ ) ⋅ B K d x which constitute the system
[ P Q S T ] [ Δ c Δ μ ] = [ r u ] \begin{equation}
\begin{bmatrix}
\BP & \BQ \\
\BS & \BT
\end{bmatrix}
\begin{bmatrix}
\Var \Bc \\ \Var \Bmu
\end{bmatrix}
=
\begin{bmatrix}
\Br \\ \Bu
\end{bmatrix}
\tag{39}\end{equation} [ P S Q T ] [ Δ c Δ μ ] = [ r u ] ( 39 ) Solution yields the update values Δ c \Var \Bc Δ c Δ μ \Var \Bmu Δ μ
c ← c ˉ + Δ c μ ← μ ˉ + Δ μ . \begin{equation}
\begin{alignedat}{3}
\Bc &\leftarrow \bar{\Bc} &&+ \Var\Bc \\
\Bmu &\leftarrow \bar{\Bmu} &&+ \Var\Bmu
.
\end{alignedat}
\tag{40}\end{equation} c μ ← c ˉ ← μ ˉ + Δ c + Δ μ . ( 40 ) The system is solved for c n + 1 c_{n+1} c n + 1 μ n + 1 \mu_{n+1} μ n + 1 t = t n t=t_n t = t n