Convergence Problem simpleFoam with reaction
 

Dear Foamers,

I hope this is the correct Forum for my Problem...

I'm trying to implement reaction into simpleFoam to simulate the steady-state of a radical polymerisation.

As a first test-setup I simplify my case:
- the case is laminar and incompressible
- density and viscosity of the fluid do not change with reaction
- no diffusion
- convection happens at constant velocity
- 2D-flow with zeroGradient walls --> species concentration changes 1D

In the simple-loop I first calculate the reaction rates of the various reactions and then the source term for each species. The component balance is then solved in YEqn.H (according to reactingFoam). However, in my solver, Yi denotes the species concentration in [mole/m³] instead of the mass fraction.
Calculations of UEqn.H and pEqn.H were not changed.

I started with simplified reactions:
Initiator --> 2*Radical //Initiator decomposition
Radical + Monomer --> Radical //Propagation
Radical + Radical --> Dead chain //Termination
The chain length of the growing radical is not important, since I'm doing all my calculations with concentration and not mass fractions.

the rate constants are first assumed constant, they are in the correct order of magnitude:
k_d = 7.8e-4 1/s
k_p = 27 m³/mol/s
k_t = 3e4 m³/mol/s

the Source terms are:
r_ini = -k_d * c(ini);
r_rad = 2*k_d*c(ini) - 2*k_t*c(rad)²
r_mon = -k_p*c(mon)*c(rad)

Without termination (if I set k_t=0), the solver works fine, and I checked the initiator concentration is correct (I have a predici curve for comparison) so there seems to be no major (theoretical) fault in the set-up.

However, if I turn on the termination reaction, the solution does not converge.
Here are the things I tried so far:

- reducing the time step
- refining the mesh
- change the convective Schemes for Yi (even upwind did not work)

I am grateful for any help or hints to stabilize my solver!
Thank you!

I attached the plots of the initial residuals for both cases (with and without termination).

I have now gained further understanding of my problem:

With termination, the r_rad reaction becomes very small (in the order of magnitude of -10e20) which leads to unphysical negative concentrations of this species. Since the Monomer and radical are coupled, the monomer does also not converge.

I defined the reaction rate of the radical species as follows (without units):

r_rad = 2* 7.8e-4 * c(ini) - 2 * 3e4 *c(rad)²

The first term is always in the order of magnitude of 1e-3 and c(ini) converges.
The second part however is sensitive to changes of rad because of the square.

So I'm thinking, some way of "dampening" the solution progress could help, what do you think?
How can I do that?