View Revisions: Issue #3385

Summary 0003385: DPMFoam - "SemiImplicit" scheme for the source term gives unphysical fluid velocity distribution around the parcel
Revision 2019-11-13 17:17 by henry
Description Using DPMFoam to investigate a solid particle moving in fluid flow the "SemiImplicit" scheme for the momentum source gives unphysical velocity distribution around the parcel.
To show the phenomenon, a simplified simulation case is attached.
In the simulation a 2D flow flows from the left to the right in a channel. At the beginning one solid particle is positioned with zero velocity at the middle of the channel. Note that the parcel = particle in this case by the "nParticle 1" setting. Due to the small dimensions and fluid velocity, the flow will be laminar. As the particle diameter is also small the fluid drag will be the dominant force on it. This means that over the time the particle has to move with the same speed as the fluid in the channel.
Although the local cell fluid velocity and particle velocity will be nearly the same indeed, this velocity is usually smaller than the surrounding velocity of the fluid flow.
The phenomenon does not exist with the "explicit" momentum source term setting. However, the implicit scheme would be favorable, as the explicit scheme seemed to be unstable at higher particle volume fractions, and worked only with significantly decreased time steps.
The problem can be checked by running the attached simulation with the "SemiImplicit" and "explicit" schemes which can be set at the constant/kinematicCloudProperties file.
The difference between the two technique can be seen by the source code. In the followings I try to summarise it based on my understanding:

In DPMFoam.C the line
    fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
calls the SU() function of the particle cloud, which can be found in the KinematicCloudI.H. Here the return values of the two options (explicit/implicit) are:

       if (solution_.semiImplicit("U"))
        {
            const volScalarField::Internal
                Vdt(mesh_.V()*this->db().time().deltaT());

            return UTrans()/Vdt - fvm::Sp(UCoeff()/Vdt, U) + UCoeff()/Vdt*U;
        }
        else
        {
            tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
            fvVectorMatrix& fvm = tfvm.ref();

            fvm.source() = -UTrans()/(this->db().time().deltaT());

            return tfvm;
        }

where the volScalarField::Internal (DimensionedField) type dUTrans field is the momentum transferred from the parcels to the fluid in the given time step.
1. In the explicit scheme the force in each cell is explicitly given by dividing this momentum with the time step. Later in DPFoam.C it is further divided by the mesh volume to get the volumetric force density with the term "cloudVolSUSu". This term then will be included in the "phicForces" term and added to HbyA (see UEqn.H and pEqn.H). As in this case the "cloudSU" term contains only the explicit source term, which will be set to zero at DPMFoam.C, the cloudSU has no effect in the momentum equation.

2. In the semiImplicit case the second and the third term in the return line "should" cancel out each other in a steady case. I assume that these terms are given to increase the stability based on https://www.cfd-online.com/Forums/openfoam-solving/149682-dpmfoam-interpretation-uceqn-h.html
The return line is given with the type of volScalarField::Internal (DimensionedField). However the function return type is fvVectorMatrix, and at the conversion between the two at fvMatrix.C adds a multiplication with mesh.V() for the source. This is the reason that here the return terms are also divided by the volume.
Later in UcEqn.H the cloudSU() adds only the second (implicit) term to the momentum equation, while the phicForces contains only the first and the third (explicit) terms.
Revision 2019-11-12 15:51 by peth
Description Using DPMFoam to investigate a solid particle moving in fluid flow the "SemiImplicit" scheme for the momentum source gives unphysical velocity distribution around the parcel.
To show the phenomenon, a simplified simulation case is attached.
In the simulation a 2D flow flows from the left to the right in a channel. At the beginning one solid particle is positioned with zero velocity at the middle of the channel. Note that the parcel = particle in this case by the "nParticle 1" setting. Due to the small dimensions and fluid velocity, the flow will be laminar. As the particle diameter is also small the fluid drag will be the dominant force on it. This means that over the time the particle has to move with the same speed as the fluid in the channel.
Although the local cell fluid velocity and particle velocity will be nearly the same indeed, this velocity is usually smaller than the surrounding velocity of the fluid flow.
The phenomenon does not exist with the "explicit" momentum source term setting. However, the implicit scheme would be favorable, as the explicit scheme seemed to be unstable at higher particle volume fractions, and worked only with significantly decreased time steps.
The problem can be checked by running the attached simulation with the "SemiImplicit" and "explicit" schemes which can be set at the constant/kinematicCloudProperties file.
The difference between the two technique can be seen by the source code. In the followings I try to summarise it based on my understanding:

In DPMFoam.C the line
    fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
calls the SU() function of the particle cloud, which can be found in the KinematicCloudI.H. Here the return values of the two options (explicit/implicit) are:

       if (solution_.semiImplicit("U"))
        {
            const volScalarField::Internal
                Vdt(mesh_.V()*this->db().time().deltaT());

            return UTrans()/Vdt - fvm::Sp(UCoeff()/Vdt, U) + UCoeff()/Vdt*U;
        }
        else
        {
            tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
            fvVectorMatrix& fvm = tfvm.ref();

            fvm.source() = -UTrans()/(this->db().time().deltaT());

            return tfvm;
        }

where the volScalarField::Internal (DimensionedField) type dUTrans field is the momentum transferred from the parcels to the fluid in the given time step.
1. In the explicit scheme the force in each cell is explicitly given by dividing this momentum with the time step. Later in DPFoam.C it is further divided by the mesh volume to get the volumetric force density with the term "cloudVolSUSu". This term then will be included in the "phicForces" term and added to HbyA (see UEqn.H and pEqn.H). As in this case the "cloudSU" term contains only the explicit source term, which will be set to zero at DPMFoam.C, the cloudSU has no effect in the momentum equation.

2. In the semiImplicit case the second and the third term in the return line "should" cancel out each other in a steady case. I assume that these terms are given to increase the stability based on https://www.cfd-online.com/Forums/openfoam-solving/149682-dpmfoam-interpretation-uceqn-h.html
The return line is given with the type of volScalarField::Internal (DimensionedField). However the function return type is fvVectorMatrix, and at the conversion between the two at fvMatrix.C adds a multiplication with mesh.V() for the source. This is the reason that here the return terms are also divided by the volume.
Later in UcEqn.H the cloudSU() adds only the second (implicit) term to the momentum equation, while the phicForces contains only the first and the third (explicit) terms.