CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   Kelvin-Helmholtz Instability Solver (

wschosta June 7, 2011 16:47

Kelvin-Helmholtz Instability Solver
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!


lazersos June 8, 2011 09:37

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 June 8, 2011 13:38

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 June 8, 2011 13:52

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 June 8, 2011 14:17

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 June 8, 2011 14:41

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

wschosta June 8, 2011 14:57

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 June 8, 2011 15:30

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 June 8, 2011 15:37

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 February 14, 2012 06:58

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.



wschosta February 14, 2012 10:36

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.



All times are GMT -4. The time now is 14:59.