|
[Sponsors] |
October 23, 2014, 13:22 |
Struggling with standingg waves
|
#1 |
New Member
Remi Carmi
Join Date: Jul 2014
Posts: 15
Rep Power: 11 |
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; 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; } 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) ); Code:
patches ( patch atmosphere ( (5 6 7 4) ) empty frontAndBack ( //Back (0 4 7 3) //Front (1 2 6 5) ) wall defaultFaces ( ) ); 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; } } 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; } 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; } } 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; } 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 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 |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
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 |