1 #include "TransportEquation.hxx"
\r
2 #include "DiffusionEquation.hxx"
\r
6 #define PI 3.14159265
\r
8 void power_field_CoupledTransportDiffusionTest(Field & Phi){
\r
13 Mesh M = Phi.getMesh();
\r
14 int nbCells = M.getNumberOfCells();
\r
15 for (int j = 0; j < nbCells; j++) {
\r
17 Phi(j) = phi*cos(PI*(x-L/2)/(L+lambda));
\r
21 int main(int argc, char** argv)
\r
23 //Preprocessing: mesh and group creation
\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
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
39 // Boundary conditions
\r
40 map<string, LimitFieldDiffusion> boundaryFieldsDiffusion;
\r
41 map<string, LimitFieldTransport> boundaryFieldsTransport;
\r
43 // Boundary conditions for the solid
\r
44 LimitFieldDiffusion limitNeumann;
\r
45 limitNeumann.bcType=NeumannDiffusion;
\r
46 boundaryFieldsDiffusion["Neumann"] = limitNeumann;
\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
57 //Set the fluid transport velocity
\r
58 vector<double> transportVelocity(1,5);//Vitesse du fluide
\r
61 double cp_ur=300;//Uranium specific heat
\r
62 double rho_ur=10000;//Uranium density
\r
65 TransportEquation myTransportEquation(Water, around155bars600K,transportVelocity);
\r
67 bool FECalculation=false;
\r
68 DiffusionEquation myDiffusionEquation(spaceDim,FECalculation,rho_ur, cp_ur, lambda_ur);
\r
72 cout << "Construction de la condition initiale " << endl;
\r
74 Vector VV_Constant(1);
\r
75 VV_Constant(0) = 623;//Rod clad temperature nucleaire
\r
76 myDiffusionEquation.setInitialFieldConstant(diffusionMesh,VV_Constant);
\r
78 VV_Constant(0) = 1.3e6;
\r
79 myTransportEquation.setInitialFieldConstant(transportMesh,VV_Constant);
\r
81 Field solidTemp("Solid temperature", CELLS, diffusionMesh, 1);
\r
82 Field fluidTemp("Fluid temperature", CELLS, transportMesh, 1);
\r
84 double heatTransfertCoeff=10000;//fluid/solid heat exchange coefficient
\r
85 myTransportEquation.setHeatTransfertCoeff(heatTransfertCoeff);
\r
86 myDiffusionEquation.setHeatTransfertCoeff(heatTransfertCoeff);
\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
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
98 // set the numerical method
\r
99 myDiffusionEquation.setTimeScheme( Explicit);
\r
100 myTransportEquation.setTimeScheme( Explicit);
\r
102 // name result file
\r
103 string fluidFileName = "1DFluidEnthalpy";
\r
104 string solidFileName = "1DSolidTemperature";
\r
106 // parameters calculation
\r
107 unsigned MaxNbOfTimeStep =3;
\r
110 double maxTime = 1000000;
\r
111 double precision = 1e-6;
\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
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
127 // loop on time steps
\r
128 myDiffusionEquation.initialize();
\r
129 myTransportEquation.initialize();
\r
131 double time=0,dt=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
137 while(!stop && !(myDiffusionEquation.isStationary() && myTransportEquation.isStationary()) &&time<maxTime && nbTimeStep<MaxNbOfTimeStep)
\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
147 cout << "Failed computing time step "<<nbTimeStep<<", time = " << time <<", dt= "<<dt<<", stopping calculation"<< endl;
\r
150 // Loop on the time interval tries
\r
151 while (!ok && !stop )
\r
153 stop_transport=!myTransportEquation.initTimeStep(dt);
\r
154 stop_diffusion=!myDiffusionEquation.initTimeStep(dt);
\r
155 stop=stop_diffusion && stop_transport;
\r
157 // Prepare the next time step
\r
159 cout << "Failed initializing time step "<<nbTimeStep<<", time = " << time <<", dt= "<<dt<<", stopping calculation"<< endl;
\r
162 // Solve the next time step
\r
163 ok=myDiffusionEquation.solveTimeStep()&& myTransportEquation.solveTimeStep();
\r
165 if (!ok) // The resolution failed, try with a new time interval.
\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
173 else // The resolution was successful, validate and go to the next time step.
\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
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
190 cout<<"Error problem wants to stop!"<<endl;
\r
192 cout << "End of calculation time t= " << time << " at time step number "<< nbTimeStep << endl;
\r
194 cout << "Coupled simulation "<<fluidFileName<<" and "<<solidFileName<<" was successful !" << endl;
\r
196 cout << "Coupled simulation "<<fluidFileName<<" and "<<solidFileName<<" failed ! " << endl;
\r
198 cout << "------------ End of calculation -----------" << endl;
\r
199 myDiffusionEquation.terminate();
\r
200 myTransportEquation.terminate();
\r
202 return EXIT_SUCCESS;
\r