CoreFlows/examples/C/CoupledTransportDiffusionEquations_1DHeatedChannel.cxx

   1 #include "TransportEquation.hxx"\r
   2 #include "DiffusionEquation.hxx"\r
   3 \r
   4 using namespace std;\r
   5 \r
   6 #define PI 3.14159265\r
   7 \r
   8 void power_field_CoupledTransportDiffusionTest(Field & Phi){\r
   9         double L=4.2;\r
  10         double lambda=0.2;\r
  11         double phi=1e5;\r
  12         double x;\r
  13         Mesh M = Phi.getMesh();\r
  14         int nbCells = M.getNumberOfCells();\r
  15         for (int j = 0; j < nbCells; j++) {\r
  16                 x=M.getCell(j).x();\r
  17                 Phi(j) = phi*cos(PI*(x-L/2)/(L+lambda));\r
  18         }\r
  19 }\r
  20 \r
  21 int main(int argc, char** argv)\r
  22 {\r
  23         //Preprocessing: mesh and group creation\r
  24         double xinf=0.0;\r
  25         double xsup=4.2;\r
  26         int nx=100;\r
  27         double eps=1.E-6;\r
  28         cout << "Building of the diffusion mesh with "<<nx<<" cells" << endl;\r
  29         Mesh diffusionMesh(xinf,xsup,nx);\r
  30         diffusionMesh.setGroupAtPlan(xsup,0,eps,"Neumann");\r
  31         diffusionMesh.setGroupAtPlan(xinf,0,eps,"Neumann");\r
  32 \r
  33         cout << "Building of the transport mesh with "<<nx<<" cells" << endl;\r
  34         Mesh transportMesh(xinf,xsup,nx);\r
  35         transportMesh.setGroupAtPlan(xsup,0,eps,"Neumann");\r
  36         transportMesh.setGroupAtPlan(xinf,0,eps,"Inlet");\r
  37         int spaceDim = 1;\r
  38 \r
  39         // Boundary conditions \r
  40         map<string, LimitFieldDiffusion> boundaryFieldsDiffusion;\r
  41         map<string, LimitFieldTransport> boundaryFieldsTransport;\r
  42 \r
  43         // Boundary conditions for the solid\r
  44         LimitFieldDiffusion limitNeumann;\r
  45         limitNeumann.bcType=NeumannDiffusion;\r
  46         boundaryFieldsDiffusion["Neumann"] = limitNeumann;\r
  47 \r
  48         // Boundary conditions for the fluid\r
  49         LimitFieldTransport limitInlet;\r
  50         limitInlet.bcType=InletTransport;\r
  51         limitInlet.h =1.3e6;//Inlet water enthalpy\r
  52         LimitFieldTransport limitTransportNeumann;\r
  53         limitTransportNeumann.bcType=NeumannTransport;\r
  54         boundaryFieldsTransport["Inlet"] = limitInlet;\r
  55         boundaryFieldsTransport["Neumann"] = limitTransportNeumann;\r
  56 \r
  57         //Set the fluid transport velocity\r
  58         vector<double> transportVelocity(1,5);//Vitesse du fluide\r
  59 \r
  60         //Solid parameters\r
  61         double cp_ur=300;//Uranium specific heat\r
  62         double rho_ur=10000;//Uranium density\r
  63         double lambda_ur=5;\r
  64 \r
  65         TransportEquation  myTransportEquation(LiquidPhase, around155bars600KTransport,transportVelocity);\r
  66 \r
  67         bool FECalculation=false;\r
  68     DiffusionEquation  myDiffusionEquation(spaceDim,FECalculation,rho_ur, cp_ur, lambda_ur);\r
  69 \r
  70 \r
  71         //Set initial field\r
  72         cout << "Construction de la condition initiale " << endl;\r
  73 \r
  74         Vector VV_Constant(1);\r
  75         VV_Constant(0) = 623;//Rod clad temperature nucleaire\r
  76         myDiffusionEquation.setInitialFieldConstant(diffusionMesh,VV_Constant);\r
  77 \r
  78         VV_Constant(0) = 1.3e6;\r
  79         myTransportEquation.setInitialFieldConstant(transportMesh,VV_Constant);\r
  80 \r
  81         Field solidTemp("Solid temperature", CELLS, diffusionMesh, 1);\r
  82         Field fluidTemp("Fluid temperature", CELLS, transportMesh, 1);\r
  83 \r
  84         double heatTransfertCoeff=10000;//fluid/solid heat exchange coefficient\r
  85         myTransportEquation.setHeatTransfertCoeff(heatTransfertCoeff);\r
  86         myDiffusionEquation.setHeatTransfertCoeff(heatTransfertCoeff);\r
  87 \r
  88         //Set heat source in the solid\r
  89         Field Phi("Heat power field", CELLS, diffusionMesh, 1);\r
  90         power_field_CoupledTransportDiffusionTest(Phi);\r
  91         myDiffusionEquation.setHeatPowerField(Phi);\r
  92         Phi.writeVTK("1DheatPowerField");\r
  93 \r
  94         //set the boundary conditions\r
  95         myTransportEquation.setBoundaryFields(boundaryFieldsTransport);//Neumann and Inlet BC will be used\r
  96         myDiffusionEquation.setBoundaryFields(boundaryFieldsDiffusion);//Only Neumann BC will be used\r
  97 \r
  98         // set the numerical method\r
  99         myDiffusionEquation.setTimeScheme( Explicit);\r
 100         myTransportEquation.setTimeScheme( Explicit);\r
 101 \r
 102         // name result file\r
 103         string fluidFileName = "1DFluidEnthalpy";\r
 104         string solidFileName = "1DSolidTemperature";\r
 105 \r
 106         // parameters calculation\r
 107         unsigned MaxNbOfTimeStep =3;\r
 108         int freqSave = 10;\r
 109         double cfl = 0.5;\r
 110         double maxTime = 1000000;\r
 111         double precision = 1e-6;\r
 112 \r
 113         myDiffusionEquation.setCFL(cfl);\r
 114         myDiffusionEquation.setPrecision(precision);\r
 115         myDiffusionEquation.setMaxNbOfTimeStep(MaxNbOfTimeStep);\r
 116         myDiffusionEquation.setTimeMax(maxTime);\r
 117         myDiffusionEquation.setFreqSave(freqSave);\r
 118         myDiffusionEquation.setFileName(solidFileName);\r
 119 \r
 120         myTransportEquation.setCFL(cfl);\r
 121         myTransportEquation.setPrecision(precision);\r
 122         myTransportEquation.setMaxNbOfTimeStep(MaxNbOfTimeStep);\r
 123         myTransportEquation.setTimeMax(maxTime);\r
 124         myTransportEquation.setFreqSave(freqSave);\r
 125         myTransportEquation.setFileName(fluidFileName);\r
 126 \r
 127         // loop on time steps\r
 128         myDiffusionEquation.initialize();\r
 129         myTransportEquation.initialize();\r
 130 \r
 131         double time=0,dt=0;\r
 132         int nbTimeStep=0;\r
 133         bool stop=false, stop_transport=false, stop_diffusion=false; // Does the Problem want to stop (error) ?\r
 134         bool ok; // Is the time interval successfully solved ?\r
 135 \r
 136         // Time step loop\r
 137         while(!stop && !(myDiffusionEquation.isStationary() && myTransportEquation.isStationary()) &&time<maxTime && nbTimeStep<MaxNbOfTimeStep)\r
 138         {\r
 139                 ok=false; // Is the time interval successfully solved ?\r
 140                 fluidTemp=myTransportEquation.getFluidTemperatureField();\r
 141                 solidTemp=myDiffusionEquation.getRodTemperatureField();\r
 142                 myDiffusionEquation.setFluidTemperatureField(fluidTemp);\r
 143                 myTransportEquation.setRodTemperatureField(solidTemp);\r
 144                 // Guess the next time step length\r
 145                 dt=min(myDiffusionEquation.computeTimeStep(stop),myTransportEquation.computeTimeStep(stop));\r
 146                 if (stop){\r
 147                         cout << "Failed computing time step "<<nbTimeStep<<", time = " << time <<", dt= "<<dt<<", stopping calculation"<< endl;\r
 148                 break;\r
 149                 }\r
 150                 // Loop on the time interval tries\r
 151                 while (!ok && !stop )\r
 152                 {\r
 153                         stop_transport=!myTransportEquation.initTimeStep(dt);\r
 154                         stop_diffusion=!myDiffusionEquation.initTimeStep(dt);\r
 155                         stop=stop_diffusion && stop_transport;\r
 156 \r
 157                         // Prepare the next time step\r
 158                         if (stop){\r
 159                                 cout << "Failed initializing time step "<<nbTimeStep<<", time = " << time <<", dt= "<<dt<<", stopping calculation"<< endl;\r
 160                         break;\r
 161                         }\r
 162                         // Solve the next time step\r
 163                         ok=myDiffusionEquation.solveTimeStep()&& myTransportEquation.solveTimeStep();\r
 164 \r
 165                         if (!ok)   // The resolution failed, try with a new time interval.\r
 166                         {\r
 167                                 myDiffusionEquation.abortTimeStep();\r
 168                                 myTransportEquation.abortTimeStep();\r
 169                                         cout << "Failed solving time step "<<nbTimeStep<<", time = " << time<<" dt= "<<dt<<", cfl= "<<cfl <<", stopping calculation"<< endl;\r
 170                                         stop=true; // Impossible to solve the next time step, the Problem has given up\r
 171                                         break;\r
 172                         }\r
 173                         else // The resolution was successful, validate and go to the next time step.\r
 174                         {\r
 175                                 cout << "Time step = "<< nbTimeStep << ", dt = "<< dt <<", time = "<<time << endl;\r
 176                                 myDiffusionEquation.validateTimeStep();\r
 177                                 myTransportEquation.validateTimeStep();\r
 178                                 time=myDiffusionEquation.presentTime();\r
 179                                 nbTimeStep++;\r
 180                         }\r
 181                 }\r
 182         }\r
 183         if(myDiffusionEquation.isStationary() && myTransportEquation.isStationary())\r
 184                 cout << "Stationary state reached" <<endl;\r
 185         else if(time>=maxTime)\r
 186                 cout<<"Maximum time "<<maxTime<<" reached"<<endl;\r
 187         else if(nbTimeStep>=MaxNbOfTimeStep)\r
 188                 cout<<"Maximum number of time steps "<<MaxNbOfTimeStep<<" reached"<<endl;\r
 189         else\r
 190                 cout<<"Error problem wants to stop!"<<endl;\r
 191 \r
 192         cout << "End of calculation time t= " << time << " at time step number "<< nbTimeStep << endl;\r
 193         if (ok)\r
 194                 cout << "Coupled simulation "<<fluidFileName<<" and "<<solidFileName<<" was successful !" << endl;\r
 195         else\r
 196                 cout << "Coupled simulation "<<fluidFileName<<" and "<<solidFileName<<"  failed ! " << endl;\r
 197 \r
 198         cout << "------------ End of calculation -----------" << endl;\r
 199         myDiffusionEquation.terminate();\r
 200         myTransportEquation.terminate();\r
 201 \r
 202         return EXIT_SUCCESS;\r
 203 }\r