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 supensions

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

$\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = \frac{\partial A^2(u)}{\partial x^2}, \quad u(0,x)=u_0(), \quad 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) = \left \{ \begin{matrix} & u_{\infty} u (1-u)^C & \mbox{for } 0 \le u < u_{\max} \\ & 0 & \mbox{for } u \le 0 \mbox{ and } u \ge u_{\max} \end{matrix} \right.$

and the diffusive flux is given by

$A(u)=\int_0^{u} a(s) ds, \quad a(u)= a_0 (1-u)^C u^{k-1} ,$

which results from the insertion of the power law

$\sigma_{\mathrm{e}}(u) = \left \{ \begin{matrix} \sigma_0 \left ( \left ( \frac{u}{u_{\mathrm{c}}} \right )^k - 1 \right ) & \mbox{for } u_{\mathrm{c}}< u \\ 0 & \mbox{for } u \le u_{\mathrm{c}} \end{matrix} \right.$

into

$a(u)= \frac{f(u) \sigma_{\mathrm{e}}'(u)}{\Delta \varrho g u}, \quad a_0=\frac{u_\infty \sigma_0 k }{\Delta \varrho g u_{\mathrm{c}}^k} .$

In the closure, the constants $u_{\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{\Delta t}{\Delta x} \left ( f_{j+1/2}^{n+1} - f_{j-1/2}^{n+1} \right ) + \frac{\Delta t}{\Delta x^2} \left ( A(u_{j-1}^{n+1}) - 2 A(u_{j}^{n+1}) + A(u_{j+1}^{n+1}) \right ), \quad j = 2, \dots, J-1$

and at the bondaries ("boundary scheme") as

$u_1^{n+1}=u_1^n - \frac{\Delta t}{\Delta x} f_{j+1/2}^{n+1} + \frac{\Delta t}{\Delta x^2} \left ( A(u_{j+1}^{n+1}) - A(u_{j}^{n+1}) \right ), \quad u_J^{n+1}=u_J^n + \frac{\Delta t}{\Delta x} f_{J-1/2}^{n+1} - \frac{\Delta t}{\Delta x^2} \left ( A(u_{J}^{n+1}) - A(u_{J-1}^{n+1}) \right )$

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{EO}}(u,v) = \left \{ \begin{matrix} & f(u), & \mbox{for } u \le u_{\mathrm{m}}, v \le u_{\mathrm{m}}, \\ & f(u)+f(v)-f(u_{\mathrm{m}}), & \mbox{for } u \le u_{\mathrm{m}}, v > u_{\mathrm{m}}, \\ & f(u_{\mathrm{m}}), & \mbox{for } u > u_{\mathrm{m}}, v \le u_{\mathrm{m}}, \\ & f(v), & \mbox{for } u > u_{\mathrm{m}}, v > u_{\mathrm{m}} . \end{matrix} \right.$

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}), \quad 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}), \quad j = 1, \dots, J$

are inserted. Abbreviating the time evolution step as

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

the linearized marching formula for the interior scheme becomes

$\Delta u_j^n = - \frac{\Delta t}{\Delta x} \left ( f_{j+1/2}^n + \mathcal{J}_{j+1/2}^- \Delta u_j^n + \mathcal{J}_{j+1/2}^+ \Delta u_{j+1}^n + f_{j-1/2}^n + \mathcal{J}_{j-1/2}^- \Delta u_{j-1}^n + \mathcal{J}_{j-1/2}^+ \Delta u_j^n \right )$$+\frac{\Delta t}{\Delta x^2} \left ( A(u_{j-1}) - 2A(u_j) + A(u_{j+1}) + a(u_{j-1}) \Delta u_{j-1} -2 a(u_j) \Delta u_j + a(u_{j+1}) \Delta u_{j+1} \right ) ,$

where

$J_{j+1/2}^+ = \frac{\partial f(u_j^n, u_{j+1}^n)}{\partial u_{j+1}^n}, \quad 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

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

which is of the form

$m_{j,j-1} \Delta u_{j-1}^n + m_{j,j} \Delta u_j^n + m_{j,j+1} \Delta u_{j+1} = b_j ,$

or, in more compact notation,

$M \Delta u^n = b, \quad \begin{bmatrix} m_{11} & m_{12} & 0 & \cdots & \cdots & \cdots & \cdots & 0 \\ m_{21} & m_{22} & m_{23} & 0 & \cdots & \cdots & \cdots & 0 \\ \vdots & & & & & & & \vdots \\ \vdots & & & & & & \ddots & \vdots \\ 0 & \cdots & \cdots & \cdots & \cdots & 0 & m_{J,J-1} & m_{J,J} \end{bmatrix} , \quad b= \begin{bmatrix} b_1 \\ \vdots \\ b_J \end{bmatrix} , \quad \Delta u^n = \begin{bmatrix} \Delta u_1^n \\ \vdots \\ \Delta u_J^n \end{bmatrix} .$

Parameter identification as Optimization

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

$\min_e q(e), \quad q(e) = \| h(e) - H \|^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) consiting of two parameters. The optimization can be iteratively done by employing the Newton method as

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

where

$Q = \begin{bmatrix} \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{bmatrix} ,$

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

$0 \approx Q(e^*) = Q(e) + \nabla_e Q(e) \Delta e + \frac{1}{2} \Delta e^{\mathrm{T}} Q \Delta e, \quad \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