# Optimized Analytical-Numerical Procedure for Ultrasonic Sludge Treatment for Agricultural Use
## 6 The Numerical Model
### Algebraic Formulation  
For computational needs, we define the mass matrix and stiffness matrix for \(T\) as:  
\[(\mathbf{M}_{T})_{ij}=\int_{\Omega}\rho C_{p}\phi_{j}\phi_{i}\,d\mathbf{x}. \tag{28}\]  
and:  
\[(\mathbf{A}_{T})_{ij}=\int_{\Omega}k\
abla\phi_{j}\cdot\
abla\phi_{i}\,d \mathbf{x}+\int_{\Omega}\rho C_{p}(\mathbf{u}\cdot\
abla\phi_{j})\phi_{i}\,d \mathbf{x}. \tag{29}\]  
respectively, while, for \(W\), the mass and stiffness matrix (dependent on \(T\)), take the following forms:  
\[(\mathbf{M}_{W})_{ij}=\int_{\Omega}\phi_{j}\phi_{i}\,d\mathbf{x}, \tag{30}\]\[(\mathbf{A}_{W}(T))_{ij}=\int_{\Omega}k_{e}(T)\big{(}p_{v}(T)-p_{v,\text{amb}} \big{)}\
abla\phi_{j}\cdot\
abla\phi_{i}\,d\mathbf{x}. \tag{31}\]  
Finally, considering the source term:  
\[(\mathbf{F}_{T})_{i}=\int_{\Omega}Q\phi_{i}\,d\mathbf{x},\quad(\mathbf{F}_{W})_ {i}=\int_{\Omega}\phi_{i}W^{n}\,d\mathbf{x}, \tag{32}\]  
our coupled linear system becomes:  
\[\begin{cases}\mathbf{M}_{T}\frac{\mathbf{T}^{n+1}-\mathbf{T}^{n}}{\Delta t}+ \mathbf{A}_{T}\mathbf{T}^{n+1}=\mathbf{F}_{T}+\lambda\mathbf{M}_{T}\frac{ \mathbf{W}^{n+1}-\mathbf{W}^{n}}{\Delta t}\\ \mathbf{M}_{W}\frac{\mathbf{W}^{n+1}-\mathbf{W}^{n}}{\Delta t}-\mathbf{A}_{W} (\mathbf{T}^{n+1})\mathbf{W}^{n+1}=0.\end{cases} \tag{33}\]  
The resolution of the described coupled algebraic system requires numerical methods since it includes nonlinear dependencies and couplings between the variables \(\mathbf{T}^{n+1}\) and \(\mathbf{W}^{n+1}\). The system (33) can be written in the following form:  
\[\begin{bmatrix}\mathbf{M}_{T}/\Delta t+\mathbf{A}_{T}&-\lambda \mathbf{M}_{T}/\Delta t\\ 0&\mathbf{M}_{W}/\Delta t-\mathbf{A}_{W}(\mathbf{T}^{n+1})\end{bmatrix}\begin{bmatrix} \mathbf{T}^{n+1}\\ \mathbf{W}^{n+1}\end{bmatrix}\] \[=\begin{bmatrix}\mathbf{F}_{T}+\mathbf{M}_{T}\mathbf{T}^{n}/ \Delta t+\lambda\mathbf{M}_{T}\mathbf{W}^{n}/\Delta t\\ \mathbf{M}_{W}\mathbf{W}^{n}/\Delta t\end{bmatrix}. \tag{34}\]  
It is nonlinear, since \(\mathbf{A}_{W}(\mathbf{T}^{n+1})\) depends on \(\mathbf{T}^{n+1}\)). Then, using an iterative approach based on Newton-Raphson, we define the residual as:  
\[\mathbf{R}(\mathbf{T},\mathbf{W})=\begin{bmatrix}\mathbf{M}_{T}\frac{\mathbf{T }-\mathbf{T}^{n}}{\Delta t}+\mathbf{A}_{T}\mathbf{T}-\lambda\mathbf{M}_{T} \frac{\mathbf{W}-\mathbf{W}^{n}}{\Delta t}-\mathbf{F}_{T}\\ \mathbf{M}_{W}\frac{\mathbf{W}-\mathbf{W}^{n}}{\Delta t}-\mathbf{A}_{W}( \mathbf{T})\mathbf{W}\end{bmatrix}, \tag{35}\]  
such that we find \((\mathbf{T}^{n+1},\mathbf{W}^{n+1})\) such that \(\mathbf{R}(\mathbf{T}^{n+1},\mathbf{W}^{n+1})=0\). Then, we expand \(\mathbf{R}(\mathbf{T},\mathbf{W})\) around an initial estimate \((\mathbf{T}_{k},\mathbf{W}_{k})\) at step \(k\) using a Taylor series:  
\[\mathbf{R}(\mathbf{T}_{k+1},\mathbf{W}_{k+1})\approx\mathbf{R}(\mathbf{T}_{k},\mathbf{W}_{k})+\frac{\partial\mathbf{R}}{\partial\mathbf{T}}\Delta\mathbf{T }+\frac{\partial\mathbf{R}}{\partial\mathbf{W}}\Delta\mathbf{W}, \tag{36}\]  
where:  
\[\Delta\mathbf{T}=\mathbf{T}_{k+1}-\mathbf{T}_{k},\quad\Delta\mathbf{W}= \mathbf{W}_{k+1}-\mathbf{W}_{k}. \tag{37}\]