[stefsolo]

Hello all!

I'm simulating two-phase (air-water) flow using the interFoam solver to model the filling of a microfluidic channel connected to a microchamber.

The microfluidic device material is permeable to air but impermeable to water. The device is bonded to a non-permeable surface, therefore I must take into account the permeability only for the surfaces of the device.

To achieve this I have implemented a custom boundary condition on the air-permeable surface, using a fixedCodedValue that calculates the normal velocity to the permeable surface using the Darcy law, based on the internal and external pressure, to calculate the velocity with which air escapes the device. Below you can see my fixedCodedValue. Note that the velocity calculation is only active when the volume fraction of water, α, is less than 0.01, therefore only when air is present.

wall_pdms
{
type            codedFixedValue;
value           uniform (0 0 0);
name            darcyporouswall1;
codeInclude
#{
#include "fvCFD.H"
#include "alphaFixedPressureFvPatchScalarField.H"
#};
codeOptions
#{
-I$(WM_PROJECT_DIR)/src/finiteVolume/lnInclude \
-I$(WM_PROJECT_DIR)/src/meshTools/lnInclude \
-I$(WM_PROJECT_DIR)/src/transportModels/twoPhaseProperties/lnInclude
#};
code
#{
// Constants
const scalar mu = 1.813e-5; // Dynamic viscosity (Pa·s)
const scalar K = 1.0e-16; // Permeability (m^2)
const scalar Pout = 0; // Outlet pressure (Pa)
const scalar L = 1e-3; // Characteristic length (m)
// Access to the current fvPatch
const fvPatch& fvP = this->patch();
// Access to the face normal vectors
const vectorField& n = fvP.nf();
// Access to the entire pressure and alpha.water fields
const volScalarField& p = this->db().lookupObject<volScalarField>("p");
const volScalarField& alpha = this->db().lookupObject<volScalarField>("alpha.water");
// Access the patch values of pressure and alpha
const fvPatchScalarField& pPatch = p.boundaryField()[fvP.index()];
const fvPatchScalarField& alphaPatch = alpha.boundaryField()[fvP.index()];
// Loop through each face on the patch to set the velocity
forAll(fvP, faceI)
{
if (alphaPatch[faceI] < 0.01)
{
// Apply Darcy's Law only when alpha < 0.01
scalar Un = (K / mu) * ((pPatch[faceI] - Pout) / L);
this->operator[](faceI) = Un * n[faceI];
}
else
{
// Otherwise, set velocity to zero
this->operator[](faceI) = vector(0, 0, 0);
}
}
#};
   }

What I have tried:

1. For low K values (less than 1e-16) (so low velocity values) simulation proceeds. When I try to increase the K values and therefore the velocity the Courant number explodes and the simulation stops. The higher the value the faster the simulation stops.
   
2. When completely removing the codedValue for another type (noSlip) the simulation also proceeds.
   
3. For low problematic K values (10^-15, 10^-14) the simulation seems to stops at the same point of where the fluid reaches a specific place in the device where it encounters the second sharp edge of the microchamber inlet. I have tried to refine around that edge using a refinementBox but the problem persisted.
   
4. I tried to fillet the sharp edges of the inlet of the microchamber but still had the same problem.
   
5. I removed the permeability of the side walls and kept it active only at the top surface. Still I had the same results.
   
6. Regarding the mesh I tried to increase the cell number at the thin direction from 1 --> 5 --> 8, but still had the same problems
   
7. I tried to decrease the maxCo from 1 to 0.5. No change.
   
8. I tried to change div(rhoPhi,U) Gauss linearUpwind grad(U); => div(rhoPhi,U) bounded Gauss linearUpwind grad(U); and div(phirb,alpha) Gauss linear; => div(phirb,alpha) bounded Gauss linear;, still the same results.

All the above seem to happen no matter the mesh (coarse or fine).

Please note that the thickness of the device is very small (0.015e-3 m).

Please note that my timeStep is adjustable and despite for it becoming extremely small loosing all its physical meaning Courant still remains high, that is why the simulation stops.

I am also attaching my U, p_rgh, fvSchemes and fvSolutions files as well as a screenshot of the geometry.

I am using openfoam2506,

Any suggestions to fix the problem will be very much appreciated.

Thanks in advance.