[wschosta]

Hello Foamers! I've been working on a solver that joins the vof navier-stokes equations and the magnetic field equation for a few weeks now and I'm trying to model the Kelvin-Helmoltz instability test. I've run into a few questions that I was hoping someone might be able to point me in the right direction on:

How do you set different initial velocities for different parts of the domain? I understand how to "ramp up" velocities over time but not how to have two different velocities to begin with.

What's a good, Open Source / Free Mesh builder? I've been trying to edit meshes by myself in OpenFOAM and it's getting to the point where I just can't take it anymore. A good meshbuilder would be very helpful.

In OpenFOAM how do you perturb the boundary condition? As I mentioned before, I get how to adjust initial velocities and even "ramp up" velocities over time but adjusting them over the problem domain... not so much.

With periodic boundary conditions, are there any special considerations that you have to take into account within the solver itself or does it not affect the solver? I don't think it would have an affect on the solver itself (other than increasing runtime due to more complicated calculations) but the first few times I tried to run my solver I got some funky errors that I managed to resolve somehow, not really sure.

Sorry if these are basic/easy questions, I'm still getting used to OpenFOAM.

Thanks in advance for your help!

Walter

[lazersos]

This is not a trivial problem in practice, I've been trying to do similar things (not KH yet in OpenFOAM).

First let's start with the grid. I've been using blockMesh to let OF do my work for me. If you want a cartesian box then it's very simple. I've managed to mesh rather complex toroidal domains without issue. So I'd suggest starting with blockMesh as a preprocessor utility and deprecate any notions of doing your own meshes for right now.

Second, let's touch on initial conditions. You'll have to code up your own solver at this point but that's not bad (and it sounds like you may have already started this). Let me give you a snippet of code

forAll(mech.C, celli)
{
     vector c = mech.C()[celli];
     scalar x = c.x();    // X-Position of cell center
     scalar y = c.y();    // Y-Position
     scalar z = c.z();    // Z-Position

     // How to set something like density
     rho[celli] = 0.5*Foam::tanh(z);

    // How to set a vector field like (U)
    U[celli] = vector(0,0.5*Foam::tanh(z),0);
}

Here I've assumed you have a volScalarField rho and a volVectorField U. I hope this helps a bit. I'd add this type of code after you include createFields.H so all your variables exist but before you enter your while(runTime.loop()).

[wschosta]

Thanks for your quick response!

So I've been using blockMesh to actually create the mesh but I've been running into some problems with stating the parameters of the mesh, particularly the initial boundary conditions of the fluid of the (there ends up being something like 60,000 points in my alpha1 file and it seems the only way to adjust the initial boundary conditions is to go in by hand and change them around).

I do have a solver that seems to work fairly well, I'll have to add some functionality to it but it's a work in process so that's not a huge deal.

I'm afraid that the purpose of the code snippet you've provided isn't exactly clear to me. I see that it's setting the density and velocity but could you use something like that to set velocity below a certain domain, for example:

forAll(mech.C, celli)
{
vector c = mech.C()[celli];
scalar x = c.x(); // X-Position of cell center
scalar y = c.y(); // Y-Position
scalar z = c.z(); // Z-Position

// set density
if (y < 0.5)
rho[celli] = 1
if (y >= 0.5)
rho[celli] = 2

// set a vector field like (U)
if (y < 0.5)
U[celli] = vector(-5,0,0);
if (y >= 0.6)
U[celli] = vector(5,0,0);
}

is the "mech.C" something that's built in?

Thanks again, I really appreciate all your help.

[lazersos]

Yes from what I undersand x,y and z are the centers of the cells. Your code snippet should work. You may want to wrap your 1 and 2 in a call to scalar (in case you run into type conversion but I'd think what you have would work).

As for your boundary conditions, you'll want to search around a bit. If you're doing a K-H simulation in a cartesian box then perhaps you want to just set the BC to be cyclic in the direction of the initial velocity field.

Now to explain what my code snippet is doing. First mesh.C I believe are the centers of all the cells. So vector c is just what you'd think a vector variable with the location of a given cell is located (indexed by celli). The c.x references the x position, etc. So in essence you're looping over all the cells setting the values of a variable at each cell location.

Now here is a subtly that I'm not sure of. How do other variables update? For example if I declare phi to be a face flux field (rho*U or just U if incompressible) do I need to add code which updates phi after every iteration? Or is this handled by openFOAM implicitly when one variable is updated?

[wschosta]

As far as I can tell, you're correct about mesh.C being the cell center. I'm having some trouble in terms of implementation, it's returning this error:

error: ‘mesh.Foam::fvMesh::C’ does not have class type

and it's line number is the line that the 'forAll(mesh.C, celli)' is on.

have you tried this code and experienced similar errors?

What I think is the best course of action is to put the correction of values immediately after the declaration of U within createFields.H, this should allow the correct values of U to be used in the creation of other variables and prevent contradictions within the code. As far as I can tell variables that aren't explicitly changed later on in the code, either by redeclaration or by solving, stay constant throughout. I'm not sure if that's right though.

[lazersos]

Yeah it compiles for me. Are you including createMesh.H?

[wschosta]

Yeah, I figured it out - it has to be:

forAll(mesh.C(), celli)

I'm running a test case now but I think that will work. I'm going to use that same method to set the phase boundary, I think that will work well.

Any thoughts on making a wave between the phases? I'm thinking there might be a way to set a 'ramping' change in velocity for areas where the phase is transitioning (the solver I built is a combination of mhd and interFoam so it views boundary conditions).

[lazersos]

For KH simulations it's pretty commonplace to assume some shear width for the velocity (boundary layer). I've used a shifted hyperbolic tangent in other codes. Then to apply some small velocity perturbation perpendicular to the shear layer. What I've done in the past is to use a sin wave along the layer and clamp it with a hyperbolic cosine to the -2 power. This confines the perturbation to the layer width.


As for BC I'd use cyclic in the flow direction and maybe fixed at the walls. It's really a question of just how 3D your flow field is.

[wschosta]

Yeah, what I think i'm going to use is a sinusodial function to put a perturbed boundary in the initial condition and then I'll set the velocity based on the phase condition.

I think right now, in terms of boundary conditions, i have cyclic on the left and right (the flow is in the x direction), a wall below, and atmosphere above. With cyclic boundaries do you need to define and inlet/outlet patch in addition to the cyclic patch?

[hismother]

Hi Walter,

I am trying to simulate the KH-instability using DNS solver in Openfoam. The results become unstable after some timesteps.
I just wanna ask which solver in Openfoam that you used and did you get good results for it.

Thanks in advance.

Regards,

Terrence

[wschosta]

Hi Terrence,

I actually used the mhdfoam and interfoam solvers and built a solver to combine their functionality. I was getting good results but the funding for my position (I'm an undergraduate) ran out so I wasn't able to finish my work.

Regards,

Walter