diff --git a/docs/assets/T4NAFEMS--T3-solution.png b/docs/assets/T4NAFEMS--T3-solution.png new file mode 100644 index 0000000..6ea1716 Binary files /dev/null and b/docs/assets/T4NAFEMS--T3-solution.png differ diff --git a/docs/src/guide/guide.md b/docs/src/guide/guide.md index b8342c6..07c8976 100644 --- a/docs/src/guide/guide.md +++ b/docs/src/guide/guide.md @@ -39,4 +39,7 @@ Algorithms typically (not always) accept a single argument, `modeldata`, a dicti There is an implementation of an algorithm for steady-state heat conduction. +## Example + +![Alt Visualization of the temperature field](T4NAFEMS--T3-solution.png "Temperature field raised as a surface above the mesh") diff --git a/examples/steady_state/2-d/T4NAFEMS_examples.jl b/examples/steady_state/2-d/T4NAFEMS_examples.jl index 45136fe..977299a 100644 --- a/examples/steady_state/2-d/T4NAFEMS_examples.jl +++ b/examples/steady_state/2-d/T4NAFEMS_examples.jl @@ -5,10 +5,9 @@ using FinEtoolsHeatDiff.AlgoHeatDiffModule using FinEtools.AlgoBaseModule: richextrapol function T4NAFEMS_T3_algo() - ## Two-dimensional heat transfer with convection: convergence study - # - - ## Description + # ## Two-dimensional heat transfer with convection: convergence study + + # ### Description # # Consider a plate of uniform thickness, measuring 0.6 m by 1.0 m. On one # short edge the temperature is fixed at 100 °C, and on one long edge the @@ -19,79 +18,71 @@ function T4NAFEMS_T3_algo() # There is no internal generation of heat. Calculate the temperature 0.2 m # along the un-insulated long side, measured from the intersection with the # fixed temperature side. The reference result is 18.25 °C. - - ## + # # The reference temperature at the point A is 18.25 °C according to the - # NAFEMS publication ( hich cites the book Carslaw, H.S. and J.C. Jaeger, + # NAFEMS publication (which cites the book Carslaw, H.S. and J.C. Jaeger, # Conduction of Heat in Solids. 1959: Oxford University Press). - - ## + # # The present tutorial will investigate the reference temperature and it # will attempt to estimate the limit value more precisely using a # sequence of meshes and Richardson's extrapolation. - ## Solution - # - + # ### Solution + println(""" NAFEMS benchmark. Two-dimensional heat transfer with convection: convergence study. Solution with linear triangles. - Version: 05/29/2017 + Version: 02/23/2024 """) - - kappa = [52.0 0; 0 52.0] * phun("W/(M*K)") # conductivity matrix - h = 750 * phun("W/(M^2*K)")# surface heat transfer coefficient - Width = 0.6 * phun("M")# Geometrical dimensions + # Conductivity matrix + kappa = [52.0 0; 0 52.0] * phun("W/(M*K)") + # Surface heat transfer coefficient + h = 750 * phun("W/(M^2*K)") + # Geometrical dimensions + Width = 0.6 * phun("M") Height = 1.0 * phun("M") HeightA = 0.2 * phun("M") Thickness = 0.1 * phun("M") tolerance = Width / 1000 + # Create a material model. m = MatHeatDiff(kappa) + # Five progressively refined models will be created and solved. modeldata = nothing - resultsTempA = FFlt[] - for nref = 1:5 - t0 = time() - - # The mesh is created from two triangles to begin with + resultsTempA = Float64[] + params = Float64[] + for nref in 2:6 + # The mesh is created from two rectangular blocks to begin with. fens, fes = T3blockx([0.0, Width], [0.0, HeightA]) fens2, fes2 = T3blockx([0.0, Width], [HeightA, Height]) + # The meshes are then glued into a single entity. fens, newfes1, fes2 = mergemeshes(fens, fes, fens2, fes2, tolerance) fes = cat(newfes1, fes2) - # Refine the mesh desired number of times - for ref = 1:nref + # Refine the mesh desired number of times. + for ref in 1:nref fens, fes = T3refine(fens, fes) end + # The boundary is extracted. bfes = meshboundary(fes) - - # Define boundary conditions - - ## # The prescribed temperature is applied along edge 1 (the bottom - # edge in Figure 1).. - - l1 = selectnode(fens; box = [0.0 Width 0.0 0.0], inflate = tolerance) - essential1 = FDataDict("node_list" => l1, "temperature" => 100.0) - - ## - # The convection boundary condition is applied along the edges - # 2,3,4. The elements along the boundary are quadratic line - # elements L3. The order-four Gauss quadrature is sufficiently - # accurate. - l2 = selectelem(fens, bfes; box = [Width Width 0.0 Height], inflate = tolerance) - l3 = selectelem(fens, bfes; box = [0.0 Width Height Height], inflate = tolerance) + # edge in Figure 1). + list1 = selectnode(fens; box = [0.0 Width 0.0 0.0], inflate = tolerance) + essential1 = FDataDict("node_list" => list1, "temperature" => 100.0) + # The convection (surface heat transfer) boundary condition is applied + # along the edges 2,3,4. + list2 = selectelem(fens, bfes; box = [Width Width 0.0 Height], inflate = tolerance) + list3 = selectelem(fens, bfes; box = [0.0 Width Height Height], inflate = tolerance) + # The boundary integrals are evaluated using a surface FEMM. cfemm = FEMMHeatDiffSurf( - IntegDomain(subset(bfes, vcat(l2, l3)), GaussRule(1, 3), Thickness), + IntegDomain(subset(bfes, vcat(list2, list3)), GaussRule(1, 3), Thickness), h, ) convection1 = FDataDict("femm" => cfemm, "ambient_temperature" => 0.0) - - # The interior + # The interior integrals are evaluated using a volume FEMM. femm = FEMMHeatDiff(IntegDomain(fes, TriRule(3), Thickness), m) region1 = FDataDict("femm" => femm) - # Make the model data modeldata = FDataDict( "fens" => fens, @@ -99,86 +90,43 @@ function T4NAFEMS_T3_algo() "essential_bcs" => [essential1], "convection_bcs" => [convection1], ) - # Call the solver modeldata = AlgoHeatDiffModule.steadystate(modeldata) - - println("Total time elapsed = ", time() - t0, "s") - - l4 = selectnode(fens; box = [Width Width HeightA HeightA], inflate = tolerance) - - geom = modeldata["geom"] + # Locate the node at the point A [coordinates (Width,HeightA)]. + list4 = selectnode(fens; box=[Width Width HeightA HeightA], inflate=tolerance) + # Collect the temperature at the point A. Temp = modeldata["temp"] - - ## - # Collect the temperature at the point A [coordinates - # (Width,HeightA)]. - println("$(Temp.values[l4][1])") - push!(resultsTempA, Temp.values[l4][1]) + println("$(Temp.values[list4][1])") + push!(resultsTempA, Temp.values[list4][1]) + push!(params, 1.0 / 2^nref) end - ## # These are the computed results for the temperature at point A: println("$( resultsTempA )") + # Richardson extrapolation can be used to estimate the limit. + solnestim, beta, c, residual = + richextrapol(resultsTempA[(end-2):end], params[(end-2):end]) + println("Solution estimate = $(solnestim)") + println("Convergence rate estimate = $(beta)") # Postprocessing geom = modeldata["geom"] Temp = modeldata["temp"] regions = modeldata["regions"] vtkexportmesh( - "T4NAFEMS--T3.vtk", + "T4NAFEMS--T3-solution.vtk", connasarray(regions[1]["femm"].integdomain.fes), - [geom.values Temp.values / 100], + [geom.values (Temp.values / 100)], FinEtools.MeshExportModule.VTK.T3; scalars = [("Temperature", Temp.values)], ) vtkexportmesh( - "T4NAFEMS--T3--base.vtk", + "T4NAFEMS--T3-mesh.vtk", connasarray(regions[1]["femm"].integdomain.fes), - [geom.values 0.0 * Temp.values / 100], + geom.values, FinEtools.MeshExportModule.VTK.T3, ) - - # ## - # # Richardson extrapolation is used to estimate the true solution from the - # # results for the finest three meshes. - # [xestim, beta] = richextrapol(results(end-2:end),mesh_sizes(end-2:end)); - # disp(['Estimated true solution for temperature at A: ' num2str(xestim) ' degrees']) - - # ## - # # Plot the estimated true error. - # figure - # loglog(mesh_sizes,abs(results-xestim)/xestim,'bo-','linewidth',3) - # grid on - # xlabel('log(mesh size)') - # ylabel('log(|estimated temperature error|)') - # set_graphics_defaults - - # ## - # # The estimated true error has a slope of approximately 4 on the log-log - # scale. - # ## - # # Plot the absolute values of the approximate error (differences of - # # successive solutions). - # figure - # loglog(mesh_sizes(2:end),abs(diff(results)),'bo-','linewidth',3) - # Thanksgrid on - # xlabel('log(mesh size)') - # ylabel('log(|approximate temperature error|)') - # set_graphics_defaults - - ## Discussion - # - ## - # The last segment of the approximate error curve is close to the slope of - # the estimated true error. Nevertheless, it would have been more - # reassuring if the three successive approximate errors were located more - # closely on a straight line. - - ## - # The use of uniform mesh-size meshes is sub optimal: it would be more - # efficient to use graded meshes. The tutorial pub_T4NAFEMS_conv_graded - # addresses use of graded meshes in convergence studies. + true end # T4NAFEMS_T3_algo function T4NAFEMS_T6_algo() @@ -228,8 +176,8 @@ function T4NAFEMS_T6_algo() m = MatHeatDiff(kappa) modeldata = nothing - resultsTempA = FFlt[] - params = FFlt[] + resultsTempA = Float64[] + params = Float64[] for nref = 3:7 t0 = time() @@ -251,18 +199,18 @@ function T4NAFEMS_T6_algo() # The prescribed temperature is applied along edge 1 (the bottom # edge in Figure 1).. - l1 = selectnode(fens; box = [0.0 Width 0.0 0.0], inflate = tolerance) - essential1 = FDataDict("node_list" => l1, "temperature" => 100.0) + list1 = selectnode(fens; box = [0.0 Width 0.0 0.0], inflate = tolerance) + essential1 = FDataDict("node_list" => list1, "temperature" => 100.0) ## # The convection boundary condition is applied along the edges # 2,3,4. The elements along the boundary are quadratic line - # elements L3. The order-four Gauss quadrature is sufficiently + # elements list3. The order-four Gauss quadrature is sufficiently # accurate. - l2 = selectelem(fens, bfes; box = [Width Width 0.0 Height], inflate = tolerance) - l3 = selectelem(fens, bfes; box = [0.0 Width Height Height], inflate = tolerance) + list2 = selectelem(fens, bfes; box = [Width Width 0.0 Height], inflate = tolerance) + list3 = selectelem(fens, bfes; box = [0.0 Width Height Height], inflate = tolerance) cfemm = FEMMHeatDiffSurf( - IntegDomain(subset(bfes, vcat(l2, l3)), GaussRule(1, 3), Thickness), + IntegDomain(subset(bfes, vcat(list2, list3)), GaussRule(1, 3), Thickness), h, ) convection1 = FDataDict("femm" => cfemm, "ambient_temperature" => 0.0) @@ -284,7 +232,7 @@ function T4NAFEMS_T6_algo() println("Total time elapsed = ", time() - t0, "s") - l4 = selectnode(fens; box = [Width Width HeightA HeightA], inflate = tolerance) + list4 = selectnode(fens; box = [Width Width HeightA HeightA], inflate = tolerance) geom = modeldata["geom"] Temp = modeldata["temp"] @@ -292,7 +240,7 @@ function T4NAFEMS_T6_algo() ## # Collect the temperature at the point A [coordinates # (Width,HeightA)]. - push!(resultsTempA, Temp.values[l4][1]) + push!(resultsTempA, Temp.values[list4][1]) push!(params, 1.0 / 2^nref) end