nonlinear_heat_transfer.py 6.27 KB
 Jeremy BLEYER committed Mar 17, 2020 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 #!/usr/bin/env python # coding: utf-8 # # Stationnary non-linear heat transfer # # ## Description of the non-linear constitutive heat transfer law # # ## FEniCS implementation # # We consider a rectanglar domain with imposed temperatures Tl (resp. Tr) on the left (resp. right boundaries). We want to solve for the temperature field T inside the domain using a $P^1$-interpolation. We initialize the temperature at value Tl throughout the domain. We finally load the material library with a plane_strain hypothesis. # In[1]: import matplotlib.pyplot as plt from dolfin import * import mfront_wrapper as mf length = 30e-3 width = 5.4e-3 mesh = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10) V = FunctionSpace(mesh, "CG", 1) T = Function(V, name="Temperature") def left(x, on_boundary): return near(x[0], 0) and on_boundary def right(x, on_boundary): return near(x[0], length) and on_boundary Tl = 300 Tr = 800 T.interpolate(Constant(Tl)) bc = [DirichletBC(V, Constant(Tl), left), DirichletBC(V, Constant(Tr), right)] material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so", "StationaryHeatTransfer", hypothesis="plane_strain") # The MFront behaviour implicitly declares the temperature as an external state variable called "Temperature". We must therefore associate this external state variable to a known mechanical field. This can be achieved explicitly using the register_external_state_variable method. In the present case, this can be donc automatically since the name of the unknown temperature field matches the [TFEL Glossary](http://tfel.sourceforge.net/glossary.html) name "Temperature". In this case, the following message is printed: #  # Automatic registration of 'Temperature' as an external state variable. #  # For problems in which the temperature only acts as a parameter (no jacobian blocks with respect to the temperature), the temperature can be automatically registered as a constant value ($293.15 \text{ K}$ by default) or to any other (dolfin.Constant or float) value using the register_external_state_variable method. # # In the FEniCS interface, we instantiate the main mechanical unknown, here the temperature field T which has to be named "Temperature" in order to match MFront's predefined name. Using another name than this will later result in an error saying: #  # ValueError: 'Temperature' could not be associated with a registered gradient or a known state variable. #  # # The MFront behaviour declares the field "TemperatureGradient" as a Gradient variable, with its associated Flux called "HeatFlux". We can check that the material object retrieves MFront's gradient and flux names, as well as the different tangent operator blocks which have been defined, namely dj_ddgT and dj_ddT in the present case: # In[2]: print(material.get_gradient_names()) print(material.get_flux_names()) print(["d{}_d{}".format(*t) for t in material.get_tangent_block_names()]) # When defining the non-linear problem, we will specify the boundary conditions and the requested quadrature degree which will control the number of quadrature points used in each cell to compute the non-linear constitutive law. Here, we specify a quadrature of degree 2 (i.e. 3 Gauss points for a triangular element). Finally, we need to associate to MFront gradient object the corresponding UFL expression as a function of the unknown field T. To do so, we use the register_gradient method linking MFront "TemperatureGradient" object to the UFL expression grad(T). Doing so, the corresponding non-linear variational problem will be automatically be built: # # # F(\widehat{T}) = \int_\Omega \boldsymbol{j}\cdot \nabla \widehat{T} \text{dx} = 0 \quad \forall \widehat{T} # # In[3]: problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=2, bcs=bc) problem.register_gradient("TemperatureGradient", grad(T)) # From the two tangent operator blocks dj_ddgT and dj_ddT, it will automatically be deduced that the heat flux $\boldsymbol{j}$ is a function of both the temperature gradient $\boldsymbol{g}=\nabla T$ and the temperature itself i.e. $\boldsymbol{j}=\boldsymbol{j}(\boldsymbol{g}, T)$. The following tangent bilinear form will therefore be used when solving the above non-linear problem: # # # J(\widehat{T},T^*) = \int_{\Omega} \nabla \widehat{T}\cdot\left(\dfrac{\partial \boldsymbol{j}}{\partial \boldsymbol{g}}\cdot \nabla T^*+\dfrac{\partial \boldsymbol{j}}{\partial T}\cdot T^*\right) \text{dx} # # # Similarly to the case of external state variables, common gradient expressions for some [TFEL Glossary](http://tfel.sourceforge.net/glossary.html) names have been already predefined which avoid calling explicitly the register_gradient method. Predefined expressions can be obtained from: # In[4]: mf.list_predefined_gradients() # We can see that the name "Temperature Gradient" is in fact a predefined gradient. Omitting calling the register_gradient method will in this case print the following message upon calling solve: #  # Automatic registration of 'TemperatureGradient' as grad(Temperature). #  # meaning that a predefined gradient name has been found and registered as the UFL expression $\nabla T$. # # We finally solve the non-linear problem using a default Newton non-linear solver. The solve method returns the number of Newton iterations (4 in the present case) and converged status . # In[5]: problem.solve(T.vector()) # We finally check that the thermal conductivity coefficient $k$, computed from the ratio between the horizontal heat flux and temperature gradient matches the temperature-dependent expressions implemented in the MFront behaviour. # In[7]: j = problem.fluxes["HeatFlux"].function g = problem.gradients["TemperatureGradient"].function k_gauss = j.vector().get_local()[::2]/g.vector().get_local()[::2] T_gauss = problem.state_variables["external"]["Temperature"].function.vector().get_local() A = 0.0375; B = 2.165e-4; k_ref = 1/(A + B*T_gauss) plt.plot(T_gauss, k_gauss, 'o', label="FE") plt.plot(T_gauss, k_ref, '.', label="ref") plt.xlabel(r"Temperature $T\: (K)$") plt.ylabel(r"Thermal conductivity $k\: (W.m^{-1}.K^{-1})$") plt.legend() plt.show() # In[ ]: