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         Field fluidEnthalpy("Enthalpie", CELLS, transportMesh, 1);\r
  67         bool FECalculation=false;\r
  68     DiffusionEquation  myDiffusionEquation(spaceDim,FECalculation,rho_ur, cp_ur, lambda_ur);\r
  69 \r
  70         Field solidTemp("Solid temperature", CELLS, diffusionMesh, 1);\r
  71         Field fluidTemp("Fluid temperature", CELLS, transportMesh, 1);\r
  72 \r
  73         double heatTransfertCoeff=10000;//fluid/solid heat exchange coefficient\r
  74         myTransportEquation.setHeatTransfertCoeff(heatTransfertCoeff);\r
  75         myDiffusionEquation.setHeatTransfertCoeff(heatTransfertCoeff);\r
  76 \r
  77         //Set heat source in the solid\r
  78         Field Phi("Heat power field", CELLS, diffusionMesh, 1);\r
  79         power_field_CoupledTransportDiffusionTest(Phi);\r
  80         myDiffusionEquation.setHeatPowerField(Phi);\r
  81         Phi.writeVTK("1DheatPowerField");\r
  82 \r
  83         //Initial field creation\r
  84         Vector VV_Constant(1);\r
  85         VV_Constant(0) = 623;//Rod clad temperature nucleaire\r
  86 \r
  87         cout << "Construction de la condition initiale " << endl;\r
  88         // generate initial condition\r
  89         myDiffusionEquation.setInitialFieldConstant(diffusionMesh,VV_Constant);\r
  90 \r
  91 \r
  92         VV_Constant(0) = 1.3e6;\r
  93         myTransportEquation.setInitialFieldConstant(transportMesh,VV_Constant);\r
  94 \r
  95         //set the boundary conditions\r
  96         myTransportEquation.setBoundaryFields(boundaryFieldsTransport);//Neumann and Inlet BC will be used\r
  97         myDiffusionEquation.setBoundaryFields(boundaryFieldsDiffusion);//Only Neumann BC will be used\r
  98 \r
  99         // set the numerical method\r
 100         myDiffusionEquation.setTimeScheme( Explicit);\r
 101         myTransportEquation.setTimeScheme( Explicit);\r
 102 \r
 103         // name result file\r
 104         string fluidFileName = "1DFluidEnthalpy";\r
 105         string solidFileName = "1DSolidTemperature";\r
 106 \r
 107         // parameters calculation\r
 108         unsigned MaxNbOfTimeStep =3;\r
 109         int freqSave = 10;\r
 110         double cfl = 0.5;\r
 111         double maxTime = 1000000;\r
 112         double precision = 1e-6;\r
 113 \r
 114         myDiffusionEquation.setCFL(cfl);\r
 115         myDiffusionEquation.setPrecision(precision);\r
 116         myDiffusionEquation.setMaxNbOfTimeStep(MaxNbOfTimeStep);\r
 117         myDiffusionEquation.setTimeMax(maxTime);\r
 118         myDiffusionEquation.setFreqSave(freqSave);\r
 119         myDiffusionEquation.setFileName(solidFileName);\r
 120 \r
 121         myTransportEquation.setCFL(cfl);\r
 122         myTransportEquation.setPrecision(precision);\r
 123         myTransportEquation.setMaxNbOfTimeStep(MaxNbOfTimeStep);\r
 124         myTransportEquation.setTimeMax(maxTime);\r
 125         myTransportEquation.setFreqSave(freqSave);\r
 126         myTransportEquation.setFileName(fluidFileName);\r
 127 \r
 128         // loop on time steps\r
 129         myDiffusionEquation.initialize();\r
 130         myTransportEquation.initialize();\r
 131 \r
 132         double time=0,dt=0;\r
 133         int nbTimeStep=0;\r
 134         bool stop=false, stop_transport=false, stop_diffusion=false; // Does the Problem want to stop (error) ?\r
 135         bool ok; // Is the time interval successfully solved ?\r
 136 \r
 137         // Time step loop\r
 138         while(!stop && !(myDiffusionEquation.isStationary() && myTransportEquation.isStationary()) &&time<maxTime && nbTimeStep<MaxNbOfTimeStep)\r
 139         {\r
 140                 ok=false; // Is the time interval successfully solved ?\r
 141                 fluidTemp=myTransportEquation.getFluidTemperatureField();\r
 142                 solidTemp=myDiffusionEquation.getRodTemperatureField();\r
 143                 myDiffusionEquation.setFluidTemperatureField(fluidTemp);\r
 144                 myTransportEquation.setRodTemperatureField(solidTemp);\r
 145                 // Guess the next time step length\r
 146                 dt=min(myDiffusionEquation.computeTimeStep(stop),myTransportEquation.computeTimeStep(stop));\r
 147                 if (stop){\r
 148                         cout << "Failed computing time step "<<nbTimeStep<<", time = " << time <<", dt= "<<dt<<", stopping calculation"<< endl;\r
 149                 break;\r
 150                 }\r
 151                 // Loop on the time interval tries\r
 152                 while (!ok && !stop )\r
 153                 {\r
 154                         stop_transport=!myTransportEquation.initTimeStep(dt);\r
 155                         stop_diffusion=!myDiffusionEquation.initTimeStep(dt);\r
 156                         stop=stop_diffusion && stop_transport;\r
 157 \r
 158                         // Prepare the next time step\r
 159                         if (stop){\r
 160                                 cout << "Failed initializing time step "<<nbTimeStep<<", time = " << time <<", dt= "<<dt<<", stopping calculation"<< endl;\r
 161                         break;\r
 162                         }\r
 163                         // Solve the next time step\r
 164                         ok=myDiffusionEquation.solveTimeStep()&& myTransportEquation.solveTimeStep();\r
 165 \r
 166                         if (!ok)   // The resolution failed, try with a new time interval.\r
 167                         {\r
 168                                 myDiffusionEquation.abortTimeStep();\r
 169                                 myTransportEquation.abortTimeStep();\r
 170                                         cout << "Failed solving time step "<<nbTimeStep<<", time = " << time<<" dt= "<<dt<<", cfl= "<<cfl <<", stopping calculation"<< endl;\r
 171                                         stop=true; // Impossible to solve the next time step, the Problem has given up\r
 172                                         break;\r
 173                         }\r
 174                         else // The resolution was successful, validate and go to the next time step.\r
 175                         {\r
 176                                 cout << "Time step = "<< nbTimeStep << ", dt = "<< dt <<", time = "<<time << endl;\r
 177                                 myDiffusionEquation.validateTimeStep();\r
 178                                 myTransportEquation.validateTimeStep();\r
 179                                 time=myDiffusionEquation.presentTime();\r
 180                                 nbTimeStep++;\r
 181                         }\r
 182                 }\r
 183         }\r
 184         if(myDiffusionEquation.isStationary() && myTransportEquation.isStationary())\r
 185                 cout << "Stationary state reached" <<endl;\r
 186         else if(time>=maxTime)\r
 187                 cout<<"Maximum time "<<maxTime<<" reached"<<endl;\r
 188         else if(nbTimeStep>=MaxNbOfTimeStep)\r
 189                 cout<<"Maximum number of time steps "<<MaxNbOfTimeStep<<" reached"<<endl;\r
 190         else\r
 191                 cout<<"Error problem wants to stop!"<<endl;\r
 192 \r
 193         cout << "End of calculation time t= " << time << " at time step number "<< nbTimeStep << endl;\r
 194         if (ok)\r
 195                 cout << "Coupled simulation "<<fluidFileName<<" and "<<solidFileName<<" was successful !" << endl;\r
 196         else\r
 197                 cout << "Coupled simulation "<<fluidFileName<<" and "<<solidFileName<<"  failed ! " << endl;\r
 198 \r
 199         cout << "------------ End of calculation -----------" << endl;\r
 200         myDiffusionEquation.terminate();\r
 201         myTransportEquation.terminate();\r
 202 \r
 203         return EXIT_SUCCESS;\r
 204 }\r