CoreFlows/examples/Python/SinglePhase/SinglePhase_2DPoiseuilleFlow.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 import CoreFlows as cf
   5 import cdmath as cm
   6
   7 def SinglePhase_2DPoiseuilleFlow():
   8         spaceDim = 2;
   9
  10         print("Building the mesh" );
  11     # Prepare the mesh data
  12         xinf = 0 ;
  13         xsup = 1.0;
  14         yinf = 0.0;
  15         ysup = 4;
  16         nx   = 10;
  17         ny   = 40;
  18
  19         my_mesh=cm.Mesh(xinf,xsup,nx,yinf,ysup,ny)
  20         # set the boundary names for each boundary
  21         eps=1e-6;
  22         my_mesh.setGroupAtPlan(xsup,0,eps,"wall")
  23         my_mesh.setGroupAtPlan(xinf,0,eps,"wall")
  24         my_mesh.setGroupAtPlan(ysup,1,eps,"neumann")
  25         my_mesh.setGroupAtPlan(yinf,1,eps,"neumann")
  26
  27     # physical constants
  28         viscosity=0.025
  29         viscosite=[viscosity];
  30         Vy_max             = 1.5
  31         a = -8*viscosity*Vy_max / ((ysup-yinf)*(ysup-yinf) ) #pressure slope
  32
  33         initialTemperature = 573
  34         outletPressure     = 155e5
  35
  36         initial_field=cm.Field("Initial field", cm.CELLS, my_mesh, 4)
  37         for i in range( 0 , my_mesh.getNumberOfCells() ):
  38                 Ci=my_mesh.getCell(i)
  39                 x=Ci.x()
  40                 y=Ci.y()
  41                 initial_field[i,0] =  outletPressure + a*(y - ysup )
  42                 initial_field[i,1] =  0  #x component of the velocity
  43                 initial_field[i,2] =  a/(2*viscosity)*( (x-(xsup+xinf)/2)*(x-(xsup+xinf)/2) - (xsup-xinf)*(xsup-xinf)/4)  #y component of the velocity
  44                 initial_field[i,3] =  initialTemperature
  45
  46
  47     # set the limit field for each boundary
  48         wallVelocityX=0;
  49         wallVelocityY=0;
  50         wallTemperature=573;
  51
  52         myProblem = cf.SinglePhase(cf.Liquid,cf.around155bars600K,spaceDim);
  53         nVar =myProblem.getNumberOfVariables();
  54
  55     #Initial field creation
  56         print("Setting initial data" );
  57         myProblem.setInitialField(initial_field)
  58
  59     # the boundary conditions
  60         myProblem.setWallBoundaryCondition("wall", wallTemperature, wallVelocityX, wallVelocityY);
  61         myProblem.setNeumannBoundaryCondition("neumann");
  62
  63     # set physical parameters
  64         myProblem.setViscosity(viscosite);
  65
  66         # set the numerical method
  67         myProblem.setNumericalScheme(cf.upwind, cf.Implicit);
  68         #myProblem.setLinearSolver(cf.GMRES,cf.LU,True);
  69
  70         # name file save
  71         fileName = "2DPoiseuilleFlow";
  72
  73         # parameters calculation
  74         MaxNbOfTimeStep = 10000 ;
  75         freqSave = 1;
  76         cfl = 0.5;
  77         maxTime = 5000;
  78         precision = 1e-6;
  79
  80         myProblem.setCFL(cfl);
  81         myProblem.setPrecision(precision);
  82         myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
  83         myProblem.setTimeMax(maxTime);
  84         myProblem.setNewtonSolver(float('inf'),20);
  85         myProblem.setFreqSave(freqSave);
  86         myProblem.setFileName(fileName);
  87         myProblem.saveConservativeField(True);
  88         if(spaceDim>1):
  89                 myProblem.saveVelocity();
  90                 pass
  91
  92         # evolution
  93         myProblem.initialize();
  94
  95         ok = myProblem.run();
  96         if (ok):
  97                 print( "Simulation python " + fileName + " is successful !" );
  98                 pass
  99         else:
 100                 print( "Simulation python " + fileName + "  failed ! " );
 101                 pass
 102
 103         print( "------------ End of calculation !!! -----------" );
 104
 105         myProblem.terminate();
 106         return ok
 107
 108 if __name__ == """__main__""":
 109     SinglePhase_2DPoiseuilleFlow()