7 #===============================================================================================================================
8 # Name : Simulation of a 2D heat equation
9 # Description : Test solving the diffusion of the temperature T in a solid
10 # \rho cp dT/dt-\lambda\Delta T=\Phi + \lambda_{sf} (T_{fluid}-T_{solid})
11 # Neumann or Dirichlet boundary conditions
12 # Finite elements or finite volumes
13 # Author : Michaël Ndjinga
14 # Copyright : CEA Saclay 2021
15 #================================================================================================================================
18 def DiffusionEquation_2DSpherical(FECalculation, fileName):
20 """ Description : Test solving the diffusion of the temperature T in a solid (default is Uranium).
21 Equation : Thermal diffusion equation \rho cp dT/dt-\lambda\Delta T=\Phi + \lambda_{sf} (T_{fluid}-T_{solid})
22 Heat capacity cp, density \rho, and conductivity \lambda of the solid MUST be defined
23 The solid may be refrigerated by a fluid with temperature T_{solid} transfer coefficient \lambda_{sf} using functions setFluidTemperature and setHeatTransfertCoeff
24 The solid may receive some extra heat power due to nuclear fissions using function setHeatSource
26 #Space dimension of the problem
29 # Mandatory physical values
30 solid_specific_heat=300# specific heat capacity, default value 300
31 solid_density=10000# density, default value 10000
32 solid_conductivity=5# conductivity, default value 5
34 myProblem = solverlab.DiffusionEquation(spaceDim,FECalculation,solid_density,solid_specific_heat,solid_conductivity);
36 # Definition of field support parameter
38 supportOfField=solverlab.NODES
40 supportOfField=solverlab.CELLS
42 # Set the mesh and initial data
43 initial_data_inputfile="../resources/meshSquare";
44 initial_data_fieldName="Solid temperature";
45 print("Loading unstructured mesh and initial data", " in file ", initial_data_inputfile )
46 initial_data_time_iteration=0# default value is 0
47 initial_data_time_sub_iteration=0# default value is 0
48 initial_data_time_meshLevel=0# default value is 0
49 myProblem.setInitialField(initial_data_inputfile, initial_data_fieldName, initial_data_time_iteration, initial_data_time_sub_iteration, initial_data_time_meshLevel, supportOfField)
51 #### Optional physical values (default value is zero) : fluid temperature field, heat transfert coefficient, heat power field
52 # Loading and setting fluid temperature field
53 fluid_temperature_inputfile="../resources/meshSquare";
54 fluid_temperature_fieldName="Fluid temperature";
55 fluid_temperature_time_iteration=0# default value is 0
56 fluid_temperature_time_sub_iteration=0# default value is 0
57 fluid_temperature_meshLevel=0# default value is 0
58 print("Loading field :", fluid_temperature_fieldName, " in file ", fluid_temperature_inputfile)
59 fluidTemperatureField=solverlab.Field(fluid_temperature_inputfile, supportOfField, fluid_temperature_fieldName, fluid_temperature_time_iteration, fluid_temperature_time_sub_iteration, fluid_temperature_meshLevel)
60 myProblem.setFluidTemperatureField(fluidTemperatureField)
61 # Setting heat transfert coefficient
62 heatTransfertCoeff=1000.;#fluid/solid exchange coefficient, default value is 0
63 myProblem.setHeatTransfertCoeff(heatTransfertCoeff);
64 # Loading heat power field
65 heat_power_inputfile="../resources/meshSquare";
66 heat_power_fieldName="Heat power";
67 heat_power_time_iteration=0# default value is 0
68 heat_power_time_sub_iteration=0# default value is 0
69 heat_power_meshLevel=0# default value is 0
70 print("Loading field :", heat_power_fieldName, " in file ", heat_power_inputfile)
71 heatPowerField=solverlab.Field(heat_power_inputfile, supportOfField, heat_power_fieldName, heat_power_time_iteration, heat_power_time_sub_iteration, heat_power_meshLevel)
72 myProblem.setHeatPowerField(heatPowerField)
74 # the boundary conditions :
76 boundaryGroupNames=myProblem.getMesh().getNameOfNodeGroups()
77 print(len(boundaryGroupNames), " Boundary Node Group detected : ", boundaryGroupNames)
79 boundaryGroupNames=myProblem.getMesh().getNameOfFaceGroups()
80 print(len(boundaryGroupNames), " Boundary Face Group detected : ", boundaryGroupNames)
82 # for each boundary we load the boundary field (replace by a loop over the boundaries)
83 boundary1_type=solverlab.NeumannDiffusion
84 boundary1_inputfile="../resources/meshSquare";
85 boundary1_fieldName="Left temperature";
86 boundary1_time_iteration=0# default value is 0
87 boundary1_time_sub_iteration=0# default value is 0
88 boundary1_meshLevel=0# default value is 0
89 print("Boundary ", boundaryGroupNames[3], ", loading field :", boundary1_fieldName, " in file ", boundary1_inputfile)
90 boundary1Field=solverlab.Field(boundary1_inputfile, supportOfField, boundary1_fieldName, boundary1_time_iteration, boundary1_time_sub_iteration, boundary1_meshLevel)
91 boundary2_type=solverlab.DirichletDiffusion
92 boundary2_inputfile="../resources/meshSquare";
93 boundary2_fieldName="Right temperature";
94 boundary2_time_iteration=0# default value is 0
95 boundary2_time_sub_iteration=0# default value is 0
96 boundary2_meshLevel=0# default value is 0
97 print("Boundary ", boundaryGroupNames[2], ", loading field :", boundary2_fieldName, " in file ", boundary2_inputfile)
98 boundary2Field=solverlab.Field(boundary2_inputfile, supportOfField, boundary2_fieldName, boundary2_time_iteration, boundary2_time_sub_iteration, boundary2_meshLevel)
99 boundary3_type=solverlab.NeumannDiffusion
100 boundary3_inputfile="../resources/meshSquare";
101 boundary3_fieldName="Top temperature";
102 boundary3_time_iteration=0# default value is 0
103 boundary3_time_sub_iteration=0# default value is 0
104 boundary3_meshLevel=0# default value is 0
105 print("Boundary ", boundaryGroupNames[4], ", loading field :", boundary3_fieldName, " in file ", boundary3_inputfile)
106 boundary3Field=solverlab.Field(boundary3_inputfile, supportOfField, boundary3_fieldName, boundary3_time_iteration, boundary3_time_sub_iteration, boundary3_meshLevel)
107 boundary4_type=solverlab.DirichletDiffusion
108 boundary4_inputfile="../resources/meshSquare";
109 boundary4_fieldName="Bottom temperature";
110 boundary4_time_iteration=0# default value is 0
111 boundary4_time_sub_iteration=0# default value is 0
112 boundary4_meshLevel=0# default value is 0
113 print("Boundary ", boundaryGroupNames[1], ", loading field :", boundary4_fieldName, " in file ", boundary4_inputfile)
114 boundary4Field=solverlab.Field(boundary4_inputfile, supportOfField, boundary4_fieldName, boundary4_time_iteration, boundary4_time_sub_iteration, boundary4_meshLevel)
116 # for each boundary we need to know if we want a Neumann or a Dirichlet boundary condition
117 if boundary1_type==solverlab.NeumannDiffusion :
118 myProblem.setNeumannBoundaryCondition("Left", boundary1Field)
119 elif boundary1_type==solverlab.DirichletDiffusion :
120 myProblem.setDirichletBoundaryCondition("Left", boundary1Field)
121 if boundary2_type==solverlab.NeumannDiffusion :
122 myProblem.setNeumannBoundaryCondition("Right", boundary2Field)
123 elif boundary2_type==solverlab.DirichletDiffusion :
124 myProblem.setDirichletBoundaryCondition("Right", boundary2Field)
125 if boundary3_type==solverlab.NeumannDiffusion :
126 myProblem.setNeumannBoundaryCondition("Top", boundary3Field)
127 elif boundary3_type==solverlab.DirichletDiffusion :
128 myProblem.setDirichletBoundaryCondition("Top", boundary3Field);
129 if boundary4_type==solverlab.NeumannDiffusion :
130 myProblem.setNeumannBoundaryCondition("Bottom", boundary4Field)
131 elif boundary4_type==solverlab.DirichletDiffusion :
132 myProblem.setDirichletBoundaryCondition("Bottom", boundary4Field);
134 # set the numerical method
135 myProblem.setTimeScheme( solverlab.Explicit)# Otherwise solverlab.Implicit
136 max_nb_its_lin_solver = 50
137 myProblem.setLinearSolver(solverlab.GMRES, solverlab.ILU, max_nb_its_lin_solver );
139 # computation parameters
140 MaxNbOfTimeStep = 3 ;# default value is 10
141 freqSave = 1;# default value is 1
142 cfl = 0.95;# default value is 1
143 maxTime = 100000000;# default value is 10
144 precision = 1e-6;# default value is 1e-6
145 result_directory="."# default value = current directory
147 myProblem.setCFL(cfl);
148 myProblem.setPrecision(precision);
149 myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
150 myProblem.setTimeMax(maxTime);
151 myProblem.setFreqSave(freqSave);
152 myProblem.setFileName(fileName);
153 myProblem.setResultDirectory(result_directory)
154 myProblem.setSaveFileFormat(solverlab.MED)#default value is solverlab.VTK
157 myProblem.initialize();
158 print("Running python test "+ fileName );
160 ok = myProblem.run();
162 print( "Python simulation " + fileName + " is successful !" );
165 print( "Python simulation " + fileName + " failed ! " );
168 print( "------------ End of simulation !!! -----------" );
170 myProblem.terminate();
173 if __name__ == """__main__""":
174 if len(sys.argv) >1 :
175 FECalculation = bool(int(sys.argv[1]))
176 # name of result file
178 fileName = "2DSpherical_FE";# default value is ""
180 fileName = "2DSpherical_FV";# default value is ""
181 DiffusionEquation_2DSpherical(FECalculation, fileName)
183 raise ValueError("DiffusionEquation_2DSpherical : missing one argument")