Need help in plotting graph of velocity (ux) vs time

General Discussion
1 
	///////////////////////////////////////////////////
	///  UNSTEADY 2D NAVIER--STOKES                ////
	///////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////
/////////// MESH /////////////

// Mesh Trial //
real DD = 1.0;
real WW = 10.0;
// CAVITY Sipp LEBEDEV 2007 :: Trial //

///Thick line object
border upstream1 (t =-0.6,0.0) {x = t; y = 0.0; label = 2;}; //solid thick line
border cavityleft1 (t =0.0,-0.1) {x = 0.0; y = t; label = 2;}; //solid thick line
border cavitybottom1 (t =0.0,1.0) {x=t;y=-0.1;label =2;}; //solid thick line
border cavityright1 (t =-0.1,0.0) {x = 1.0; y = t; label = 2;}; // solid thick line
border downstream1 (t =1.0, 2.50) {x = t; y = 0.0; label = 2;}; //solid thick line
border outlet1 (t =0.0,0.1) {x = 2.5; y = t; label = 3;}; // solid thick line
border top1 (t =2.5,-0.6) {x=t;y=0.1; label = 999;}; //solid thick line
border inlet1 (t =0.1,0.0) {x=-0.6;y=t; label = 1;}; // solid thick line

//thin line bottom
border cavityleft2 (t =-0.1,-0.2) {x = 0.0; y = t; label = 2;}; //solid thin line
border cavitybottom2 (t =0.0,1.0) {x=t;y=-0.20;label =2;}; //solid thin line
border cavityright2 (t =-0.2,-0.1) {x = 1.0; y = t; label = 2;}; //solid thin line

//dashed bottom

border cavityleft3 (t =-0.2,-0.35) {x = 0.0; y = t; label = 2;}; // dashed
border cavitybottom3 (t =0.0,1.0) {x=t;y=-0.35;label = 2;}; //dashed
border cavityright3 (t =-0.35,-0.2) {x = 1.0; y = t; label = 2;}; //dashed

//dotted bottom
border cavityleft4 (t =-0.35,-1.0) {x = 0.0; y = t; label = 2;}; // dotted
border cavitybottom4 (t =0.0,1.0) {x = t; y = -1.0; label = 2;}; //dotted
border cavityright4 (t =-1.0,-0.35) {x = 1.0; y = t; label = 2;}; //dotted

//thin line top
border outlet5 (t =0.1,0.15) {x = 2.5; y = t; label = 3;}; //solid thin line
border top5 (t =2.5,-1.2) {x=t;y=0.15;label = 4;}; //solid thin line
border inlet5 (t =0.15,0.0) {x=-1.2;y=t;label=1;}; // solid thin line
border upstream5 (t =-1.2,-0.6) {x = t; y = 0.0; label = 2;}; //solid thin line

//dashed top
border outlet6 (t =0.15,0.30) {x = 2.5; y=t; label = 3;}; //dashed
border top6 (t =2.5,-1.2) {x=t;y=0.3; label = 4;}; //dashed
border inlet6 (t =0.3,0.15) {x=-1.2; y=t; label=1;}; //dashed

//dotted top
border outlet7 (t = 0.30,0.50) {x = 2.5; y=t; label=3;}; //dotted
border top7 (t = 2.5, -1.2) {x = t; y = 0.5; label = 4;}; //dotted
border inlet7 (t = 0.50,0.30) {x = -1.2; y = t; label = 1;}; //dotted

// For refinement while considering a large domain size

mesh th = buildmesh( upstream1(350)
+ cavityleft1(350)
+ cavitybottom1(350)
+ cavityright1(350)
+ downstream1(350)
+ outlet1(350)
+ top1(350)
+ inlet1(350)

                + cavityleft2(150)
                + cavitybottom2(150)
                + cavityright2(150)

                + cavityleft3(100)
                + cavitybottom3(100)
                + cavityright3(100)

                + cavityleft4(50)
                + cavitybottom4(50)
                + cavityright4(50)

                + outlet5(200)
                + top5(200)
                + inlet5(200)
                + upstream5(200)

                + outlet6(100)
                + top6(100)
                + inlet6(100)

                + outlet7(100)
                + top7(100)
                + inlet7(100)
                
                );

plot(th, value = 1, wait = 1);

//////////////////////////////////////////////////////////////////
//////////// FINITE ELEMENT SPACES ////////////////

/////////*

fespace Uh(th, P2); //velocities
fespace Ph(th, P1); //pressure
fespace UhxUhxPh(th,[P2,P2,P1]); // 2 velocity components and 1 pressure field

//////////////////////////////////////////////////////////////////
////////// NUMERICAL PARAMETERS ////////////////

real u0 = 1.0; // free-stream velocity
real diam = 1.0; // diameter of the cylinder, length scale

real Re = 100.0; // Reynolds number
real dt = 0.01; // time-step
real nu = u0*diam/Re;

int rep=0; // To re-start the simulation using a known flow field
int stepvisu=1000; // To visualise the flow fields
int stepsave=1000; // To save the flow fields
int vis=0; // 0: don’t wait; 1: wait
int iter=0; // if rep neq 0, the value of the iteration corresponding to the starting flow field

////////////////////////////////////////////////////////////////////
/////// INLET CONDITIONS (if needed) //////////

// Thickness of the boundary layer at the inlet
real delta=0.2;

func real uinit(real x, real y)
{
if(y<0)
return 0;
else
return u0*tanh((y)/delta);
};

//////////////////////////////////////////////////////////////////////
/////// VELOCITY FIELDS //////////

UhxUhxPh [ux,uy,p]; //Flow field for the present time step

UhxUhxPh [vx,vy,q]; // Test functions

UhxUhxPh [upx,upy,pp]; // Flow field for the previous time step

UhxUhxPh [rhs1,rhs2,rhs3]; // Right hand side of the NS formed by the terms from the previous time step

//////////////////////////////////////////////////////////////////////
//////// VARIATIONAL FORMULATION ////////////////

// LINEAR TERMS ///
varf NS ([ux,uy,p],[vx,vy,q])=
int2d(th)(
(uxvx+uyvy)/dt
+ nu*( dx(ux)*dx(vx)+dy(ux)*dy(vx)
+ dx(uy)dx(vy)+dy(uy)dy(vy) )
- p( dx(vx) + dy(vy) )
+ q( dx(ux) + dy(uy) ) )

// NO-STRESS CONDITION
+int1d(th,3)(nu*( (N.ydy(ux)+N.xdx(ux))vx)
- ( (N.yvy + N.x*vx )*p) )

  • on(1,uy=0.0,ux=u0)
  • on(2,uy=0.0, ux = 0.0)
  • on(3,uy=0.0)
  • on(4,uy=0.0);
    // Boundary cond.

// NONLINEAR TERMS

// The nonlinear terms: upx upy stand for flow fields from the previous time step //

// you should also go through the manual where they have discussed “varf” in detail…

varf NSRITE ([ux,uy,p],[vx,vy,q])=
int2d(th)( (upxvx+upyvy)/dt
- ( upy*dy(upy)vy + upxdx(upy)vy
+ upydy(upx)vx + upxdx(upx)*vx )
)

+int1d(th,3)(nu*( (N.ydy(ux)+N.xdx(ux))vx)
- ( (N.yvy + N.x*vx )*p) )

// same boundary conditions of course apply to the convective nonlinear terms too…

  • on(1,uy=0.0,ux=u0)
  • on(2,uy=0.0, ux = 0.0)
  • on(3,uy=0.0)
  • on(4,uy=0.0);
    // Boundary cond.

//////////////////////////////////////////////////////////////////////////
///////// INITIAL FLOW FIELDS //////////////////

real tps;

// Start from a scratch solution
tps = 0.0;

if(rep == 0)
{
[ux,uy,p] = [0.0,0.0,0.0]; // good starting flow field in general!
plot(ux,value=1,wait=vis);
plot(uy,value=1,wait=vis);
plot(p,value=1,wait=vis);
}
else
{
// Start from a previous solution
tps=rep*dt;
{
ifstream file(“sol-20-2000_2D_DNS_CYL_dt_0p01_Re_20_mesh1.txt”);
file >> ux[];
};
plot(ux,value=1,wait=vis);
plot(uy,value=1,wait=vis);
plot(p,value=1,wait=vis);
};

////////////////////////////////////////////////////////
//////// MATRIX FORMULATION //////////////

// varf NS ([ux,uy,p],[vx,vy,q])

matrix NSMAT = NS(UhxUhxPh,UhxUhxPh,solver=UMFPACK);

////////////////////////////////////////////////////////
//////// TIME MARCHING //////////////////

for (int i=1; i<=10000; i++)
{

[upx,upy,pp] = [ux,uy,p];

rhs1[] = NSRITE(0,UhxUhxPh);
ux[] = NSMAT^-1*rhs1[];


tps=tps+dt;
iter = iter + 1;

cout << " i = " << i  << "\n";
cout << " tps = " << tps  << "\n";
cout << " UX = " << ux(0.5, 0.5)  << "\n";
cout << " UY = " << uy(0.5, 0.5)  << "\n";

// Vizualise the flow fields
if (!(i%stepvisu))
{
    plot(p,value=1,wait=vis);
    plot(ux,value=1,wait=vis);
    plot(uy,value=1,wait=vis);
  //  plot([UX,tps],value=1,wait=vis);    //
};

// Save the flow fields
if(!(i%stepsave))
{
    ofstream file("sol-"+Re+"-"+iter+"_2D_DNS_CYL_dt_0p01_Re_46p6_Mesh1.txt");
    file << ux[] << endl;
}
///// For time traces at certain points within the domain /////

ofstream uxvel1("ux_2D_DNS_CAV_TEST.dat",append);
ofstream uyvel1("uy_2D_DNS_CAV_TEST.dat",append);


uxvel1<< tps << " " << ux(0.5,0.5) <<endl ;
uyvel1<< tps << " " << uy(0.5,0.5) <<endl ;

};
// End of the time loop


plot(p,value=1,wait=vis);
plot(ux,value=1,wait=vis);
plot(uy,value=1,wait=vis);

///////////////////////////////////////////////
/////// POST-PROCESSING ////////////

fespace PPH(th,P1);
PPH UX, UY, VORT;
UX = ux;
UY = uy;

Uh vort;
vort = dx(uy) - dy(ux);
VORT = vort;

//////

///// MATLAB ///////

include "../../ffmatlib/ffmatlib.idp"

savemesh(th,"th_Mesh1.msh");
ffSaveVh(th,PPH,"fespace_PPH_P1.txt");
ffSaveData(UX,"UX_2D_DNS_CYL_"+iter+"_dt_0p01_Re_46p6_Mesh1.txt");
ffSaveData(UY,"UY_2D_DNS_CYL_"+iter+"_dt_0p01_Re_46p6_Mesh1.txt");
ffSaveData(VORT,"VORT_2D_DNS_CYL_"+iter+"_dt_0p01_Re_46p6_Mesh1.txt");

//////////////////////////
2

Use load "PETSc" + savevtk(..., append = true);, see this example.
It will generate a .pvd file that you can open with ParaView, and you’ll have a movie.

4

What did you write and what is not working?

6

Sorry, I don’t use MATLAB. Again, I’d suggest using ParaView.

1 Like