[blue8803]

I'm trying to solve Exner equation (mass conservation equation of sediment phase) using FAM and mesh move is implemented using dynamic mesh motion. My test case is U-shape open channel flow composed of bottom wall (patch name: bed), side wall (patch name:sidewall), inlet, outlet. First, only flow is simulated until flow is fully developed (it is not problem). After that, turn on my implemented code, morphological changes start. At the beginning, it seems good, but some time later problem occur at the vicinity of the sidewall near the bed (please refer to attached picture). It seems that the sidewall is correctly moved in the mesh motion stage. What causes the problem..?
Thank you for any questions or answers!

Here is boundary condition for motionU.

dimensions [0 1 -1 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{
inlet
{
type fixedValue;
value uniform (0 0 0);
}
outlet
{
type fixedValue;
value uniform (0 0 0);
}
sidewall
{
type slip;
value uniform (0 0 0);
} /// I think it is right..
bed
{
type fixedValue;
value uniform (0 0 0);
}
atmosphere
{
type fixedValue;
value uniform (0 0 0);
}
}




I'm now suspicious of relTol....



solvers
{
motionU
{
solver PCG;
preconditioner DIC;
tolerance 1e-09;
relTol 0.001;
}
}






In the dictionary file dynamicMeshDict, what means the twoDMotion??


dynamicFvMeshLibs ("libdynamicTopoFvMesh.so");

//- Select the type of dynamicFvMesh
dynamicFvMesh dynamicMotionSolverFvMesh; //dynamicTopoFvMesh;

twoDMotion yes;

//- Select the type of motionSolver
//solver laplace;
/*
RBFMotionSolver /dynamicMeshDict/movingPatches
displacementComponentLaplacian
displacementInterpolation /dynamicMeshDict/displacementInterpolationCoeffs
displacementLaplacian /0/pointDisplacement
displacementSBRStress /0/pointDisplacement
laplace
mesquiteMotionSolver /dynamicMeshDict/mesquiteMotionSolver
pseudoSolid /dynamicMeshDict/pseudoSolid
refVelocityLaplacian /0/pointMotionU
velocityComponentLaplacian
velocityLaplacian /0/pointMotionU
*/
solver laplace;
diffusivity quadratic inverseDistance;

distancePatches 1 (bed);

frozenDiffusion off;

//





I reseachred some source code but,, i didn't reach understanding.


src/dynamicMesh/dynamicMesh/motionsolver/motionSolver.H
// Member Functions
//- Return true for 2-D motion
bool twoDMotion() const
{
return twoDPointCorrector_.required();
}


src/meshTools/twoDPointCorrector/twoDPointCorrector.C
/*Description
Class applies a two-dimensional correction to mesh motion point field.
The correction guarantees that the mesh does not get twisted during motion
and thus introduce a third dimension into a 2-D problem.*/


// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
:
mesh_(mesh),
required_(mesh_.nGeometricD() == 2), ///??????????????
planeNormalPtr_(nullptr),
normalEdgeIndicesPtr_(nullptr)
{}

[Santiago]

which solver is this?

[blue8803]

I used the solver I made based on the interDyMFoam or pimpleDyMFoam.

[Santiago]

Originally Posted by blue8803

src/meshTools/twoDPointCorrector/twoDPointCorrector.C
/*Description
Class applies a two-dimensional correction to mesh motion point field.
The correction guarantees that the mesh does not get twisted during motion
and thus introduce a third dimension into a 2-D problem.*/


// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
:
mesh_(mesh),
required_(mesh_.nGeometricD() == 2), ///??????????????
planeNormalPtr_(nullptr),
normalEdgeIndicesPtr_(nullptr)
{}

I think the description of the class is rather clear: the grid nodes will not twist for a 2D problem.

The statement:

 mesh_.nGeometricD() == 2

Is asking whether the GRID is 2D. Obviously the 2D correction algorithm doesnt make much sense in 3D grids, specially in bed erosion problems. I assume you are following Xiafeng Liu's work then, unless you are running a 2D case (where you have wedges/empty patches), setting

twoDMotion yes;

is not correct.

The motion solver seems correct to me, not very convinced of you BCs for motionU but, in order for me to give you an informed opinion, I need to see the implementation of the finiteArea model + the faMesh...

[blue8803]

Originally Posted by Santiago

I think the description of the class is rather clear: the grid nodes will not twist for a 2D problem.

The statement:

 mesh_.nGeometricD() == 2

Is asking whether the GRID is 2D. Obviously the 2D correction algorithm doesnt make much sense in 3D grids, specially in bed erosion problems. I assume you are following Xiafeng Liu's work then, unless you are running a 2D case (where you have wedges/empty patches), setting

twoDMotion yes;

is not correct.

The motion solver seems correct to me, not very convinced of you BCs for motionU but, in order for me to give you an informed opinion, I need to see the implementation of the finiteArea model + the faMesh...

Many thanks for your reply.

This is the piece of my code.

(bedload calculation)
...
(solving Exner equation)
...


(solving mesh motion)

tetPointVectorField& motionU = 
const_cast<tetPointVectorField&>
(
  mesh.lookupObject<tetPointVectorField>
  (
    "motionU"
  )
);

fixedValueTetPolyPatchVectorField& motionUaPatch =
refCast<fixedValueTetPolyPatchVectorField>
(
  motionU.boundaryField()[aPatchID]
);

tetPolyPatchInterpolation tppiAPatch
(
  refCast<const faceTetPolyPatch>
  (
    motionUaPatch.patch()
  )
);

primitivePatchInterpolation patchInterpolator (mesh.boundaryMesh()[aPatchID]);
scalarField faceVal = patchInterpolator.faceToPointInterpolate(deltaH);

vectorField dH (Psize);
forAll(dH, i)
{
  dH[index].x() = 0;
  dH[index].y() = 0;
  dH[index].z() = faceVal[i];
}

motionUaPatch == tppiAPatch.pointToPointInterpolate (dH/mesh.time().deltaT().value());
//motionUaPatch == tppiAPatch.pointToPointInterpolate (dH);

meshChanged1 = mesh.update();
reduce(meshChanged1, orOp<bool>());

(sand-sliding)
...
(interFoam)
...

Also, I have tested some factors.
(i) relTol in dynamicMeshDict
0, 1e-3
(ii) diffusivity in dynamicMeshDict
inverseDistance, inverseFaceDistance, inversePointDistance with linear manipulator
inverseDistance, inverseFaceDistance, inversePointDistance with qudratic manipulator
(iii) twoDMotion in dynamicMeshDict
yes, no
In fact, there is no significant impact on final results considering above different factors..

Now I made a case rather simple. only one face are upward in constant motionU velocity (in code deltaH(somewhere on face) = 0.2), and found that at the sidewall upper layer of the cell moving upward converged each others (see picture).

[blue8803]

I really don't know why they are converged
--------------------------------------

[blue8803]

With turning off the sand-sliding features, mesh motion is works fine. However, I think it is not the problem of the sand-sliding model I made, but the iterative positive/negative values of motionU are the major reasons.
Because there is no steep slope that activate sand-sliding model in the first simulation (U-shape channel flow). So, i dig into this perspective

[blue8803]

Many thanks.