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.

IntroductionEdit

Modeling of flocculated supensionsEdit

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 schemeEdit

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 OptimizationEdit

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.

WeblinksEdit

Last modified on 21 August 2007, at 20:55