CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Struggling with standingg waves

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

Reply
 
LinkBack Thread Tools Display Modes
Old   October 23, 2014, 13:22
Default Struggling with standingg waves
  #1
New Member
 
Remi Carmi
Join Date: Jul 2014
Posts: 13
Rep Power: 2
rcarmi is on a distinguished road
Hi all,

I am a beginner in OpenFoam. I am trying to simulate waves and as a convergence and stability test I try to run the simple case of a standing wave in a 1x1 tank.

I tried to increase the Mesh accuracy and I measure the convergence two ways :
with the Kinematic energy : I measure it at each time step by modifying the solver .C file (ihFoam in my case but it should give the same result with interFoam) by adding the lines :
Code:
OFstream energyTotIntegral("energyTotIntegral.dat");
  while (runTime.run())
    {
 scalar QQ = fvc::domainIntegrate(rho*mag(U)*mag(U)*.5*alpha1).value();
	energyCinIntegral << runTime.timeName() << " "<< QQ << endl;
and a L2 convergence test on the velocity profile :
Code:
double ux (double a, double om, double k, double h, double time, double x, double z)
{
  double res =0;
  res = a*om*Foam::cosh(k*(h+z))/Foam::sinh(k*h);
  res = res*Foam::sin(k*x)*Foam::sin(om*time);
 // Info<< x << " " << z << " "<< time  <<" "  << res << endl;
  return res;
}
to calculate the standing waves Airy (Stokes I) solution for small perturbations and I measure the error like that :
Code:
// Calculate the error in velocity 
	//First Construct the volscalar field of Ux-uAnalytic
        volScalarField Ux = U.component(1);
        volScalarField Uz = U.component(2);
        volScalarField xcoord = mesh.C().component(1); 
        volScalarField zcoord = mesh.C().component(2); 
        // Info<< "Here" << endl;
        forAll (mesh.C(), celli)
        {
               // Info<< "Here" << endl;
		Ux[celli]=alpha1[celli]*(Ux[celli]-ux(a,om,k,h,runTime.value(),xcoord[celli],zcoord[celli]));
		Ux[celli]=Ux[celli]*Ux[celli];
		Uz[celli]=alpha1[celli]*(Uz[celli]-uz(a,om,k,h,runTime.value(),xcoord[celli],zcoord[celli]));
		Uz[celli]=Uz[celli]*Uz[celli];
        }
        scalar er1 = fvc::domainIntegrate(Ux).value();
        errorUx << runTime.timeName() << " "<< er1 << endl;
        scalar er2 = fvc::domainIntegrate(Uz).value();
        errorUz << runTime.timeName() << " "<< er2 << endl;

My mesh is the following :
Code:
convertToMeters 1;

vertices
(
    (0.500	-0.5   	-1   )
    (0.505      -0.5    -1   )
    (0.505	 0.5	-1   )
    (0.500 	 0.5	-1   )
    (0.500	-0.5   	 0.2 )
    (0.505      -0.5     0.2 )
    (0.505	 0.5	 0.2 )
    (0.500 	 0.5	 0.2 )

);

blocks
(
 hex (0 1 2 3 4 5 6 7) (1 100 120) simpleGrading (1 1 1)  
);
The BC are :
Code:
patches
(
    patch atmosphere
        (
            (5 6 7 4)
        )
    empty frontAndBack
    (
	//Back
         (0 4 7 3)
        //Front
         (1 2 6 5)
    )
    wall defaultFaces
    (
    )
);
and I set up
for alpha1 :
Code:
dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    atmosphere
    {
        type            inletOutlet;
        inletValue      uniform 0;
        value           uniform 0;
    }
    defaultFaces
    {
        type            zeroGradient;
    }
        frontAndBack
    {
        type            empty;
    }

}
for p_rgh
Code:
dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
        atmosphere 
    {
         //type            fixedValue;
         //value uniform 0;
	type            totalPressure;
        U               U;
        phi             phi;
        rho             rho;
        psi             none;
        gamma           1;
        p0              uniform 0;
        value           uniform 0;
    }
     defaultFaces
    {
        type            slip;
    }
       frontAndBack
    {
        type            empty;
    }
and for U
Code:
dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
        atmosphere
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }
     defaultFaces
    {
        type            slip;
    }
        frontAndBack
    {
        type            empty;
    }

}
Finally the schemes I use are the following :
Code:
ddtSchemes
{
    default         CrankNicolson/Euler or backward;
}

gradSchemes
{
     default     cellMDLimited Gauss linear 1;
}

divSchemes
{
    div((rho*phi|interpolate(porosity)),U) Gauss limitedLinearV 1;
    div(phi,alpha)  Gauss vanLeer01;
    div(phirb,alpha) Gauss interfaceCompression;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha1;
}
Ok now Initial Conditions : I set up the field using the command :
Code:
echo Setting the fields...
funkySetFields -time 0 -field alpha1 -expression '0.5-tanh((pos().z-(0.+0.02*cos(2*pi*(pos().y))))/(3/4*1/N))*0.5' -time 0 > funkySetField.log
the tanh should spread the interface on 3 cells and avoid "loss" of mass at the beginning of the simulation.
All the rest of the fields are set to zero.

So I have :
h=1m
A=0.02 m << h
L=1
Dispersion relation gives om = 7.85096 rad/s so f = 1.24952 Hz
I run the simulations for 7 seconds

and I tried the cases
dt = {0.01,0.005, 0.001}
and L/Dy = L/Dz = {50,100,200,500}=N
I use the scheme Euler first and I get rather good results. It is stable for all the cases but dt=0.01 and N=500.

Unfortunately the case backward does not work...
I played with the case CrankNicolson and it is unstable if I put a too big Psi.


Can somebody help me out? Where am I wrong? Wrong BCs? Wrong Schemes?

If you are interested I can send you the case by email for let's say N=100 and dt =0.005 so it is faster for you to try the case.

Thanks

Remi
rcarmi 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
WAVES GENERATION IN A PORT OR BRAKE WAVES. important tecnicaf1 STAR-CCM+ 1 May 27, 2011 07:32
CFD Animations of breaking waves Dommermuth Main CFD Forum 0 June 17, 2009 11:47
Papers request-Generation of water waves by source Mehdi BEN HAJ Main CFD Forum 0 June 11, 2007 12:55
Following waves generation Nico Main CFD Forum 0 October 4, 2006 08:28
Terrible Mistake In Fluid Dynamics History Abhi Main CFD Forum 12 July 8, 2002 09:11


All times are GMT -4. The time now is 20:18.