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