Im trying to solve the following (transport & diffusion) nonlinear PDE numerically (via finite volume on a cuboid region. Some Material gets cooled down, s.t. in some areas the material becomes solid ($T<T_{sol}$) and in some other areas its still liquid ($T>T_{liq}$):

$0 = \lambda(T)\Delta T(x,t) + \nabla \cdot (u\rho h(T))$, where $u$ is a velocity vector. I guess the boundary conditions are not that important(on some sides dirichlet and on some other sides neumann; I could give more information, but i dont think thats necessary).

To solve that, i can linearize $h$ with a first order Taylor expansion: $h(T)=h(T_{(i-1)})+h'(T_{(i-1)})*(T-T_{(i-1)})$, plug it in and discretize to get a nonlinear system of equations: $A(T)T=b(T)$. To solve this system i use a fixpoint(picard) iteration with a very small $1 > \omega > 0$ relaxation parameter: Initial guess $T_{(0)}$, set $T_{(i)} = T_{(i-1)} + \omega (A(T_{(i-1)})b(T_{(i-1)}) - T_{(i-1)})$. Actually the $\omega$ is not fixed, its chosen in way such that the maximum entry in $|T_{(i)} - T_{(i-1)}|$ is decreasing in every step.

Now i have the nonlinear coefficients $\lambda(T_{(i-1)})$ and the coefficients $h'(T_{(i-1)})$ in $A(T_{(i-1)})$, and something like $-h(T_{(i-1)})-h'(T_{(i-1)})T_{(i-1)}$ in $b(T_{(i-1)})$.

Depending on the behaviour of $h$ (It's the specific enthalpy an depends on the modelled material), the method works fast, slow or doesnt converge at all. $h(T)$ is a continious function which is always strictly monotonically increasing, one could describe the function as follows:

$h'(T)=\alpha T$ for $T \notin [T_{sol},T_{liq}]$,

$h'(T) = \beta T$ for $T\in [T_{sol},T_{liq}]$, where $\beta >> \alpha$ and $[T_{sol},T_{liq}]$ is just a small pretty small interval. The smaller this Interval is, the higher is $\beta$. And if this Interval is too small and/or $\beta$ too large in comparism to $\alpha$, my method doesn't converge.

How can one explain this behaviour?

If i use another strategy to linearize the $h$:

$h(T)=h(T_0) + \frac{h(T_{(i-1)}) - h(T_0)}{T_{(i-1)} - T_0}(T-T_0)$, with some fixed $T_0 < T_{liq}$ , the method is more likely to converge, but it needs way more iterations.

Maybe one could think of a faster strategy or could explain to me when and why i can expect a convergence? I know that a fixpoint iteration converges if its contractive, but i guess thats impossible too proof in such cases.