Sedimentation/Parameter Identification

< Sedimentation

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.

Contents

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