CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Moving boundary problem based on calculated data

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree15Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   August 22, 2013, 04:52
Default Moving boundary problem based on calculated data
  #1
New Member
 
Aleksander Sandro Grm
Join Date: Feb 2010
Location: Slovenia, Kanal ob Soči
Posts: 9
Rep Power: 7
sandro.grm is on a distinguished road
Send a message via Skype™ to sandro.grm
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.
sandro.grm is offline   Reply With Quote

Old   August 22, 2013, 15:03
Default
  #2
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
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:

Code:
OpenFOAM-1.6-ext/applications/solvers/surfaceTracking/freeSurface
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

Code:
OpenFOAM-1.6-ext/applications/solvers/surfaceTracking/freeSurface/freeSurface.C
on line 647.

Good luck,

Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Old   August 22, 2013, 16:28
Default
  #3
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 8
fredo490 is on a distinguished road
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:
Code:
    // Get the patch ID where impingement occurs: patch name must be "wall"
    label patchWallID = mesh.boundaryMesh().findPatchID("wall");
    const fvPatch& patchWallFaces = mesh.boundary()[patchWallID];

    Info<< "#--------------- Displacement ---------------# " << endl;

//Find the reference to the location of pointDisplacement field
pointVectorField& PointDisplacement = const_cast<pointVectorField&>
(
	mesh.objectRegistry::lookupObject<pointVectorField>
	(
	"pointDisplacement"
	)
);

//Get the vector field of the patch
vectorField &pDisp=refCast<vectorField>(PointDisplacement.boundaryField()[patchWallID]);

//Find the relevant size of the vector and declare a vectorField.
int Psize= pDisp.size();
vectorField dispVals(Psize);

// ########################################################################## //

//- set-up interpolator
primitivePatchInterpolation patchInterpolator (mesh.boundaryMesh()[patchWallID] );

scalarField MimpPatch = Mimp.boundaryField()[patchWallID];

//- perform interpolation
scalarField faceValues = patchInterpolator.faceToPointInterpolate(MimpPatch);

vectorField &PointPointer = refCast<vectorField>(PointDisplacement.boundaryField()[patchWallID]);
vectorField PointNormalVector = mesh.boundaryMesh()[patchWallID].pointNormals();

// loop over points to move the nodes
forAll(dispVals, index)
{
	if(faceValues[index] > 0.001)
	{
	dispVals[index].x() = PointPointer[index].x() 
	 - faceValues[index] * PointNormalVector[index].x() * runTime.deltaT().value() / rho_mat.value();
	dispVals[index].y() = PointPointer[index].y() 
	 - faceValues[index] * PointNormalVector[index].y() * runTime.deltaT().value() / rho_mat.value();

		if(Problem2D.value() == 1)
		{
		dispVals[index].z() = 0.0;
		}
		else
		{
		dispVals[index].z() = PointPointer[index].z() 
		 - faceValues[index] * PointNormalVector[index].z() * runTime.deltaT().value() / rho_mat.value();
		}
	}
	else
	{
	dispVals[index].x() = PointPointer[index].x();
	dispVals[index].y() = PointPointer[index].y();

		if(Problem2D.value() == 1)
		{
		dispVals[index].z() = 0.0;
		}
		else
		{
		dispVals[index].z() = PointPointer[index].z();
		}
	}
}

// ########################################################################## //

//Once the values have been assigned to dispVals, assign them to cellDisplacement boundaryField
PointDisplacement.boundaryField()[patchWallID] == dispVals;

mesh.update();
mesh.moving(false); //set moving = false to allow localEuler ddt - only for negligible displacement compared to the cell size
		
// ########################################################################## //
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.
deepinheart, Tobi, giack and 3 others like this.
fredo490 is offline   Reply With Quote

Old   August 22, 2013, 16:33
Default
  #4
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 8
fredo490 is on a distinguished road
Here is my dynamicMeshDict file that is in Constant

Code:
dynamicFvMesh      dynamicMotionSolverFvMesh;

motionSolverLibs ("libfvMotionSolvers.so");

solver            displacementLaplacian;

displacementLaplacianCoeffs
{
    diffusivity     quadratic inversePointDistance (wall);
}
and an example of "pointDisplacement" file that is in 0:

Code:
FoamFile
{
    version     2.0;
    format      ascii;
    class       pointVectorField;
    location    "0";
    object      pointDisplacement;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 0 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    sides
    {
        type            empty;
    }
    wall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}
fredo490 is offline   Reply With Quote

Old   August 23, 2013, 06:45
Default
  #5
New Member
 
Aleksander Sandro Grm
Join Date: Feb 2010
Location: Slovenia, Kanal ob Soči
Posts: 9
Rep Power: 7
sandro.grm is on a distinguished road
Send a message via Skype™ to sandro.grm
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
sandro.grm is offline   Reply With Quote

Old   August 23, 2013, 11:12
Default dissolving walls
  #6
Member
 
Tony Ladd
Join Date: Aug 2013
Posts: 34
Rep Power: 4
tladd is on a distinguished road
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
tladd is offline   Reply With Quote

Old   August 23, 2013, 16:35
Default
  #7
Member
 
Tony Ladd
Join Date: Aug 2013
Posts: 34
Rep Power: 4
tladd is on a distinguished road
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
tladd is offline   Reply With Quote

Old   August 23, 2013, 16:39
Default
  #8
Member
 
Tony Ladd
Join Date: Aug 2013
Posts: 34
Rep Power: 4
tladd is on a distinguished road
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
tladd is offline   Reply With Quote

Old   August 23, 2013, 16:50
Default
  #9
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 8
fredo490 is on a distinguished road
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.
giack and PicklER like this.
fredo490 is offline   Reply With Quote

Old   August 23, 2013, 16:52
Default
  #10
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 8
fredo490 is on a distinguished road
I don't have my finaly source code on this computer but I think this can solve your problem:

options files:
Code:
EXE_INC = \
    -I../rhoPorousMRFPimpleFoam \
    -I.. \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
    -I$(LIB_SRC)/finiteVolume/cfdTools \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/fvOptions/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/fvOptions/lnInclude

EXE_LIBS = \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -ldynamicMesh \
    -lfluidThermophysicalModels \
    -lspecie \
    -lcompressibleTurbulenceModel \
    -lcompressibleRASModels \
    -lcompressibleLESModels \
    -lfiniteVolume \
    -lmeshTools \
    -lsampling \
    -lfvOptions
and my includes
Code:
#include "meshTools.H"
#include "primitiveMesh.H"
#include "primitivePatch.H"
#include "dynamicFvMesh.H"
#include "primitivePatchInterpolation.H"
#include "pointMesh.H" //Perhaps not needed..?
#include "pointFields.H" //Perhaps not needed..?
#include "pointPatchField.H"
#include "volPointInterpolation.H"
#include "fixedValuePointPatchField.H"

#include "fvCFD.H"
#include "psiThermo.H"
#include "turbulenceModel.H"
#include "fvIOoptionList.H"
#include "fvcSmooth.H"
#include "pimpleControl.H"
#include "bound.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"

    pimpleControl pimple(mesh);

    #include "setInitialrDeltaT.H"
    #include "createFields.H"
    #include "createFieldsIcing.H"

    #include "createFvOptions.H"
    #include "initContinuityErrs.H"
ngj and Fanfei like this.
fredo490 is offline   Reply With Quote

Old   August 23, 2013, 18:27
Default
  #11
Member
 
Tony Ladd
Join Date: Aug 2013
Posts: 34
Rep Power: 4
tladd is on a distinguished road
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.
tladd is offline   Reply With Quote

Old   August 23, 2013, 23:27
Default
  #12
Member
 
Tony Ladd
Join Date: Aug 2013
Posts: 34
Rep Power: 4
tladd is on a distinguished road
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
tladd is offline   Reply With Quote

Old   August 24, 2013, 04:24
Default
  #13
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 8
fredo490 is on a distinguished road
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
fredo490 is offline   Reply With Quote

Old   August 24, 2013, 13:18
Default
  #14
Member
 
Tony Ladd
Join Date: Aug 2013
Posts: 34
Rep Power: 4
tladd is on a distinguished road
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
tladd is offline   Reply With Quote

Old   August 25, 2013, 02:57
Default
  #15
New Member
 
Aleksander Sandro Grm
Join Date: Feb 2010
Location: Slovenia, Kanal ob Soči
Posts: 9
Rep Power: 7
sandro.grm is on a distinguished road
Send a message via Skype™ to sandro.grm
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.
sandro.grm is offline   Reply With Quote

Old   August 25, 2013, 05:55
Default
  #16
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 8
fredo490 is on a distinguished road
It is a C++ trick. It is an overloaded operator that force the dispVals to go in the pointdisplacement

Last edited by fredo490; August 25, 2013 at 07:15.
fredo490 is offline   Reply With Quote

Old   August 25, 2013, 07:21
Default
  #17
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 8
fredo490 is on a distinguished road
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.
fredo490 is offline   Reply With Quote

Old   August 25, 2013, 13:08
Default
  #18
New Member
 
Aleksander Sandro Grm
Join Date: Feb 2010
Location: Slovenia, Kanal ob Soči
Posts: 9
Rep Power: 7
sandro.grm is on a distinguished road
Send a message via Skype™ to sandro.grm
Dear Frederic,

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

Besta, ASG
sandro.grm is offline   Reply With Quote

Old   October 9, 2013, 03:27
Default
  #19
New Member
 
Aleksander Sandro Grm
Join Date: Feb 2010
Location: Slovenia, Kanal ob Soči
Posts: 9
Rep Power: 7
sandro.grm is on a distinguished road
Send a message via Skype™ to sandro.grm
And complete example is here. Code for myMoveMesh and test. Test is based on mesh motion due to dissolution at the interface.

Cheers, Sandro.
Attached Files
File Type: gz moveMeshExample.tar.gz (60.5 KB, 219 views)
sandro.grm is offline   Reply With Quote

Old   October 10, 2013, 05:30
Default
  #20
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 8
fredo490 is on a distinguished road
Nice contribution Have you tried to run it in parallel ?
fredo490 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
An error has occurred in cfx5solve: volo87 CFX 5 June 14, 2013 17:44
domain imbalance for enrgy equation happy CFX 14 September 6, 2012 01:54
Framework for moving mesh based on 2D-computation as a boundary condition Arnoldinho OpenFOAM Running, Solving & CFD 0 May 17, 2011 12:48
RPM in Wind Turbine Pankaj CFX 9 November 23, 2009 05:05
How to update polyPatchbs localPoints liu OpenFOAM Running, Solving & CFD 6 December 30, 2005 18:27


All times are GMT -4. The time now is 08:09.