View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0003385OpenFOAMBugpublic2019-11-12 15:512020-08-05 10:39
Reporterpeth Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
Platformx86-64OSArch LinuxOS Version(please specify)
Product Versiondev 
Fixed in Versiondev 
Summary0003385: DPMFoam - "SemiImplicit" scheme for the source term gives unphysical fluid velocity distribution around the parcel
DescriptionUsing 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.
Steps To Reproduce1. Download and extract the attached simulation folder.
2. Generate the mesh with the ./mesh_gen script.
3. Run the simulation with the ./run script.
4. Check the fluid velocity cell values around the moving parcel.
Additional InformationThe simulation was run with OpenFOAM-5.x (manually compiled). The corresponding code for the source term handling in OpenFOAM-7 looks similar, so I assume that the problem will also exist there.
TagsNo tags attached.

Relationships

related to 0003464 closedwill Unstable packing model in MPPICFoam without the isotropy model 

Activities

peth

2019-11-12 15:51

reporter  

dpmfoam_source_report.tar.gz (5,012 bytes)
magU.png (16,037 bytes)   
magU.png (16,037 bytes)   

henry

2019-11-13 11:18

manager   ~0010893

Have you tried running with outer correctors to update the semi-implicit sources? e.g. set

   nOuterCorrectors 5;

for example.

henry

2019-11-13 11:22

manager   ~0010894

I tried to run the case but it fails at start-up:

--> FOAM FATAL IO ERROR:
keyword alpha.water is undefined in dictionary

Could you provide an updated version for OpenFOAM-dev?

peth

2019-11-13 16:08

reporter   ~0010898

I have compiled the current version of OpenFOAM-dev from GitHub, and corrected the simulation folder accordingly. The new simulation is attached as "dpmfoam_source_report_dev.tar.gz".
Initially the nOuterCorrectors was 1 and now I have checked it with 5 and 10. The fluid velocity is decreased around the parcel in all cases.
dpmfoam_source_report_dev.tar.gz (4,909 bytes)

henry

2019-11-13 17:01

manager   ~0010899

Thanks, I will investigate further.

henry

2019-11-13 17:17

manager   ~0010900

Try this:

commit b22545b9315cd72034905474b28e88d90b012583 (HEAD -> master, origin/master, origin/HEAD)
Author: Henry Weller <http://openfoam.org>
Date: Wed Nov 13 17:15:12 2019 +0000

    DPMFoam: Changed the cloud source splitting to handle symmetric semi-implicit sources consistently
    
    Resolves bug-report https://bugs.openfoam.org/view.php?id=3385

peth

2019-11-14 19:37

reporter   ~0010902

Thanks for the correction, I have tried it. The velocity field looks much more smooth now, see the attached pictures with some injected particles from the inlet with the original and the corrected solver. (The particles are magnified compared to their real sizes just for the visualization).
If I interpret the modification well, as
    cloudSU.diag()= - Ucoeff/dt
the modified code
     cloudVolSUSu.primitiveFieldRef() = (cloudSU.diag()*Uc() - cloudSU.source())/mesh.V();
means that the explicit (Ucoeff/Vdt * Uc_old) term will be added by the second part (.source()) and also subtracted by the first part (.diag()), so finally in cloudVolSUsu only the (UTrans/Vdt) part remains. I think this term is much smaller than the (Ucoeff/Vdt * Uc_old) part (and it should be almost zero) in the investigated case.
Then the line
    cloudSU.source() = cloudSU.diag()*Uc();
results that in UcEqn besides the implicit source term -fvm::Sp(UCoeff/Vdt, Uc_new) the similar, but explicit source (Ucoeff/Vdt * Uc_old) is also included.
Can you please tell me why this method gives better result than the previous?
mag_Uf_injected_particles_orig.png (127,827 bytes)   
mag_Uf_injected_particles_orig.png (127,827 bytes)   
mag_Uf_injected_particles_corrected.png (102,847 bytes)   
mag_Uf_injected_particles_corrected.png (102,847 bytes)   

henry

2019-11-14 19:57

manager   ~0010903

This approach avoids the cell-face splitting error and ensures that the cell parts and the face parts of the source are consistent.
How do the results of the explicit and semi-implicit approach compare for this case?

peth

2019-11-16 13:27

reporter   ~0010915

I have run the simulation with the corrected code and the explicit scheme. The velocity field and the particle distribution looks similar to the previous one with the semiImplicit scheme, see on the attached picture. The fluid velocity profiles have been compared at x=L/4 (where L is the whole channel length), and they are almost the same (it is shown at the attached pdf).
The two method anyway can provide different results, e.g. when a lot of particles are injected at t=0 to a specific part of the channel with zero velocity (this simulation is also attached). At the beginning the particles have to be accelerated up, which results - for a short time - a locally decreased velocity field. The semiImplicit technique is able to capture this phenomenon, while the explicit (with the given time step size) does not.
mag_Uf_particle_inject_corrected_explicit.png (105,747 bytes)   
mag_Uf_particle_inject_corrected_explicit.png (105,747 bytes)   
velocity_profile_compare.pdf (9,522 bytes)
dpmfoam_source_report_dev_centrecloud.tar.gz (11,243 bytes)

henry

2019-11-18 08:06

manager   ~0010919

Resolved by commit b22545b9315cd72034905474b28e88d90b012583

will

2020-07-31 09:49

manager   ~0011442

Note that commit https://github.com/OpenFOAM/OpenFOAM-dev/commit/43d66b5e7c566e087db187df023f134422ee8a4b has reverted this change by default. The change to the splitting caused significant problems with MPPIC simulations. There is now a keyword in "system/fvSolution.PIMPLE" called "cloudMomentumSplit" which can be set to the following values:

    faceExplicitCellImplicit
    faceExplicitCellLagged
    faceImplicit

The first is the default and causes the artefacts in the velocity field detailed in this report. The second is the fix from commit b22545b9 which removes the artefacts but causes issues in MPPIC simulations. The third is another approach which removes the artefacts and is less troublesome for most MPPIC cases, but the cyclone tutorial still doesn't run. See the following file for more information:

     https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/applications/solvers/lagrangian/denseParticleFoam/denseParticleFoam.C

These choices may be removed in the future if MPPIC is developed to the point of being tolerant to these sort of details. We may also judge that the default needs to be different, so if this is important to you then I recommend you set the control explicitly.

So, this report is still considered fixed, it's just that an additional entry (faceExplicitCellLagged or faceImplicit) is now needed to resolve the issue.

Issue History

Date Modified Username Field Change
2019-11-12 15:51 peth New Issue
2019-11-12 15:51 peth File Added: dpmfoam_source_report.tar.gz
2019-11-12 15:51 peth File Added: magU.png
2019-11-13 11:18 henry Note Added: 0010893
2019-11-13 11:22 henry Note Added: 0010894
2019-11-13 16:08 peth File Added: dpmfoam_source_report_dev.tar.gz
2019-11-13 16:08 peth Note Added: 0010898
2019-11-13 17:01 henry Note Added: 0010899
2019-11-13 17:17 henry Note Added: 0010900
2019-11-13 17:17 henry Product Version => dev
2019-11-13 17:17 henry Description Updated
2019-11-13 17:17 henry Steps to Reproduce Updated
2019-11-14 19:37 peth File Added: mag_Uf_injected_particles_orig.png
2019-11-14 19:37 peth File Added: mag_Uf_injected_particles_corrected.png
2019-11-14 19:37 peth Note Added: 0010902
2019-11-14 19:57 henry Note Added: 0010903
2019-11-16 13:27 peth File Added: mag_Uf_particle_inject_corrected_explicit.png
2019-11-16 13:27 peth File Added: velocity_profile_compare.pdf
2019-11-16 13:27 peth File Added: dpmfoam_source_report_dev_centrecloud.tar.gz
2019-11-16 13:27 peth Note Added: 0010915
2019-11-18 08:06 henry Assigned To => henry
2019-11-18 08:06 henry Status new => resolved
2019-11-18 08:06 henry Resolution open => fixed
2019-11-18 08:06 henry Fixed in Version => dev
2019-11-18 08:06 henry Note Added: 0010919
2020-07-31 09:49 will Note Added: 0011442
2020-08-05 10:39 will Relationship added related to 0003464