[jennyjian]

Hi guys!
I am working on a scalar transport equation in OpenFOAM, which is a little bit complex to make it converge. The equation is attached in the attachment equation 1:


where alpha, Df, Db, Hb have constant values in the whole fluid domain, while pO2 is the unknown oxygen partial pressure and SO2 is a function of pO2.

Firstly, without any treatment, the terms with Hb SO2 are treated as source terms:
The code is
fvm::div(phi,pO2) - fvm::laplacian(Df,pO2) ==
- fvc::SuSp(Hb/alpha*fvc::div(phi,SO2)/pO2,pO2)
+ fvc::Su(1/alpha*fvc::laplacian(Db,SO2))

Since the term Hb/alpha*fvc::div(phi,SO2) is about 7 times the value of div(phi,pO2). This solver diverges very quickly.

Later to solve this, I try to transform the term fvc::div(phi, SO2)
div(u) = 0 when the fluid equation converges. (incompressible)
The equation is listed in the attachment (equation2)
The code is
fvm::div(phi,pO2)-1.0/(1+Hb/alpha*inter1)*fvm::laplacian(Df,pO2)
==
1.0/(1+Hb/alpha*inter1)*fvc::laplacian(Drbc*Hb/alpha,SO2)
inter1 is grad(SO2)/grad(PO2)

To verify this solver, I use the steady flow across the cylinder Re = 5 as a case, where pO2 at the surface of the cylinder is 760 mmHg, and the inflow is 56mmHg.
However, the result (in attachment result and the correct result in literature Re5_full) seems the scalar pO2 is transported to only the first layer meshes near the boundary, which is a sign of convection only, without the diffusion. I wonder if the term 1.0/(1+Hb/alpha*inter1) can't be written outside the fvm::laplacian, is there any idea or suggestion on how to write this equation in OpenFOAM?

Thanks a lot!

[jennyjian]

The first figure is uploaded again