[ziemowitzima]

Dear OFoamers,

I am working on some problem with time dependent domain. Namely one of the borders of the domain is changing with time.
I was able to make my mesh time dependent and it works fine, but I encountered some other difficulty which I cant understand/solve.

Namely, after mesh is changed, solution behaves wrong. Some energy is added to the system. It happens only during first change of the domain. Consecutive changes do not show this strange behavior any more (solution behaves as it should)

You ca see it on the animation here:
http://youtu.be/wIexKK2vQ-8

from time 0 - 0.5 everything is ok
for t =0.5 there is change of the domain geometry
from time 0.5 - 1 we can see this "strange" behaviour, some extra energy is added and blue patch is moving (in fact it should not)
at time t = 1 there is change of the geometry again, but now (and later as well) everything seems to be ok (no extra/phantom energy is added).

Summing up, there is some problem after first change of the mesh.

Equation I am solving is just:

fvScalarMatrix omEqn
        (
        dimDt*fvm::ddt(om)
        + dimG*fvm::div(phi, om)    
      - fvm::laplacian(dimensionedScalar("1",dimensionSet(0, 2, 0, 0, 0),1), om)      
        );        
        omEqn.solve();

but the same happens as well for

fvScalarMatrix omEqn
        (
        dimDt*fvm::ddt(om)
      - fvm::laplacian(dimensionedScalar("1",dimensionSet(0, 2, 0, 0, 0),1), om)      
        );        
        omEqn.solve();

Do you have any idea why it is like that ?
Thanks
ZMM

[wyldckat]

Hi ZMM,

Disclaimer: I'm not an expert on this! I'll share as much as I know, but keep in mind that I might have some terminology confused and/or I might be mistaken .


A few possibilities come to mind:

• The first transition might give you an imperfect mesh, which induces the strange value. Run checkMesh with full check:
  

  checkMesh -allGeometry -allTopology

  Check for all of the geometry changes.
• Flux correction might be necessary, in case you aren't using it already.
• Mapping values might be incomplete, such as surface and point fields might not being mapped.

But there are at least a few details you haven't mentioned, which would make it easier to diagnose:

1. Are you running in serial/sequential mode or in parallel?
2. How exactly is your mesh being modified?
   1. Do you stop the solver, re-mesh, use mapFields and then launch it again?
   2. Or is your solver adapted to use dynamic meshing?
      1. If you are using adaptive meshing, what options are you using in "dynamicMeshDict"?
3. In which solver did you base your code on? Because each solver has its own tricks for dynamic meshes...

edit: after viewing the video a few times... it might be strange, but the solution you've got seems physical enough to me! Because you had large step change in the geometry, which could easily bump anything away from the wall !

Best regards,
Bruno

[ziemowitzima]

Hi Bruno,
checkMesh returns mesh OK.

It is simple 2D geometry, as you can see on the movie, so there is no troubles with topology.

I am not sure about flux correction, I added some commands (after pimpleDyMFoam), but it seems that they do not affect results. But truly speaking I do not understand this issue to much. I have the same strange behavior even for simple heat equations, where velocity is not needed.

Answering your questions:
1. I am running in serial model.

2. My mesh is modified by changing one of the boundaries during time integration.

My solver is based on icoStructFoam. Operation on mesh looks as follows:

forAll(fluidMesh,fluidI) {
            disp[fluidI].component(1) = 0.0;
                 } 
     fvc::makeAbsolute(phi, U);    
for (int ii=0; ii<=2; ii++)
     {
        volVectorField::GeometricBoundaryField &meshDisplacement=
        const_cast<volVectorField&>(mesh.objectRegistry::lookupObject<volVectorField>("cellDisplacement")).boundaryField();
        vectorField &mDisp=refCast<vectorField>(meshDisplacement[ss]);
        
        forAll(fluidMesh,fluidI) {                                                     
            vector here=fluidMesh.faceCentres()[fluidI];
            vector old = oldPoints[fluidI];                  
            disp[fluidI].component(1) = disp[fluidI].component(1)-0.1*(old.component(1) - ds2[fluidI]);
            //Info<<here.component(1)<<" "<< ds[fluidI]<<" "<<fabs(here.component(1) - ds[fluidI])<<endl;
            disp[fluidI].component(0) = 0.0;
            disp[fluidI].component(2) = 0.0;
            // vector neu = neu + disp;
            vector move=disp[fluidI];

            mDisp[fluidI]=move;                       
            }
            
        mesh.movePoints(motionPtr->newPoints());
        mesh.update();
        #include "correctPhi_cuette.H"
        fvc::makeRelative(phi, U);
        
    }

Which just move one of the boundary domain. To avoid to large changes I am repeating this operation several time with the fraction of the desired final state. It works fine, namely mesh is modified the way I need.

commands:

 fvc::makeAbsolute(phi, U);
        #include "correctPhi_cuette.H"
        fvc::makeRelative(phi, U);

I added recently, they do not change anything in the solution.

2. I am using dynamic meshing. dynamicMeshDict looks as follows:

twoDMotion      yes;

dynamicFvMesh dynamicMotionSolverFvMesh;

motionSolverLibs ("libfvMotionSolvers.so");

solver displacementLaplacian;

diffusivity  uniform;

3. I based it on icoStructFoam

edit: after viewing the video a few times... it might be strange, but the solution you've got seems physical enough to me! Because you had large step change in the geometry, which could easily bump anything away from the wall !

I does not because:
1. in case of just diffusion equation, it should get weaker, not stronger...
2. If I re-run starting from time at which geometry is changed (time = 0.5), solution behaves well. So it seems that there is some dynamic response after mesh is changed, which seems to be incorrect...

Thanks again for your time.
Ziemo

[ziemowitzima]

Dear Bruno,

Please see animation at:
http://youtu.be/mVWLc1LkjtQ

where I am solving heat equation only:

fvScalarMatrix omEqn
 ( dimDt*fvm::ddt(om) - 
  fvm::laplacian(dimensionedScalar("1",dimensionSet(0, 2, 0, 0, 0),1), om)
 );  
 omEqn.solve();

Mesh is updated/changed only once during the solution process at time = 0.5.
I am keeping the same BC all the time:
on the upper, left and right boundary there is Dirichlet BC: om = 0,
on the bottom boundary there is Drichlet as well, but as function of space:
om(at bottom boundary) = f(s).

You can see its shape/values e.g. on the frame for time = 0.
Its max absolute value is -13.86...

In my opinion solution is not correct after mesh is updated, because I would expect that during the solution "om" field should get weaker, and max value should be the one which is on the boundary.

Anyway my goal is to start calculations "from the beginning" after mesh is changed. Namely I do not want to have any dynamics response because of changed domain during time integration.

Best regards
Ziemo

[ziemowitzima]

Hi Bruno,
I am still "fighting" with this problem.
But in fact, it seems, that I need something simple.
I need only to change my mesh (as I am sucesfuly doing) and I want to start my calculations "from the beginning" with the new mesh.

The problem I have is some "memory" (mesh change response) which is passed to my solution after mesh is updated. I do not need this "memory".

Is there any way to update mesh, but cancel this "mesh change response" ??

Thanks and best
ZMM

[wyldckat]

Hi ZMM,

Sorry for the late reply, but this is a bit out of my skill-level

Nonetheless, I only very vaguely know of a few things that might help you:

• Some dynamic mesh solvers use "flux correction", if I'm not mistaken. I haven't managed to figure out what this exactly does, but I suspect that you either need this or need to disable this.
• mapFields might help, in the sense that it will ignore that the mesh was moved. The problem is: what will it assign to the cells that were created.
• Have you checked what are the pressure values on the cells that are moved the most?
• Do you have the "forces" function object turned on?

And without a test case, it's even harder to help you any further.

Best regards,
Bruno

[ziemowitzima]

Hi,
Thanks for your reply and help !

In my dynamicMeshDict I have:

twoDMotion      yes;

dynamicFvMesh dynamicMotionSolverFvMesh;

motionSolverLibs ("libfvMotionSolvers.so");

solver displacementLaplacian;

diffusivity  uniform;

I failed to find any info about "flux correction" option, but it sound like something useful in my case.

I am not sure what it is:

Do you have the "forces" function object turned on?

I do not have pressure in my solver.
I am not sure if my "case" would be of any help. I created my own solver to calculate some equations/problem arising from asymptotic simplification of high Re cuette flow.

In short there are three advection-diffusion equations which need to be solved on the domain which is changing as a result of a solution of one of these equations.

For now I switched to my own code in curve-linear coordinates. It works fine, but I would like to use OpenFOAM as well.

Best

[wyldckat]

Hi ZMM,

The forces function object needs at least the "U" and "p" fields to work. If you do not use the "p" field, then you'll need to create a variant of the forces function object, since you'll need to calculate the force some other way than relying on "U" and "p".

As for flux correction... see the tutorial file "./multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/dynamicMeshDict", entry "correctFluxes".

Best regards,
Bruno