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}, \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 boundaries ("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) consisting 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