[Mirza8]

Hello everyone,

Short description of my problem: I am trying to develop a multi-region solver, where Region 2 is inside the Region 1 (similar to the attached picture), and Region 2 shrinks in size during the simulation. So I need to have a dynamic mesh in my solver, but I do not know how and from where to start. Can anyone shed some light please?

More details about problem: Consider the Region 2 as a large solid particle and Region 1 as the fluid domain surrounding the solid particle. Due to changes in the properties of the solid region during simulation, it shrinks and its radius decreases. The shape of the particle is usually circular or rectangular, and I hope to find a solution that works for both shapes.

I am not familiar with dynamic mesh in OpenFOAM, and I could not find my answer over the web in the past few days. I have found some tutorials which were mainly for very old versions of OF (I am using OF6), and they were mainly describing how to move the mesh as a rigid body, or how to move one boundary in a specific direction! Also I tried to read and understand the code of coldEngineFoam, which again moves the cylinder patch in the z direction! So again it was not very helpful for my problem.

Ideally, I would like to move (deform/resize) the interface patch between the two regions. But if it is not easily possible, changing the definition of the interface between two regions (without modifying the mesh) might also work with a fine enough mesh. What I am trying to say is to find a way to detach some elements from Region 2 and attach them to Region 1 and define the new patch in its new place and continue the simulation! But I have no idea how to start working on the second approach as well, so I would highly appreciate any comments or hints on how I should start working on this problem!


Best regards,
Morteza

[Mirza8]

OK here is what I have been trying in the past few days. As I said before, I want to model two regions with two different meshes. One region (solid region) shrinks over time, and therefore the other region (fluid region) which is around the first region, expands over time. I read the layeredEngineMesh code, and I found that simply I have to get the mesh points, modify them, and use the fvMesh::movePoints function to move the points. So, I have to define a simple function like the following:

void Foam::myMesh::move(scalar dTheta)
{
        // dTheta is the volume change for the solid region
        // get the mesh points
	pointField newPoints(points());

	forAll(newPoints, pointi)
	{
		// Here I calculate the new coordinate of every point
                ....
        }
        // Move the mesh points to the new coordinates
        Info<<"Text 1"<<endl;
        movePoints(newPoints);
        Info<<"Text 2"<<endl;
}

Based on my case setup, I managed to calculate the points new coordinates. For a very simple case, I checked the newPoints and they are correct for both regions. However, I still have problems with the movePoints(newPoints) function (which is derived from the fvMesh class)! The code of this function is difficult for me to understand I have the following questions about it, which I hope some of you might be able to help:

1- What exactly this function does? Does it only move the points of the mesh, or does it also fix the cell volumes and everything related to mesh? What about field variables such as density or pressure of each cell?
2- Does it require special treatment for parallel cases?
3- Is there any alternative to using this function?

I have tried the same function in 3 different solver/case setups, and for each case I get different results:
a- If I don't solve any equations in the solid region, trying to shrink the solid region over time, the function works fine both in serial and parallel. I can see solid region shrinking and fluid region expanding, no problem at all!
b- If I run in serial, and try to solve governing equations in the solid region, even a slight change in the mesh of the solid region causes the solution to diverge!!! I do not understand how for instance a tiny displacement of 1e-6 m (which cannot even be seen compared to the my cell sizes which are in order of 1e-3 m) can cause a previously stable solution to diverge?! Am I missing something here?
c- If I run the (case b) in parallel, the solver stucks in movePoints(newPoints) function, without any errors! i.e. I can see "Text 1" in the output but not "Text 2", even after several hours! Never seen such thing before, and I have no idea what to do!!

I have been struggling with this for few weeks now, and I highly appreciate any advice on this matter!

[einstein_zee]

Originally Posted by Mirza8

OK here is what I have been trying in the past few days. As I said before, I want to model two regions with two different meshes. One region (solid region) shrinks over time, and therefore the other region (fluid region) which is around the first region, expands over time. I read the layeredEngineMesh code, and I found that simply I have to get the mesh points, modify them, and use the fvMesh::movePoints function to move the points. So, I have to define a simple function like the following:

void Foam::myMesh::move(scalar dTheta)
{
        // dTheta is the volume change for the solid region
        // get the mesh points
	pointField newPoints(points());

	forAll(newPoints, pointi)
	{
		// Here I calculate the new coordinate of every point
                ....
        }
        // Move the mesh points to the new coordinates
        Info<<"Text 1"<<endl;
        movePoints(newPoints);
        Info<<"Text 2"<<endl;
}

Based on my case setup, I managed to calculate the points new coordinates. For a very simple case, I checked the newPoints and they are correct for both regions. However, I still have problems with the movePoints(newPoints) function (which is derived from the fvMesh class)! The code of this function is difficult for me to understand I have the following questions about it, which I hope some of you might be able to help:

1- What exactly this function does? Does it only move the points of the mesh, or does it also fix the cell volumes and everything related to mesh? What about field variables such as density or pressure of each cell?
2- Does it require special treatment for parallel cases?
3- Is there any alternative to using this function?

I have tried the same function in 3 different solver/case setups, and for each case I get different results:
a- If I don't solve any equations in the solid region, trying to shrink the solid region over time, the function works fine both in serial and parallel. I can see solid region shrinking and fluid region expanding, no problem at all!
b- If I run in serial, and try to solve governing equations in the solid region, even a slight change in the mesh of the solid region causes the solution to diverge!!! I do not understand how for instance a tiny displacement of 1e-6 m (which cannot even be seen compared to the my cell sizes which are in order of 1e-3 m) can cause a previously stable solution to diverge?! Am I missing something here?
c- If I run the (case b) in parallel, the solver stucks in movePoints(newPoints) function, without any errors! i.e. I can see "Text 1" in the output but not "Text 2", even after several hours! Never seen such thing before, and I have no idea what to do!!

I have been struggling with this for few weeks now, and I highly appreciate any advice on this matter!

Hii Morteza,

I just recently got involved with dynamic meshes in OF and I thought I might add some stuff here. I cannot answer all of your questions but I'll post here things I noticed so far. As I had a look into this movePoints I figured out it goes deep inside OF. fvMesh::movePoints seems to be the tip of the iceberg. polyMesh::movePoints() && primitiveMesh::movePoints are the other layers. You may find them at fvMesh.C:716, polyMesh.C:1291 & primitiveMesh.C:309 respectively. On the lowest level (primitiveMesh::movePoints) the swept volumes are calculated based on old and new faces (this will fill sweptVols). Then there is some forced recalculation of the geometric data(boundries, pointZones, faceZones , ...). The final product at this stage will be a tmp scalarField which holds the swept volume because of the motion. Out of this OF constructs the mesh flux(phiMesh). Regarding your last point also, I should say in btw there are some mpi communications (when you trace the funcs deeply) and the reason you see that it hangs in the air is because of a hanging mpi communication which could be caused by the things you added. I hope others also get involved in this for more clarification.

[Mirza8]

Hi Hosein,

Thank you very much for your response. So here is what I understand from your explanation:
1- The movePoints function does all the mesh-related calculations such as volume change, face displacements, or phiMesh.
2- But movePoints, does not change any other field variable, such as density field, which will be changed due to volume change. Therefore, I have to use phiMesh as source term in different governing equations to fix them.
Did I understand correctly? At the moment, I have not used phiMesh anywhere in my code, could that be the reason for divergence?

And regarding the last part of your response:

Originally Posted by einstein_zee

...Regarding your last point also, I should say in btw there are some mpi communications (when you trace the funcs deeply) and the reason you see that it hangs in the air is because of a hanging mpi communication which could be caused by the things you added. I hope others also get involved in this for more clarification.

I cannot understand clearly what part of the things I added could cause mpi problems? The only code that I added is the modification of points coordinates:

    pointField newPoints(points());    
    forAll(newPoints, pointi)
    {
        point& p = newPoints[pointi];

	// Move x coordinates
	if (p.x() >= x1 && p.x() <= x2)    		// Particle
        {
            p.x() -= deltaX / 2.0 * (p.x() - xc) / (x2 - xc);
        }
		else if (p.x() > x0 && p.x() < x1)		// Left side of particle
        {
            p.x() += deltaX / 2.0 * (p.x() - x0) / (x1 - x0);
        }        
        else if (p.x() > x2 && p.x() < x3)		// Right side of particle
        {
            p.x() -= deltaX / 2.0 * (x3 - p.x()) / (x3 - x2);
        }
        
        // Move y coordinates
	if (p.y() > y1 && p.y() <= y2)		// Particle
        {
            p.y() -= deltaY * (p.y() - yc) / (y2 - yc);
        }
        else if (p.y() > y2 && p.y() < y3)		// Top side of particle
        {
            p.y() -= deltaY * (y3 - p.y()) / (y3 - y2);
        }
        
        // Move z coordinates
        if (p.z() >= z1 && p.z() <= z2)			// Particle
        {
            p.z() -= deltaZ / 2.0 * (p.z() - zc) / (z2 - zc);
        }
        else if (p.z() > z0 && p.z() < z1)		// Back side of particle
        {
            p.z() += deltaZ / 2.0 * (p.z() - z0) / (z1 - z0);
        }
        else if (p.z() > z2 && p.z() < z3)		// Front side of particle
        {
            p.z() -= deltaZ / 2.0 * (z3 - p.z()) / (z3 - z2);
        }
    }
    
    // Update coordinates of the particle, y1 = yc is stationary
    x1 += deltaX / 2.0;
    x2 -= deltaX / 2.0;
    
    y2 -= deltaY;
    
    z1 += deltaZ / 2.0;
    z2 -= deltaZ / 2.0;
	
    Info<<"Text 1"<<endl;
    Foam::fvMesh::movePoints(newPoints);
    Info<<"Text 2"<<endl;

[einstein_zee]

Hi Morteza,

Point 1 seems reasonable. Point 2 also makes sense but you may notice that for the solvers that OF states explicitly that they support movingMesh such a thing should already be considered so you don't need to take care of it. Well, I cannot judge at this point that the divergence issue is because of it.
And finally for the last point, things happen . You may try the hanging case again by adding this to your controlDict and see if anything changes.

OptimizationSwitches
{
    commsType        blocking;
}