Moving boundary problem based on calculated data
 

Dear Foamer's,

I'm working on a problem, where bottom of the channel is dissolving due to some physical mechanism ( surface velocity is proportional to element concentration difference/gradient). I would like to implement the mesh move solver (velocity Laplacian solver) where the "bottom" patch movement is based on a calculated data.

Until now, I was able to fix moveMesh utility to work with constant speed on a bottom boundary.

The idea is couple moveMesh utility and scalarTransportLaw solver. Somehow to read all inputs from dynamicMeshdict 0/pointMotionU (here the speed at bottom is 0 and the rest of boundary is as needed - fix/slip) and my concentration data (0/c1) solved previously by scalarTransportLaw. In the next step I would like to prescribe new data for pointMotionU on a "bottom" patch, based on calculated data ( vel = a*abs(c1 - c0) ) and to update mesh.

I was looking on a code of moveMesh utility, but have no idea how to update new values of pointMotionU to "bottom" patch to be able to update mesh.

Any help in this direction will be very appreciated.

Thanks + Best,
Sandro.

Hi Sandro,

I am not aware of which solver you are using, but the following in the version 1.6-ext would give you a lot of inspiration on your further work:

Please note that this solver is not available in any of the releases from ESI. The most noticeable part is the handling of the different mesh motion techniques. Some specify motion in the face centres, others in the mesh vertices and then some in both at the same time. Obviously, setting boundary conditions for these need individual attention. This part begin in

on line 647.

Good luck,

Niels

Just a precision, I think that what you want to do is "deforming a boundary" and now "Moving a boundary" ! It's not the same...

I did something similar for one of my project that I will publish very soon. I didn't use the point velocity but the point displacement (I give the distance vector of the displacement and not the velocity vector). I had some accretion on a wall that needed to be simulate.

I made a code that move the nodes of a wall (patchWallID) from a variable distance in the "normal" direction. The distance of displacement is given by my solver (Mimp = impinging mass flux) and the nodes are moved in the normal direction (= average vector of the faces normal vector).


I'm sure the code can be improved but at least it works well for OF 2.2.

Here is the part of my code that might interest you:
Ps. I set mesh.moving to false so that OF doesn't know that I'm actually moving the mesh. I do this because some tools I use are incompatible with the moving mesh such as the localEuler ddt (but I choose to ignore it since my mesh motion is very small at each time step).

Edit. I also do an interpolation since the value of my displacement (Mimp) is only given at the face of the wall; I use a linear interpolation to get the value at the nodes.

Edit2. I set a variable to define if the simulation is 2D or 3D. Indeed, because I have a lot of iterations, a very small error on the normal Z value can slowly move the wall nodes away from their original plane. I use an external variable to force the Z direction to remain 0.

Here is my dynamicMeshDict file that is in Constant

and an example of "pointDisplacement" file that is in 0:

Dear Niels and Frédéric,

thanks for help. If I would run in some problems I'll make a post. I hope to have nice results to post the solution.

Thanks again.

Best, Sandro.

PS. I use OF 2.2.x

dissolving walls
 

Dear Foamers

I also have a simulation where I want to move the solid boundaries in response to a reactant flux across the surface. I have just started playing with OpenFOAM and I had been looking for some guidance, in vain up until yesterday. Like Sandro, I had thought to play with the input files and the moveMesh/moveDynamicMesh utilities but I concluded (perhaps incorrectly) that this would not be as simple as it first appeared and that it would be better to move the mesh in the solver.

I was investigating the surfaceTracking extensions when I came across this thread (does anyone know why these applications have not been included in the EFI release since this seems to be a missing functionality).

When Frederic says moving the mesh is different from deforming the mesh is this the point that "moving" can include topological change where "deforming" cannot?

Thanks to Frederic for posting his code - it is really helpful to see the links to the OpenFOAM data structures. I need to study it more but I have a few questions:

1) Is the field Mimp the reactant concentration? That would make sense to me but I was not sure.
2) Was there a reason to interpolate the concentration and normals rather than calculate the displacements of the face points and then interpolate?
3) In regard to turning off the mesh motion; is this because OpenFOAM will otherwise move the mesh in response to the fluid velocity field? Or is there some other reason.

I was very pleased to see this thread. This is a long term project for us so I hope we will eventually be able to contribute something.

Tony

Frederic

Can you post the include files needed to make your application run. I have been trying to get my concentration solver to compile with your code in it, but without success. I tracked down primitivePatchInterpolation.H but I am stuck on the PointDisplacement class which does not seem to be listed in the OpenFOAM classes. I was wondering if it was a class you made.

Thanks

Tony

PS - sorry about the duplicate post

Frederic

I was wondering if you would post the include files needed to make your application run. I have been trying to compile my solver with your code in it to move the mesh, but without success. I tracked down primitivePatchInterpolation.H but I am stuck on the PointDisplacement class, which does not seem to be listed. Is it a class you defined?

Thanks

Tony

1) Mimp is a mass flux impinging the wall.
I did a freezing simulation where I had some supercooled droplet impinging a wall. I had a first solver to compute the water droplet trajectories, then I compute the water mass flux impacting the wall (= Mimp) and finally I freeze the water.

I had for example 1kg/(m2.s) of water impacting the wall. Then I freeze the water. By knowing the density of the ice (kg/m3) I get a displacement in m/s which is the ice growth rate. Then I multiply by the time step lenght to get the displacement of the point in m.

2) the normal at the node is given by the function "pointNormals", I don't interpolate it. It is already included in OpenFoam.
I interpolate the displacement from the face to the node because I compute the impingement at the face center. I don't have other choice.

3) I turn off the mesh motion because some schemes are not compatible with the mesh motion. Indeed, you should consider the motion of the boundary in your momentum equation. However, in my simulation, the displacement was so small at each time step that I could ignore it.
In my case, I used the localEuler time scheme that make a test if mesh motion is on. If I leave it on, the scheme will crash...
You should actually leave it on, but in some cases you can make the assumption that the motion is so small that it because negligeable.


Ps. a "moving boundary" is a boundary that doesn't deform. It only translate or rotate... In this case you can use some technics such as the layer addition and so on. In a "deforming boundary", the boundary can actually change shape. Each node of the boundary can be manipulated independantly.

I don't have my finaly source code on this computer but I think this can solve your problem:

options files:
and my includes

Frederic

Thanks a lot. The code compiles and executes now. I am still debugging but I can see I have the correct fluxes at the boundaries. That was massively helpful.

Tony

PS Is there a systematic way of figuring out which include files (and library modules) one needs? I was missing many of them.

Frederic

My channel is dissolving now and producing a reasonable looking shape. There is still quite a bit to do but you got me over a huge hurdle. Again many thanks.

Tony

Unfortunately no... I usually use this website to guess: http://foam.sourceforge.net/docs/cpp...7a38c0d80e0ba7

but often I had to put all the includes from openfoam and then slowly remove the useless one.

I also had to struggle during a while to make this code ;)

That's what I suspected. I guess in time I will learn the major classes and their includes.

I plan to make a dissolution/precipitation solver (i.e. Stokes flow, Steady-State convection diffusion, mesh motion from boundary flux). I have all the bits now thanks to you - I just need to assemble it into a proper program. I will post it when I get done.

Tony

Dear Frederic,

I can't find the meaning of == operator in line

PointDisplacement.boundaryField()[patchWallID] == dispVals;

I saw it someware but at the moment I can't find it.

I have in mid that it does something around interpoaltion/update of the mesh.

Thanks + best,
Sandro

PS. Code is runing great. I'll post a complete example after having clean form.

It is a C++ trick. It is an overloaded operator that force the dispVals to go in the pointdisplacement

To make the things clear:
1) PointDisplacement.boundaryField()[patchWallID] == dispVals;
it gives the new position to the boundary (only the boundary moves and not the inner domain)

2) mesh.update();
This will update the mesh by solving an equation. The boundaries are fixed and the inner domain is the unknown. The solver will move the domain nodes to match the requirements.


Be careful, we often consider that the boundary should not move of more than half of the first cell height. If your motion is of less than 1% of the first cell height, you can start to make the assumption that the mesh motion is negligeable compared to you flow motion. You can then turn off the mesh motion by "mesh.moving(false)" in order to be able to use some advanced schemes that were not coded for dynamic mesh.

Edit, you should also be careful while running the code in parallel. I'm not 100% sure it is ok. I have seen some strange behavior if the deforming boundary belong to two CPU domains. The node at the interface doesn't move smoothly.

Dear Frederic,

thanks for the explanation. I'll test for parallel too and post the results.

Besta, ASG

And complete example is here. Code for myMoveMesh and test. Test is based on mesh motion due to dissolution at the interface.

Cheers, Sandro.

Nice contribution ;) Have you tried to run it in parallel ?

Yes it works. Results are identical as for single run. I added figures and log files. I run my Gentoo Linux in VM with 4 cores, so the scaling is not linear but is significant.

Cheers, Sandro

Thank you Aleksander

Really helped me. I see this thread was last updated 2013, but I will try. 2 questions:

1. Is it possible to move the boundary in the direction of the boundary cells normal vectors? I want to have a cylinder which expands.

2. How can I specify the "expansion"? I want the boundary to expand according to a calculated variable.

So far I figured out that each boundary will move with a fixed constant (specified in pointDisplacement) if the while-loop is placed just after the autoPtr<motionSolver> in the solver and the endTime is increased, only the boundary (specified as patchName) will continue to move. This is fine, since I can allocate the outside of the cylinder as one boundary. Still figuring out how the constants alpha and Tmax work, since their constants from createField.

thank you
Vrede

Hi,
1. It is what my code on the first page does. There is one variable called pointnormal. This takes the normal vector at the node.
But in your case I would suggest to use the radial direction instead of the normal. the normal vector theoretically matches the radial direction but from my experience I would not trust it.
Simply give the location of the center of the cylinder (it can be hard coded) and after calculate the vector between the node and the center. This will be your expansion vector.

2)again look at the first page. I do an interpolation from the face center to the node. The face center is the "calculated variable"

Hi Frédéric

Your method helped, thank you. So I experimented with moveMesh solver.

Below are 3 screenshots, the first is the original mesh, which consists of a block at the bottom and a block on top with a slight skew top boundary. The next screenshot shows the top and bottom boundaries move according to die pointDisplacement file in the 0 directory. Here is a snippet:

As I have said before, the movement defined in the pointDisplacement file will only occur once, therefore all values can be set to zero especially if the speed of the movement will be determined in the solver. The 3rd screenshot show the 4th iteration where only the top boundary move.

As seen, it is only the top "layer" of cells which expand. My question: Is it possible to enlarge all the cells in the top block as the boundary moves?

Thank you
Vrede

My previous problem has become somewhat irrelevant since my mesh movement is very small and my solver takes the change in volume into consideration.

Not sure I need to start a new thread...

My next challenge/question is how to "fix" a mesh when two sides combine (see figures) after movement. Is there a way to remove the cells (or nodes) when they "enter" the mesh?

Kind regards
Vrede

Hey,

Very nice work of you and Frédéric! Thanks a bunch!

I'm just a bit wondering about a few things that I was wondering about. Why did you not use the
mesh.update() function and instead do a timeloop? Is it better?

And what is the use of the buf variable? The nondimensional alpha value?

I know this is an old thread but I hope I can get a closer understanding for this piece of code !

Thanks
/Simon

Hi,
I don't really know who you are talking to and what you are referring to.

You need to have one mesh.update at each time step because it is this command that actually move the mesh. Therefore the mesh.update has to be inside a time loop.

I don't understand what is the buf variable you mentioned. Also I don't know where is your alpha. Usually alpha is a non dimensional number that is called the collection efficiency. This number tells you how much matter impact your boundary compared to the mass flux in the free stream.
It is defined as the mass flux in the free stream (kg.m-2.s-1) divided by the mass flux at your boundary.

Thanks for the answer!

Sorry for the confusion, first post here so didn't really understand how a "quick reply" would work. My questions was mainly for Aleksander regarding the changes he did in his version of the code compared to your code. Thnaks for the clarification though! :)

Hi HECKMANN Frédéric:
I use your method to move mesh on the boundary of my solver. however, there is a error as i test a case. could you give me a hinte, where should i to modify the code. Thanks.

best Regards
Fan fei

Hi, it seems that you used the extended version of openfoam. Unfortunately I am not familiar with this version...

In your solver, did you include all the library that I used?

Edit, see the post #10 of page 1

Hi, I copy all the library you used. And I found some library are different in extend version. so I'm trying to find where is wrong.

Best regards
Fan Fei

pointMotionU
 

My postdoc and I have been working on a moving mesh code for a while now, starting with Frederic Heckmann's post as a template. Let me just make a few general points:

1) pointMotionU is better than pointDisplacement. In long runs with large mesh motion a discrepancy builds up between the positions of the cell centers calculated by averaging the vertexes (how OF determines them) and the original position + pointDisplacement. In the end the code breaks in situations where it runs fine with pointMotionU. I think pointDisplacement is intended for solid motion (small displacements). pMU works in differential mode and there is no inconsistency.

2) faceToPointInterpolate is a method of the primitivePatch class. Thus it does not work properly across coupled patches such as cyclic or processor. Connected with this is that pointNormal is a vector field (with no information about boundaries, just a list of data). Thus you cannot get OF to apply coupled patch type boundary conditions and pointNormal is always wrong across cyclic or processor patches.

3) Combining these ideas, a better strategy is to calculate the MotionU field at the face centers using the faceNormals. Then interpolate to the points to get the pointMotionU field. Since pMU is a pointVectorField (with boundaries) we think you can use OF's methods for coupled patches. Right now we are doing this by hand but I think we will have a proper solution soon.

4) One last point. OF always seems to apply the correct boundary conditions to the points field. So if pointMotionU has a component across a cyclic patch (for example), the points are still OK (ie on the cyclic boundary). But since pMU is wrong the Laplace solver gets the wrong solution and the internal points go haywire.

Hope this helps a few people

Tony

Hi HECKMANN Frédéric:
I try all the ways i know. But the problem didn't solve. the header field as below:
and options are
I don't kown where is wrong. Could you help me.

Best regards
Fan Fei

Hi HECKMANN Frédéric:
I find the reason of this error, which is caused by the incompatiable of creatDynamicFvMesh.H and creatFamesh.H. The creatFamesh is used to create a areField to use FAM. As deleted #include "creatFamesh.H", it will be work. But the warnning is still there. As dynamicMesh I has a question. in our program pointdisplacement is used, that means the point of bounday patch will move, however, in fvSolition the cellDisplacement is used, is that mean the value will interpolated from point to cell(Vol), is that right? I so sorry to trouble you again, but i hope i can get some help. Thanks.

Best regards
Fan Fei

In openfoam and most CFD software, you can only move the nodes. Indeed, the cell is built from the nodes coordinate.

Usually we only move the nodes of a boundary and after we ask a "solver" to find the position of the nodes of the domain. The solver will move the nodes of the domain (not those of the boundaries) according to a law that you specify. The Laplace law is often used because it is robust (it consider the rotation) but you can also use a proportional displacement that will act as a spring that push the nodes with a force that decrease as you go away from the boundary. There are tens of morphing solver and each one has its advantages and drawbacks. You have to choose it depending of the motion you are expecting.

Hi Frédéric,

Thanks for your nice job. I could adapt your code to my application. The only problem is that the mesh is moving as what is defined but only at the first time step.


while (runTime.loop())
{
PointDisplacement.boundaryField()[patchID] == dispVals;
Info<< "Time = " << runTime.timeName() << endl;
//mesh.movePoints(motionPtr->newPoints());
mesh.moving();
mesh.update();


runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Hi,
I think that your solver is ok. The problem might be located in your "0" folder and in the "pointdisplacement" file. From my guess the "pointdisplacement" file is not updated/forwarded to next time step.

Hi to all,
I'm doing simulations moving a wall boundary where the dispacement of the movement is chosen by a variable evaluated during the simulation. The first movement steps work very well but at certain iteration I encounter the same problem shown by Vrede in the post #25. I was thinking to use add/remove layers to solve it.

Does anyone have some ideas to solve it? Vrede, did you solve the problem?

Thanks to all
Regards

Giacomo

Hi Giacomo

No, sorry, I have not solved this problem yet. I do not know if there is any existing utility to remove points when faces overlap. You could try to identify and remove these points and then define the corresponding faces with the new point, but this will be brute force.

Sorry, I could not help more

Dear all,

I find it interesting to have a thread, where these things are being discussed. I have noticed above that it was suggested to use the faceToPointInterpolation method, which is available in the PrimitivePatchInterpolation class.

I have just finalised a paper where the mass conserving properties of this exact interpolation scheme was analysed (http://onlinelibrary.wiley.com/doi/1....4015/abstract). The conclusion was that the scheme does not conserve the volume of sand in the bed, when the bed level change is interpolated from the faces to the points. This lack of mass conservation is not uniquely related to sediment transport and morphodynamics.

There are several issues related to error, but to mention a few:

• There is not made a projection of the distance weights onto the horizontal plane. This means that the interpolation weights change over time, as the boundary deforms. This even though the displacements are perpendicular to the horizontal plane.
• In 2D simulations, the width of the computational domain in the empty direction has an effect on the mass conservation error. This is solely related to the fact that the weights are from the faces to the points. For instance, on non-equidistant grids and very wide meshes in the empty direction, the interpolation becomes a simple average (identical weights for all points). This is obviously not correct, since the mesh is defined as non-equidistant.

An alternative interpolation method is proposed in the same paper.

I hope that these finding are valuable for some of you.


Kind regards,


Niels