Sedimentation/Parameter Identification

Even though it is tried to keep this chapter on Parameter Identification of Flocculated Suspensions as self-comprehensive as possible, preliminar knowledge on numerical Methods and the Modeling of suspensions are useful. In particular, the Newton-Raphson scheme to solve nonlinear systems of equations for the optimization and Finite-Volume-Methods for the solution of partial differential equations are applied.

Modeling of flocculated suspensions

The batch settling process of flocculated suspensions is modeled as an initial value problem

$\frac{\partial u}{\partial t}+\frac{\partial f(u)}{\partial x}=\frac{\partial A^2(u)}{\partial x^2},u(0,x)=u_0(),x\in [0,L]$

where $u$ denotes the volume fraction of the dispersed solids phase. For the closure, the convective flux function is given by the Kynch batch settling function with Richardson-Zaki hindrance function

$f(u)=\{\begin{array}{ccc}& u_{\mathrm{\infty }}u(1-u{)}^C& \text{for }0\le u<u_{\max }\\ & 0& \text{for }u\le 0\text{ and }u\ge u_{\max }\end{array}$

and the diffusive flux is given by

$A(u)=\underset{0}{\overset{u}{\int }}a(s)ds,a(u)=a_0(1-u{)}^C{u}^{k-1},$

which results from the insertion of the power law

${\sigma }_{\mathrm{e}}(u)=\{\begin{array}{cc}{\sigma }_0({\left(\frac{u}{u_{\mathrm{c}}}\right)}^k-1)& \text{for }{u}_{\mathrm{c}}<u\\ 0& \text{for }u\le u_{\mathrm{c}}\end{array}$

into

$a(u)=\frac{f(u){\sigma }_{{\mathrm{e}}^{\prime }}(u)}{\mathrm{\Delta }\varrho gu},a_0=\frac{u_{\mathrm{\infty }}{\sigma }_{0}k}{\mathrm{\Delta }\varrho g{u}_{\mathrm{c}}^k}.$

In the closure, the constants $u_{\mathrm{\infty }},C,a_0,u_{\max },k$ are partly known.

Numerical scheme

The numerical scheme for the solution of the direct problem is written in conservative form as a marching formula for the interior points ("interior scheme") as $u_j^{n+1}=u_j^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}(f_{j+1/2}^{n+1}-f_{j-1/2}^{n+1})+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}(A(u_{j-1}^{n+1})-2A(u_j^{n+1})+A(u_{j+1}^{n+1})),j=2,\dots ,J-1$

and at the boundaries ("boundary scheme") as

$u_1^{n+1}=u_1^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{f}_{j+1/2}^{n+1}+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}(A(u_{j+1}^{n+1})-A(u_j^{n+1})),u_J^{n+1}=u_J^n+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{f}_{J-1/2}^{n+1}-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}(A(u_J^{n+1})-A(u_{J-1}^{n+1}))$

For a first-order scheme, the numerical flux function becomes

$f_{j+1/2}^n=f(u_j^n,u_{j+1}^n).$ If the flux function has one single maximum, denoted by

$u_{\mathrm{m}}$, the Engquist-Osher numerical flux function can be stated as

$f^{\mathrm{E}\mathrm{O}}(u,v)=\{\begin{array}{ccc}& f(u),& \text{for }u\le u_{\mathrm{m}},v\le u_{\mathrm{m}},\\ & f(u)+f(v)-f(u_{\mathrm{m}}),& \text{for }u\le u_{\mathrm{m}},v>u_{\mathrm{m}},\\ & f(u_{\mathrm{m}}),& \text{for }u>u_{\mathrm{m}},v\le u_{\mathrm{m}},\\ & f(v),& \text{for }u>u_{\mathrm{m}},v>u_{\mathrm{m}}.\end{array}$

For linearization, the Taylor formulae

$f(u_j^{n+1},u_{j+1}^{n+1})=f(u_j^n,u_{j+1}^n)+\frac{\partial f(u_j^n,u_{j+1})}{\partial u_j^n}(u_j^{n+1}-u_j^n)+\frac{\partial f(u_j^n,u_{j+1})}{\partial u_{j+1}^n}(u_{j+1}^{n+1}-u_{j+1}^n),j=1,\dots ,J$

and

$A(u_j^{n+1})=A(u_j^n)+\frac{\partial A(u_j^n)}{\partial u_j^n}(u_j^{n+1}-u_j^n),j=1,\dots ,J$

are inserted. Abbreviating the time evolution step as

$\mathrm{\Delta }{u}_j^n=u_j^{n+1}-u_j^n,j=1,\dots ,J,$

the linearized marching formula for the interior scheme becomes

$\mathrm{\Delta }{u}_j^n=-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}(f_{j+1/2}^n+{𝒥}_{j+1/2}^{-}\mathrm{\Delta }{u}_j^n+{𝒥}_{j+1/2}^{+}\mathrm{\Delta }{u}_{j+1}^n+f_{j-1/2}^n+{𝒥}_{j-1/2}^{-}\mathrm{\Delta }{u}_{j-1}^n+{𝒥}_{j-1/2}^{+}\mathrm{\Delta }{u}_j^n)$ $+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}(A(u_{j-1})-2A(u_j)+A(u_{j+1})+a(u_{j-1})\mathrm{\Delta }{u}_{j-1}-2a(u_j)\mathrm{\Delta }{u}_j+a(u_{j+1})\mathrm{\Delta }{u}_{j+1}),$

where

$J_{j+1/2}^{+}=\frac{\partial f(u_j^n,u_{j+1}^n)}{\partial u_{j+1}^n},J_{j+1/2}^{-}=\frac{\partial f(u_j^n,u_{j+1}^n)}{\partial u_j^n}.$

Rearrangement leads to a block-triangular linear system

$-(\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}({𝒥}_{j-1/2}^{-}+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}a(u_{j-1}^n))\mathrm{\Delta }{u}_{j-1}^n+(I+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}({𝒥}_{j+1/2}^{-}-{𝒥}_{j-1/2}^{+})+\frac{2\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}a(u_j^n))\mathrm{\Delta }{u}_j^n$ $+(\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}{𝒥}_{j+1/2}^{+}-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}a(u_{j+1}^n))\mathrm{\Delta }{u}_{j+1}^n$ $=\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}(f(u_j^n,u_{j+1}^n)-f(u_{j-1}^n,u_j^n))-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}(A(u_{j-1}^n)-2A(u_j^n)+A(u_{j+1}^n)),j=2,\dots ,J-2$

which is of the form

$m_{j,j-1}\mathrm{\Delta }{u}_{j-1}^n+m_{j,j}\mathrm{\Delta }{u}_j^n+m_{j,j+1}\mathrm{\Delta }{u}_{j+1}=b_j,$

or, in more compact notation,

$M\mathrm{\Delta }{u}^n=b,\left[\begin{array}{cccccccc}{m}_{11}& m_{12}& 0& \cdots & \cdots & \cdots & \cdots & 0\\ m_{21}& m_{22}& m_{23}& 0& \cdots & \cdots & \cdots & 0\\ ⋮& & & & & & & ⋮\\ ⋮& & & & & & \ddots & ⋮\\ 0& \cdots & \cdots & \cdots & \cdots & 0& m_{J,J-1}& m_{J,J}\end{array}\right],b=\left[\begin{array}{c}{b}_1\\ ⋮\\ b_J\end{array}\right],\mathrm{\Delta }{u}^n=\left[\begin{array}{c}\mathrm{\Delta }{u}_1^n\\ ⋮\\ \mathrm{\Delta }{u}_J^n\end{array}\right].$

Parameter identification as Optimization

The goal of the parameter identification by optimization is to minimize the cost function over the parameter space

$\underset{e}{\min }q(e),q(e)=\Vert h(e)-H{\Vert }^2$,

where h(e) denotes the interface that is computed by simulations and H is the measured interface. Without loss of generality we consider a parameter set e=(e_1, e_2) consisting of two parameters. The optimization can be iteratively done by employing the Newton method as

$Q(e_k)\mathrm{\Delta }{e}_k={\mathrm{\nabla }}_{e}q(e_k),e\subset \{C,v_{\mathrm{\infty }},a_0,u_{\mathrm{c}},k\},$

where

$Q=\left[\begin{array}{cc}\frac{{\partial }^{2}q}{\partial e_1^2}& \frac{{\partial }^{2}q}{\partial e_1{e}_2}\\ \frac{{\partial }^{2}q}{\partial e_2{e}_1}& \frac{{\partial }^{2}q}{\partial e_2^2}\end{array}\right],$

is the Hessian of q. The Newton method is motivated by the Taylor expansion

$0\approx Q(e^{*})=Q(e)+{\mathrm{\nabla }}_{e}Q(e)\mathrm{\Delta }e+\frac{1}{2}\mathrm{\Delta }{e}^{\mathrm{T}}Q\mathrm{\Delta }e,\mathrm{\Delta }e=e^{*}-e,$

where $e^{*}$ is the optimal parameter choice.

Weblinks

• http://de.wikipedia.org/wiki/Wikipedia:TeX Wikipedia:TeX
• http://en.wikipedia.org/wiki/Matlab