[yueshen]

Hi everyone,

I am modifying the DPMFoam solver, and coupling it with a DEM software. The UEqn is compiled as shown below. The porosity n is obtained from the DEM software.

fvVectorMatrix UEqn
(
fvm::ddt(n, U) + fvm::div(nPhi, U)
- fvm::Sp(fvc::ddt(n) + fvc::div(nPhi), U)
- n * fvm::laplacian(nu, U)
);


I create the surfaceScalarField nPhi by equation “fvc::interpolate(*n_)* (*phi_)” and write it in each writing inverval. However, the output value of nPhi is zero for every timestep.

nPhi_ = new surfaceScalarField (IOobject ("nPhi", runTime_->timeName(), *mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE),
fvc::interpolate(*n_)* (*phi_));


Questions:
1. Do I need to update nPhi in the nNonOrthogonalCorrectors, as shown below? Because the output update of nPhi is zero so I have tried to add it in this way, but the solver becomes unable to converge, i.e., the pressure blows up.

if (piso_->finalNonOrthogonalIter())
{

nf = fvc::interpolate(n);
phi = phiHbyA - pEqn.flux()/nf;
nPhi = phi*nf;

U = HbyA
+ rAUc*fvc::reconstruct
(
(phicSUSu - pEqn.flux()/nf)/rAUcf
);
U.correctBoundaryConditions();

}

2. From the DPMFoam in OpenFoam, why DPMFoam does not update nf and nPhi?

if (pimple.finalNonOrthogonalIter())
{
phic = phiHbyA - pEqn.flux()/alphacf;

p.relax();

Uc = HbyA
+ rAUc*fvc::reconstruct((phicForces - pEqn.flux()/alphacf)/rAUcf);
Uc.correctBoundaryConditions();
}

3. Does the selections of solver in fvSvheme plays a significant role in convergence?
Thanks!

[yueshen]

I think I found the bug. nPhi is updated at the beginning of the DPMFoam.C.