From 4e4ab74c1a85df209175664f7cd9c6b0e13c6f1b Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 22 Mar 2026 14:25:43 +0100 Subject: [PATCH 01/71] Rename SolidBodyThermalSurfaceConvection to SolidBodyThermalSurfaceHeatTransfer --- src/felupe/thermal/__init__.py | 4 +-- ...olidbody_thermal_surface_heat_transfer.py} | 35 +++++++++---------- 2 files changed, 19 insertions(+), 20 deletions(-) rename src/felupe/thermal/{_solidbody_thermal_convection.py => _solidbody_thermal_surface_heat_transfer.py} (80%) diff --git a/src/felupe/thermal/__init__.py b/src/felupe/thermal/__init__.py index 133a91f5e..4547c8d12 100644 --- a/src/felupe/thermal/__init__.py +++ b/src/felupe/thermal/__init__.py @@ -1,11 +1,11 @@ from ._solidbody_thermal import SolidBodyThermal -from ._solidbody_thermal_convection import SolidBodyThermalConvection +from ._solidbody_thermal_surface_heat_transfer import SolidBodyThermalSurfaceHeatTransfer from ._solidbody_thermal_heat_flux import SolidBodyThermalHeatFlux from ._time_step import TimeStep __all__ = [ "SolidBodyThermal", - "SolidBodyThermalConvection", + "SolidBodyThermalSurfaceHeatTransfer", "SolidBodyThermalHeatFlux", "TimeStep", ] diff --git a/src/felupe/thermal/_solidbody_thermal_convection.py b/src/felupe/thermal/_solidbody_thermal_surface_heat_transfer.py similarity index 80% rename from src/felupe/thermal/_solidbody_thermal_convection.py rename to src/felupe/thermal/_solidbody_thermal_surface_heat_transfer.py index 4a62e18a8..acfa1a464 100644 --- a/src/felupe/thermal/_solidbody_thermal_convection.py +++ b/src/felupe/thermal/_solidbody_thermal_surface_heat_transfer.py @@ -21,8 +21,8 @@ from ..mechanics import Assemble, Results -class SolidBodyThermalConvection: - r"""A thermal convection boundary condition for a thermal solid body. +class SolidBodyThermalSurfaceHeatTransfer: + r"""A thermal boundary condition for a thermal solid body. Parameters ---------- @@ -31,7 +31,7 @@ class SolidBodyThermalConvection: coefficient : float The convection coefficient :math:`h` in W/(m^2*K). temperature : float - The ambient temperature :math:`T_\infty` in K. + The ambient temperature :math:`T_\infty` in °C or K. Notes ----- @@ -52,36 +52,35 @@ class SolidBodyThermalConvection: >>> temperature = fem.Field(region, dim=1) >>> field = fem.FieldContainer([temperature]) >>> - >>> region_convection = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) - >>> temperature_convection = fem.Field(region_convection, dim=1) - >>> field_convection = fem.FieldContainer([temperature_convection]) + >>> region_heat_transfer = fem.RegionQuadBoundary(mesh, mask=mesh.x == 1.0) + >>> temperature_heat_transfer = fem.Field(region_heat_transfer, dim=1) + >>> field_heat_transfer = fem.FieldContainer([temperature_heat_transfer]) >>> >>> boundaries = fem.BoundaryDict( ... left=fem.Boundary(temperature, fx=0), - ... right=fem.Boundary(temperature, fx=1), ... ) >>> >>> solid = fem.thermal.SolidBodyThermal( ... field=field, - ... mass_density=1.0, # kg/m^3 - ... specific_heat_capacity=1.0, # J/(kg*K) - ... time_step=0.01, # s + ... mass_density=1400.0, # kg/m^3 + ... specific_heat_capacity=1000.0, # J/(kg*K) + ... time_step=720.0, # s ... thermal_conductivity=1.0, # W/(m*K) ... ) - >>> convection = fem.thermal.SolidBodyThermalConvection( - ... field=field_convection, - ... coefficient=1.0, # W/(m^2*K) - ... temperature=10.0, # K + >>> heat_transfer = fem.thermal.SolidBodyThermalSurfaceHeatTransfer( + ... field=field_heat_transfer, + ... coefficient=7.69, # W/(m^2 K) + ... temperature=10.0, # °C ... ) >>> time = fem.thermal.TimeStep([solid]) >>> table = fem.math.linsteps([0, 1], num=10) >>> ramp = { - ... boundaries["right"]: 10 * table, - ... time: 0.1 * table, - ... convection: 100 * table, + ... boundaries["left"]: 10 * table, # surface temperature + ... time: 18000 * table, # five hours + ... convection: 40 * table, # air temperature w/ transfer coeff. ... } >>> step = fem.Step( - ... items=[time, solid, convection], ramp=ramp, boundaries=boundaries + ... items=[time, solid, heat_transfer], ramp=ramp, boundaries=boundaries ... ) >>> job = fem.Job(steps=[step]).evaluate( ... filename="result.xdmf", # result file for Paraview From 1c0dfa60496203a6f01e417cd350c10fc9f3ac4d Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 22 Mar 2026 14:33:08 +0100 Subject: [PATCH 02/71] Add SolidBodyThermal example suggestion. WIP. --- examples/ex22_solid_body_thermal.py | 176 ++++++++++++++++++++++++++++ 1 file changed, 176 insertions(+) create mode 100644 examples/ex22_solid_body_thermal.py diff --git a/examples/ex22_solid_body_thermal.py b/examples/ex22_solid_body_thermal.py new file mode 100644 index 000000000..e077198ad --- /dev/null +++ b/examples/ex22_solid_body_thermal.py @@ -0,0 +1,176 @@ +r""" +Thermal Analysis +------------------------ + +.. topic:: Thermal analysis of simple construction. + + * use SolidBodyThermal + + * view the temperature field + + +This example describes a simple two-dimensional light-weight construction +system set up with nine SolidBodyThermal solids. +""" + +# sphinx_gallery_thumbnail_number = -1 + +import matplotlib.pyplot as plt +import numpy as np + +import felupe as fem + +mesh = fem.Mesh( + points=[ + [2.5, 0, 0], + [-1.25, 1.25, 0], + [1, 2, 0], + [-0.5, 1.5, 1.5], + [-2.5, 4.5, 2.5], + ], + cells=[[0, 3], [1, 3], [2, 3], [2, 4], [1, 4], [3, 4]], + cell_type="line", +) +region = fem.RegionTruss(mesh) +field = fem.FieldContainer([fem.Field(region, dim=3)]) + +# %% +# Beside points and cells we have to define displacement boundary conditions, external +# forces and the constitutive material formulation for the trusses. +boundaries = fem.BoundaryDict( + fixed_xyz=fem.Boundary(field[0], mask=[1, 1, 1, 0, 0]), + fixed_y=fem.Boundary(field[0], mask=[0, 0, 0, 0, 1], skip=(1, 0, 1)), +) +dof0, dof1 = fem.dof.partition(field, boundaries) + +solid = fem.TrussBody( + umat=fem.LinearElastic1D(E=1), + field=field, + area=[0.75, 1, 0.5, 0.75, 1, 1], +) + +force_3 = np.array([1, 1, -1]) +force_4 = np.array([-2, 0, -2]) + +load_3 = fem.PointLoad(field, [3], force_3) +load_4 = fem.PointLoad(field, [4], force_4) + +# %% +# The undeformed configuration is plotted in a 3d-view. +plotter = mesh.plot( + line_width=10, + render_lines_as_tubes=True, + show_edges=False, +) +plotter.add_points( + mesh.points, + color="black", + point_size=20, + render_points_as_spheres=True, +) +plotter = boundaries.plot(plotter=plotter) +plotter = load_3.plot(plotter=plotter, color="green", deformed=False) +plotter = load_4.plot(plotter=plotter, color="green", deformed=False) + +plotter.show() + + +# %% +# For the numeric continuation, the equilibrium function ``fun`` and its derivatives +# w.r.t. the displacement field ``dfun_du`` and the load-proportionality-factor +# ``dfun_dlpf`` have to be defined. Here, we're only interested in the active degrees of +# freedom. +def fun(x, lpf, *args): + field[0].values.ravel()[dof1] = x + load_3.update(force_3 * lpf) + load_4.update(force_4 * lpf) + return fem.tools.fun([solid, load_3, load_4], field)[dof1] + + +def dfun_du(x, lpf, *args): + field[0].values.ravel()[dof1] = x + K = fem.tools.jac([solid, load_3, load_4], field) + return fem.solve.partition(field, K, dof1, dof0)[2] + + +def dfun_dlpf(x, lpf, *args): + load_3.update(force_3) + load_4.update(force_4) + return fem.tools.fun([load_3, load_4], field)[dof1] + + +# %% +# Now that the model is finished, some additional settings have to be chosen. Initial +# allowed incremental system vector components for both the displacement vector and the +# load-proportionality-factor (LPF) have to be specified. We use ``dlpf = 0.005`` and +# ``du = 0.05`` (figured out after some trial and error). Both parameters can't be +# specified automatically, as they depend on the model configuration. The job will be +# limited to a total amount of 163 increments (again, the total number has been figured +# out after some job runs to get good looking plots). +res = contique.solve( + fun=fun, + jac=[dfun_du, dfun_dlpf], + x0=field[0][dof1], + lpf0=0, + dxmax=0.05, + dlpfmax=0.005, + maxsteps=163, + rebalance=True, + overshoot=1.25, + tol=1e-8, + low=1e-2, + high=4, + maxiter=8, +) +X = np.array([r.x for r in res]) + +# %% +# To visualize the deformed state of the model for increment 40 the deformed model plot +# is generated. +field[0].values.ravel()[dof1] = X[40, :-1] +force = solid.evaluate.gradient(field) * solid.area +plotter = field.view(cell_data={"Force": force}).plot( + "Force", + line_width=10, + show_undeformed=False, + view="xy", + cmap="coolwarm", + clim=[-abs(force).max(), abs(force).max()], + render_lines_as_tubes=True, + show_edges=False, +) +plotter.add_points( + mesh.points + field[0].values, + color="black", + point_size=20, + render_points_as_spheres=True, +) +plotter.show() + +# %% +# Path-tracing of the displacement-LPF curves +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# The path-tracing of the deformation process is shown as a History Plot of +# Displacement-LPF curves for all active DOF. Strong geometrical nonlinearities are +# observed for all active DOF. +fig, ax = plt.subplots() +ax.plot(*X[:, [0, -1]].T, ".-", label="Point 3") +ax.plot(*X[:, [3, -1]].T, ".-", label="Point 4") +ax.set_xlabel("Displacement X") +ax.set_ylabel("LPF") +ax.legend() + +# %% +fig, ax = plt.subplots() +ax.plot(*X[:, [1, -1]].T, ".-", label="Point 3") +ax.set_xlabel("Displacement Y") +ax.set_ylabel("LPF") +ax.legend() + +# %% +fig, ax = plt.subplots() +ax.plot(*X[:, [2, -1]].T, ".-", label="Point 3") +ax.plot(*X[:, [4, -1]].T, ".-", label="Point 4") +ax.set_xlabel("Displacement Z") +ax.set_ylabel("LPF") +ax.legend() From fe1fe6f1d4ad8f6187b78f80c1034ad61f8462e7 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 22 Mar 2026 15:19:01 +0100 Subject: [PATCH 03/71] Add actual functionality to SolidBodyThermal example. Still WIP. --- examples/ex22_solid_body_thermal.py | 263 +++++++++++++--------------- 1 file changed, 123 insertions(+), 140 deletions(-) diff --git a/examples/ex22_solid_body_thermal.py b/examples/ex22_solid_body_thermal.py index e077198ad..4c3d6d17e 100644 --- a/examples/ex22_solid_body_thermal.py +++ b/examples/ex22_solid_body_thermal.py @@ -20,157 +20,140 @@ import felupe as fem -mesh = fem.Mesh( - points=[ - [2.5, 0, 0], - [-1.25, 1.25, 0], - [1, 2, 0], - [-0.5, 1.5, 1.5], - [-2.5, 4.5, 2.5], - ], - cells=[[0, 3], [1, 3], [2, 3], [2, 4], [1, 4], [3, 4]], - cell_type="line", -) -region = fem.RegionTruss(mesh) -field = fem.FieldContainer([fem.Field(region, dim=3)]) - -# %% -# Beside points and cells we have to define displacement boundary conditions, external -# forces and the constitutive material formulation for the trusses. -boundaries = fem.BoundaryDict( - fixed_xyz=fem.Boundary(field[0], mask=[1, 1, 1, 0, 0]), - fixed_y=fem.Boundary(field[0], mask=[0, 0, 0, 0, 1], skip=(1, 0, 1)), -) -dof0, dof1 = fem.dof.partition(field, boundaries) - -solid = fem.TrussBody( - umat=fem.LinearElastic1D(E=1), - field=field, - area=[0.75, 1, 0.5, 0.75, 1, 1], -) - -force_3 = np.array([1, 1, -1]) -force_4 = np.array([-2, 0, -2]) - -load_3 = fem.PointLoad(field, [3], force_3) -load_4 = fem.PointLoad(field, [4], force_4) - # %% -# The undeformed configuration is plotted in a 3d-view. -plotter = mesh.plot( - line_width=10, - render_lines_as_tubes=True, - show_edges=False, -) -plotter.add_points( - mesh.points, - color="black", - point_size=20, - render_points_as_spheres=True, -) -plotter = boundaries.plot(plotter=plotter) -plotter = load_3.plot(plotter=plotter, color="green", deformed=False) -plotter = load_4.plot(plotter=plotter, color="green", deformed=False) - -plotter.show() +# Define material properties as lists [plasterboard, insulation, wood]. +density=[700, 20, 500] # kg/m^3 +specific_heat=[1125, 1450, 1000] # J/(kg K) +thermal_conductivity=[0.4, 0.032, 0.13] # W/(m K) +material_index = [0, 0, 0, 1, 2, 1, 0, 0, 0] # %% -# For the numeric continuation, the equilibrium function ``fun`` and its derivatives -# w.r.t. the displacement field ``dfun_du`` and the load-proportionality-factor -# ``dfun_dlpf`` have to be defined. Here, we're only interested in the active degrees of -# freedom. -def fun(x, lpf, *args): - field[0].values.ravel()[dof1] = x - load_3.update(force_3 * lpf) - load_4.update(force_4 * lpf) - return fem.tools.fun([solid, load_3, load_4], field)[dof1] - - -def dfun_du(x, lpf, *args): - field[0].values.ravel()[dof1] = x - K = fem.tools.jac([solid, load_3, load_4], field) - return fem.solve.partition(field, K, dof1, dof0)[2] - - -def dfun_dlpf(x, lpf, *args): - load_3.update(force_3) - load_4.update(force_4) - return fem.tools.fun([load_3, load_4], field)[dof1] - +# Set up one mesh per material area, split horizontally and vertically leading +# to nine mesh areas total. +mesh_list = [fem.mesh.Grid(np.linspace(0, 0.018, 6), # plasterboard + np.linspace(0, 0.47, 18)), + fem.mesh.Grid(np.linspace(0, 0.018, 6), + np.linspace(0.47, 0.53, 8)), + fem.mesh.Grid(np.linspace(0, 0.018, 6), + np.linspace(0.53, 1.0, 18)), + fem.mesh.Grid(np.linspace(0.018, 0.268, 12), # insulation + np.linspace(0, 0.47, 18)), + fem.mesh.Grid(np.linspace(0.018, 0.268, 12), # wood + np.linspace(0.47, 0.53, 8)), + fem.mesh.Grid(np.linspace(0.018, 0.268, 12), # insulation + np.linspace(0.53, 1.0, 18)), + fem.mesh.Grid(np.linspace(0.268, 0.286, 6), # plasterboard + np.linspace(0, 0.47, 18)), + fem.mesh.Grid(np.linspace(0.268, 0.286, 6), + np.linspace(0.47, 0.53, 8)), + fem.mesh.Grid(np.linspace(0.268, 0.286, 6), + np.linspace(0.53, 1.0, 18)), +] # %% -# Now that the model is finished, some additional settings have to be chosen. Initial -# allowed incremental system vector components for both the displacement vector and the -# load-proportionality-factor (LPF) have to be specified. We use ``dlpf = 0.005`` and -# ``du = 0.05`` (figured out after some trial and error). Both parameters can't be -# specified automatically, as they depend on the model configuration. The job will be -# limited to a total amount of 163 increments (again, the total number has been figured -# out after some job runs to get good looking plots). -res = contique.solve( - fun=fun, - jac=[dfun_du, dfun_dlpf], - x0=field[0][dof1], - lpf0=0, - dxmax=0.05, - dlpfmax=0.005, - maxsteps=163, - rebalance=True, - overshoot=1.25, - tol=1e-8, - low=1e-2, - high=4, - maxiter=8, -) -X = np.array([r.x for r in res]) +# Beside points and cells we have to define temperature boundary conditions, +# and the materials for the solid bodies. +mesh_container = fp.MeshContainer(mesh_list, merge=True) +regions = [fp.RegionQuad(m) for m in mesh_container] +field_list = [fp.Field(r, dim=1).as_container() for r in regions] + +# top level +mesh = mesh_container.stack() +region = fp.RegionQuad(mesh) +field = fp.Field(region, dim=1).as_container() + +temperature = field[0] # define top-level field values as temperature + +external_loc = fp.RegionQuadBoundary(mesh, mask=mesh.x == x.min()) +external_temp = fp.Field(external_loc, dim=1) +external_fld = fp.FieldContainer([external_temp]) + +internal_loc = fp.RegionQuadBoundary(mesh, mask=mesh.x == x.max()) +internal_temp = fp.Field(internal_loc, dim=1) +internal_fld = fp.FieldContainer([internal_temp]) + +boundaries = { # Heat transfer coefficients. + "external" : SolidBodySurfaceHeatTransfer( + field=external_fld, + coefficient=25.0, # W/(m^2 K) + temperature=20.0, # °C + ), + "internal" : SolidBodySurfaceHeatTransfer( + field=internal_fld, + coefficient=7.69, # W/(m^2 K) + temperature=20.0, # °C + ) +} + +materials = [] +for imat, fld in enumerate(field_list): + cur_mat = material_index[imat] + materials.append( + SolidBodyThermal( + fld, + # field_list[idx], + density[cur_mat], + specific_heat[cur_mat], + thermal_conductivity[cur_mat]) + ) # %% -# To visualize the deformed state of the model for increment 40 the deformed model plot -# is generated. -field[0].values.ravel()[dof1] = X[40, :-1] -force = solid.evaluate.gradient(field) * solid.area -plotter = field.view(cell_data={"Force": force}).plot( - "Force", - line_width=10, - show_undeformed=False, - view="xy", - cmap="coolwarm", - clim=[-abs(force).max(), abs(force).max()], - render_lines_as_tubes=True, - show_edges=False, -) -plotter.add_points( - mesh.points + field[0].values, - color="black", - point_size=20, - render_points_as_spheres=True, -) -plotter.show() +# Define helper and callback functions for surface flux value retrieval. +def boundary_flux(bmesh, material, coord_index=0, coord_value=0.): + """ + Extract boundary flux from mesh/material. Code by adtzlr.""" + points = bmesh.points[bmesh.cells] # shape(n_cells, 4, 2) + mask = np.isclose(points[..., coord_index], coord_value).T # shape(4, n_cells) + flux = material.evaluate.stress()[0][coord_index, mask] # shape(38) + return flux + +def callback(stepnumber, substepnumber, substep, flux_data): + """ + Extract stress data (flux), shape(n_cells, 4) per substep.""" + flux_data[0].append( + boundary_flux(mesh_list[0], materials[0], + coord_value=field.region.mesh.x.min()) + ) + + flux_data[1].append( + boundary_flux(mesh_list[-1], materials[-1], + coord_value=field.region.mesh.x.max()) + ) + +time_steps = fem.math.linsteps([0, 24*3600], num=int(24*3600/720)) + +t_int = np.around(20 + 1 * np.sin(2*np.pi * time_steps / 86400), 5) +t_ext = np.around(0 + 1 * np.sin(2*np.pi * time_steps / 86400), 5) # %% -# Path-tracing of the displacement-LPF curves -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# The path-tracing of the deformation process is shown as a History Plot of -# Displacement-LPF curves for all active DOF. Strong geometrical nonlinearities are -# observed for all active DOF. -fig, ax = plt.subplots() -ax.plot(*X[:, [0, -1]].T, ".-", label="Point 3") -ax.plot(*X[:, [3, -1]].T, ".-", label="Point 4") -ax.set_xlabel("Displacement X") -ax.set_ylabel("LPF") -ax.legend() +# Set up boundary values, job description and solve. +time = TimeStep(materials) + +ramp = {boundaries["internal"]: t_int, + boundaries["external"]: t_ext, + time: time_steps} + +step = fp.Step(items=[time] + materials, + ramp=ramp, + boundaries=boundaries) + +flux_data = {0: [], 1: []} + +job = fp.Job( + steps=[step], + callback=callback, + flux_data=flux_data).evaluate( + x0=field, + filename="result.xdmf", # result file for Paraview + point_data={"Temperature": lambda field, substep: temperature.values}, + point_data_default=False, + cell_data_default=False, + ) # %% -fig, ax = plt.subplots() -ax.plot(*X[:, [1, -1]].T, ".-", label="Point 3") -ax.set_xlabel("Displacement Y") -ax.set_ylabel("LPF") -ax.legend() +# Internal and external surface heat flux vs. time. # %% -fig, ax = plt.subplots() -ax.plot(*X[:, [2, -1]].T, ".-", label="Point 3") -ax.plot(*X[:, [4, -1]].T, ".-", label="Point 4") -ax.set_xlabel("Displacement Z") -ax.set_ylabel("LPF") -ax.legend() +# Temperature field at end of simulation period. + From 53fb632d19591b22d3c1137a4d4fa8f906499830 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 22 Mar 2026 17:57:28 +0100 Subject: [PATCH 04/71] Shorten 'HeatTransfer' file and class names. Debug example, runs through. Needs plots added next. --- examples/ex22_solid_body_thermal.py | 52 +++++++++---------- src/felupe/thermal/__init__.py | 4 +- ...py => _solidbody_surface_heat_transfer.py} | 12 +++-- 3 files changed, 35 insertions(+), 33 deletions(-) rename src/felupe/thermal/{_solidbody_thermal_surface_heat_transfer.py => _solidbody_surface_heat_transfer.py} (92%) diff --git a/examples/ex22_solid_body_thermal.py b/examples/ex22_solid_body_thermal.py index 4c3d6d17e..03b24337d 100644 --- a/examples/ex22_solid_body_thermal.py +++ b/examples/ex22_solid_body_thermal.py @@ -4,7 +4,7 @@ .. topic:: Thermal analysis of simple construction. - * use SolidBodyThermal + * use SolidBodyThermal and SolidBodySurfaceHeatTransfer * view the temperature field @@ -54,45 +54,42 @@ # %% # Beside points and cells we have to define temperature boundary conditions, # and the materials for the solid bodies. -mesh_container = fp.MeshContainer(mesh_list, merge=True) -regions = [fp.RegionQuad(m) for m in mesh_container] -field_list = [fp.Field(r, dim=1).as_container() for r in regions] +mesh_container = fem.MeshContainer(mesh_list, merge=True) +regions = [fem.RegionQuad(m) for m in mesh_container] +field_list = [fem.Field(r, dim=1).as_container() for r in regions] # top level mesh = mesh_container.stack() -region = fp.RegionQuad(mesh) -field = fp.Field(region, dim=1).as_container() +region = fem.RegionQuad(mesh) +field = fem.Field(region, dim=1).as_container() temperature = field[0] # define top-level field values as temperature -external_loc = fp.RegionQuadBoundary(mesh, mask=mesh.x == x.min()) -external_temp = fp.Field(external_loc, dim=1) -external_fld = fp.FieldContainer([external_temp]) +external_loc = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.min()) +external_temp = fem.Field(external_loc, dim=1) +external_fld = fem.FieldContainer([external_temp]) -internal_loc = fp.RegionQuadBoundary(mesh, mask=mesh.x == x.max()) -internal_temp = fp.Field(internal_loc, dim=1) -internal_fld = fp.FieldContainer([internal_temp]) +internal_loc = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.max()) +internal_temp = fem.Field(internal_loc, dim=1) +internal_fld = fem.FieldContainer([internal_temp]) -boundaries = { # Heat transfer coefficients. - "external" : SolidBodySurfaceHeatTransfer( +external_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( field=external_fld, coefficient=25.0, # W/(m^2 K) temperature=20.0, # °C - ), - "internal" : SolidBodySurfaceHeatTransfer( + ) +internal_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( field=internal_fld, coefficient=7.69, # W/(m^2 K) temperature=20.0, # °C ) -} materials = [] for imat, fld in enumerate(field_list): cur_mat = material_index[imat] materials.append( - SolidBodyThermal( + fem.thermal.SolidBodyThermal( fld, - # field_list[idx], density[cur_mat], specific_heat[cur_mat], thermal_conductivity[cur_mat]) @@ -128,19 +125,18 @@ def callback(stepnumber, substepnumber, substep, flux_data): # %% # Set up boundary values, job description and solve. -time = TimeStep(materials) +time = fem.thermal.TimeStep(materials) -ramp = {boundaries["internal"]: t_int, - boundaries["external"]: t_ext, +ramp = {internal_heat_transfer: t_int, + external_heat_transfer: t_ext, time: time_steps} -step = fp.Step(items=[time] + materials, - ramp=ramp, - boundaries=boundaries) +step = fem.Step(items=[time] + materials + [internal_heat_transfer, external_heat_transfer], + ramp=ramp) flux_data = {0: [], 1: []} -job = fp.Job( +job = fem.Job( steps=[step], callback=callback, flux_data=flux_data).evaluate( @@ -153,6 +149,10 @@ def callback(stepnumber, substepnumber, substep, flux_data): # %% # Internal and external surface heat flux vs. time. +q_ext = np.mean(flux_data[0], axis=1) # W/m2 => len=1 +q_int = np.mean(flux_data[1], axis=1) # W/m2 + + # %% # Temperature field at end of simulation period. diff --git a/src/felupe/thermal/__init__.py b/src/felupe/thermal/__init__.py index 4547c8d12..3ab10e0a5 100644 --- a/src/felupe/thermal/__init__.py +++ b/src/felupe/thermal/__init__.py @@ -1,11 +1,11 @@ from ._solidbody_thermal import SolidBodyThermal -from ._solidbody_thermal_surface_heat_transfer import SolidBodyThermalSurfaceHeatTransfer +from ._solidbody_surface_heat_transfer import SolidBodySurfaceHeatTransfer from ._solidbody_thermal_heat_flux import SolidBodyThermalHeatFlux from ._time_step import TimeStep __all__ = [ "SolidBodyThermal", - "SolidBodyThermalSurfaceHeatTransfer", + "SolidBodySurfaceHeatTransfer", "SolidBodyThermalHeatFlux", "TimeStep", ] diff --git a/src/felupe/thermal/_solidbody_thermal_surface_heat_transfer.py b/src/felupe/thermal/_solidbody_surface_heat_transfer.py similarity index 92% rename from src/felupe/thermal/_solidbody_thermal_surface_heat_transfer.py rename to src/felupe/thermal/_solidbody_surface_heat_transfer.py index acfa1a464..d46c941a4 100644 --- a/src/felupe/thermal/_solidbody_thermal_surface_heat_transfer.py +++ b/src/felupe/thermal/_solidbody_surface_heat_transfer.py @@ -21,8 +21,8 @@ from ..mechanics import Assemble, Results -class SolidBodyThermalSurfaceHeatTransfer: - r"""A thermal boundary condition for a thermal solid body. +class SolidBodySurfaceHeatTransfer: + r"""A surface boundary condition for a thermal solid body. Parameters ---------- @@ -67,7 +67,7 @@ class SolidBodyThermalSurfaceHeatTransfer: ... time_step=720.0, # s ... thermal_conductivity=1.0, # W/(m*K) ... ) - >>> heat_transfer = fem.thermal.SolidBodyThermalSurfaceHeatTransfer( + >>> heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( ... field=field_heat_transfer, ... coefficient=7.69, # W/(m^2 K) ... temperature=10.0, # °C @@ -77,7 +77,7 @@ class SolidBodyThermalSurfaceHeatTransfer: >>> ramp = { ... boundaries["left"]: 10 * table, # surface temperature ... time: 18000 * table, # five hours - ... convection: 40 * table, # air temperature w/ transfer coeff. + ... heat_transfer: 40 * table, # air temperature w/ transfer coeff. ... } >>> step = fem.Step( ... items=[time, solid, heat_transfer], ramp=ramp, boundaries=boundaries @@ -114,8 +114,10 @@ def __init__(self, field, coefficient, temperature): vector=self._vector, matrix=self._matrix, multiplier=-1.0 ) - def update(self, temperature): + def update(self, temperature, coefficient=None): self.results.temperature = temperature + if coefficient is not None: + self.results.coefficient = coefficient def _vector(self, field=None, **kwargs): if field is not None: From 2a5cb260b0f7b445831c9183bde1631a59dad485 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 22 Mar 2026 18:27:50 +0100 Subject: [PATCH 05/71] Add results views (not very refined). --- examples/ex22_solid_body_thermal.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/examples/ex22_solid_body_thermal.py b/examples/ex22_solid_body_thermal.py index 03b24337d..a7f01252d 100644 --- a/examples/ex22_solid_body_thermal.py +++ b/examples/ex22_solid_body_thermal.py @@ -152,8 +152,14 @@ def callback(stepnumber, substepnumber, substep, flux_data): q_ext = np.mean(flux_data[0], axis=1) # W/m2 => len=1 q_int = np.mean(flux_data[1], axis=1) # W/m2 - +fig, ax = plt.subplots() +ax.plot(time_steps.T/3600, q_int.T) +ax.plot(time_steps/3600, q_ext) +ax.set(xlabel='time (h)', ylabel='surface heat flux (W/(m^2 K))') +# fig.savefig("test.png") +plt.show() # %% # Temperature field at end of simulation period. - +view = mesh.view(point_data={"Field": temperature.values}) +view.plot("Field").show() From 4dcfa97aa715f9d6b7c8a88741e42f4027ff3765 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 22 Mar 2026 18:33:22 +0100 Subject: [PATCH 06/71] Update introduction text. Update material data values. --- examples/ex22_solid_body_thermal.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/examples/ex22_solid_body_thermal.py b/examples/ex22_solid_body_thermal.py index a7f01252d..7e77da6aa 100644 --- a/examples/ex22_solid_body_thermal.py +++ b/examples/ex22_solid_body_thermal.py @@ -10,7 +10,9 @@ This example describes a simple two-dimensional light-weight construction -system set up with nine SolidBodyThermal solids. +system set up with nine SolidBodyThermal solids. The temperature boundary +conditions have a +-1 K sinusoidal variation around their average value +with a period of 24 h. """ # sphinx_gallery_thumbnail_number = -1 @@ -22,9 +24,9 @@ # %% # Define material properties as lists [plasterboard, insulation, wood]. -density=[700, 20, 500] # kg/m^3 -specific_heat=[1125, 1450, 1000] # J/(kg K) -thermal_conductivity=[0.4, 0.032, 0.13] # W/(m K) +density=[1000, 20, 500] # kg/m^3 +specific_heat=[1125, 1450, 1600] # J/(kg K) +thermal_conductivity=[0.4, 0.035, 0.16] # W/(m K) material_index = [0, 0, 0, 1, 2, 1, 0, 0, 0] From fe498dd38ed02babfa0202d7d6655c33531e38ef Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 22 Mar 2026 18:39:08 +0100 Subject: [PATCH 07/71] Adapt class notes text. --- src/felupe/thermal/_solidbody_surface_heat_transfer.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_heat_transfer.py b/src/felupe/thermal/_solidbody_surface_heat_transfer.py index d46c941a4..c0e55f636 100644 --- a/src/felupe/thermal/_solidbody_surface_heat_transfer.py +++ b/src/felupe/thermal/_solidbody_surface_heat_transfer.py @@ -35,11 +35,11 @@ class SolidBodySurfaceHeatTransfer: Notes ----- - This class represents a thermal convection boundary condition for a thermal solid - body, which is used to model heat convection at the boundary of a solid material. - The convection coefficient is used to calculate the heat flux at the boundary based - on the difference between the temperature at the boundary and the ambient - temperature. + This class represents a boundary condition for a thermal solid body, which + is used to model heat transfer (convection, radiation) at the boundary of a + solid material. The coefficient is used to calculate the heat flux at the + boundary based on the difference between the temperature at the boundary + and the ambient temperature. Examples -------- From dc18cbfaeb20e976295e62b031322782d5848170 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 23 Mar 2026 07:23:38 +0100 Subject: [PATCH 08/71] format black, run isort --- examples/ex22_solid_body_thermal.py | 114 ++++++++++++++-------------- src/felupe/thermal/__init__.py | 2 +- 2 files changed, 56 insertions(+), 60 deletions(-) diff --git a/examples/ex22_solid_body_thermal.py b/examples/ex22_solid_body_thermal.py index 7e77da6aa..7c06d213f 100644 --- a/examples/ex22_solid_body_thermal.py +++ b/examples/ex22_solid_body_thermal.py @@ -24,33 +24,31 @@ # %% # Define material properties as lists [plasterboard, insulation, wood]. -density=[1000, 20, 500] # kg/m^3 -specific_heat=[1125, 1450, 1600] # J/(kg K) -thermal_conductivity=[0.4, 0.035, 0.16] # W/(m K) +density = [1000, 20, 500] # kg/m^3 +specific_heat = [1125, 1450, 1600] # J/(kg K) +thermal_conductivity = [0.4, 0.035, 0.16] # W/(m K) material_index = [0, 0, 0, 1, 2, 1, 0, 0, 0] # %% # Set up one mesh per material area, split horizontally and vertically leading # to nine mesh areas total. -mesh_list = [fem.mesh.Grid(np.linspace(0, 0.018, 6), # plasterboard - np.linspace(0, 0.47, 18)), - fem.mesh.Grid(np.linspace(0, 0.018, 6), - np.linspace(0.47, 0.53, 8)), - fem.mesh.Grid(np.linspace(0, 0.018, 6), - np.linspace(0.53, 1.0, 18)), - fem.mesh.Grid(np.linspace(0.018, 0.268, 12), # insulation - np.linspace(0, 0.47, 18)), - fem.mesh.Grid(np.linspace(0.018, 0.268, 12), # wood - np.linspace(0.47, 0.53, 8)), - fem.mesh.Grid(np.linspace(0.018, 0.268, 12), # insulation - np.linspace(0.53, 1.0, 18)), - fem.mesh.Grid(np.linspace(0.268, 0.286, 6), # plasterboard - np.linspace(0, 0.47, 18)), - fem.mesh.Grid(np.linspace(0.268, 0.286, 6), - np.linspace(0.47, 0.53, 8)), - fem.mesh.Grid(np.linspace(0.268, 0.286, 6), - np.linspace(0.53, 1.0, 18)), +mesh_list = [ + fem.mesh.Grid(np.linspace(0, 0.018, 6), np.linspace(0, 0.47, 18)), # plasterboard + fem.mesh.Grid(np.linspace(0, 0.018, 6), np.linspace(0.47, 0.53, 8)), + fem.mesh.Grid(np.linspace(0, 0.018, 6), np.linspace(0.53, 1.0, 18)), + fem.mesh.Grid( + np.linspace(0.018, 0.268, 12), np.linspace(0, 0.47, 18) # insulation + ), + fem.mesh.Grid(np.linspace(0.018, 0.268, 12), np.linspace(0.47, 0.53, 8)), # wood + fem.mesh.Grid( + np.linspace(0.018, 0.268, 12), np.linspace(0.53, 1.0, 18) # insulation + ), + fem.mesh.Grid( + np.linspace(0.268, 0.286, 6), np.linspace(0, 0.47, 18) # plasterboard + ), + fem.mesh.Grid(np.linspace(0.268, 0.286, 6), np.linspace(0.47, 0.53, 8)), + fem.mesh.Grid(np.linspace(0.268, 0.286, 6), np.linspace(0.53, 1.0, 18)), ] # %% @@ -76,30 +74,29 @@ internal_fld = fem.FieldContainer([internal_temp]) external_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( - field=external_fld, - coefficient=25.0, # W/(m^2 K) - temperature=20.0, # °C - ) + field=external_fld, + coefficient=25.0, # W/(m^2 K) + temperature=20.0, # °C +) internal_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( - field=internal_fld, - coefficient=7.69, # W/(m^2 K) - temperature=20.0, # °C - ) + field=internal_fld, + coefficient=7.69, # W/(m^2 K) + temperature=20.0, # °C +) materials = [] for imat, fld in enumerate(field_list): cur_mat = material_index[imat] materials.append( fem.thermal.SolidBodyThermal( - fld, - density[cur_mat], - specific_heat[cur_mat], - thermal_conductivity[cur_mat]) + fld, density[cur_mat], specific_heat[cur_mat], thermal_conductivity[cur_mat] ) + ) + # %% # Define helper and callback functions for surface flux value retrieval. -def boundary_flux(bmesh, material, coord_index=0, coord_value=0.): +def boundary_flux(bmesh, material, coord_index=0, coord_value=0.0): """ Extract boundary flux from mesh/material. Code by adtzlr.""" points = bmesh.points[bmesh.cells] # shape(n_cells, 4, 2) @@ -107,47 +104,46 @@ def boundary_flux(bmesh, material, coord_index=0, coord_value=0.): flux = material.evaluate.stress()[0][coord_index, mask] # shape(38) return flux + def callback(stepnumber, substepnumber, substep, flux_data): """ Extract stress data (flux), shape(n_cells, 4) per substep.""" flux_data[0].append( - boundary_flux(mesh_list[0], materials[0], - coord_value=field.region.mesh.x.min()) + boundary_flux(mesh_list[0], materials[0], coord_value=field.region.mesh.x.min()) ) flux_data[1].append( - boundary_flux(mesh_list[-1], materials[-1], - coord_value=field.region.mesh.x.max()) + boundary_flux( + mesh_list[-1], materials[-1], coord_value=field.region.mesh.x.max() + ) ) -time_steps = fem.math.linsteps([0, 24*3600], num=int(24*3600/720)) -t_int = np.around(20 + 1 * np.sin(2*np.pi * time_steps / 86400), 5) -t_ext = np.around(0 + 1 * np.sin(2*np.pi * time_steps / 86400), 5) +time_steps = fem.math.linsteps([0, 24 * 3600], num=int(24 * 3600 / 720)) + +t_int = np.around(20 + 1 * np.sin(2 * np.pi * time_steps / 86400), 5) +t_ext = np.around(0 + 1 * np.sin(2 * np.pi * time_steps / 86400), 5) # %% # Set up boundary values, job description and solve. time = fem.thermal.TimeStep(materials) -ramp = {internal_heat_transfer: t_int, - external_heat_transfer: t_ext, - time: time_steps} +ramp = {internal_heat_transfer: t_int, external_heat_transfer: t_ext, time: time_steps} -step = fem.Step(items=[time] + materials + [internal_heat_transfer, external_heat_transfer], - ramp=ramp) +step = fem.Step( + items=[time] + materials + [internal_heat_transfer, external_heat_transfer], + ramp=ramp, +) flux_data = {0: [], 1: []} -job = fem.Job( - steps=[step], - callback=callback, - flux_data=flux_data).evaluate( - x0=field, - filename="result.xdmf", # result file for Paraview - point_data={"Temperature": lambda field, substep: temperature.values}, - point_data_default=False, - cell_data_default=False, - ) +job = fem.Job(steps=[step], callback=callback, flux_data=flux_data).evaluate( + x0=field, + filename="result.xdmf", # result file for Paraview + point_data={"Temperature": lambda field, substep: temperature.values}, + point_data_default=False, + cell_data_default=False, +) # %% # Internal and external surface heat flux vs. time. @@ -155,9 +151,9 @@ def callback(stepnumber, substepnumber, substep, flux_data): q_int = np.mean(flux_data[1], axis=1) # W/m2 fig, ax = plt.subplots() -ax.plot(time_steps.T/3600, q_int.T) -ax.plot(time_steps/3600, q_ext) -ax.set(xlabel='time (h)', ylabel='surface heat flux (W/(m^2 K))') +ax.plot(time_steps.T / 3600, q_int.T) +ax.plot(time_steps / 3600, q_ext) +ax.set(xlabel="time (h)", ylabel="surface heat flux (W/(m^2 K))") # fig.savefig("test.png") plt.show() diff --git a/src/felupe/thermal/__init__.py b/src/felupe/thermal/__init__.py index 3ab10e0a5..97ea97479 100644 --- a/src/felupe/thermal/__init__.py +++ b/src/felupe/thermal/__init__.py @@ -1,5 +1,5 @@ -from ._solidbody_thermal import SolidBodyThermal from ._solidbody_surface_heat_transfer import SolidBodySurfaceHeatTransfer +from ._solidbody_thermal import SolidBodyThermal from ._solidbody_thermal_heat_flux import SolidBodyThermalHeatFlux from ._time_step import TimeStep From 19d6aa90e35198122b0a8ddfa45837654b62fe1c Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 23 Mar 2026 07:39:17 +0100 Subject: [PATCH 09/71] Update _solidbody_surface_heat_transfer.py --- src/felupe/thermal/_solidbody_surface_heat_transfer.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_heat_transfer.py b/src/felupe/thermal/_solidbody_surface_heat_transfer.py index c0e55f636..c26173edb 100644 --- a/src/felupe/thermal/_solidbody_surface_heat_transfer.py +++ b/src/felupe/thermal/_solidbody_surface_heat_transfer.py @@ -83,7 +83,7 @@ class SolidBodySurfaceHeatTransfer: ... items=[time, solid, heat_transfer], ramp=ramp, boundaries=boundaries ... ) >>> job = fem.Job(steps=[step]).evaluate( - ... filename="result.xdmf", # result file for Paraview + ... # filename="result.xdmf", # result file for Paraview ... point_data={"Temperature": lambda field, substep: temperature.values}, ... point_data_default=False, ... cell_data_default=False, @@ -114,10 +114,8 @@ def __init__(self, field, coefficient, temperature): vector=self._vector, matrix=self._matrix, multiplier=-1.0 ) - def update(self, temperature, coefficient=None): + def update(self, temperature): self.results.temperature = temperature - if coefficient is not None: - self.results.coefficient = coefficient def _vector(self, field=None, **kwargs): if field is not None: From 81462380a07ff4f22b85c0f32bdda1912be60f25 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 23 Mar 2026 07:46:39 +0100 Subject: [PATCH 10/71] Update _solidbody_surface_heat_transfer.py --- src/felupe/thermal/_solidbody_surface_heat_transfer.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_heat_transfer.py b/src/felupe/thermal/_solidbody_surface_heat_transfer.py index c26173edb..d10e3ac53 100644 --- a/src/felupe/thermal/_solidbody_surface_heat_transfer.py +++ b/src/felupe/thermal/_solidbody_surface_heat_transfer.py @@ -106,9 +106,7 @@ def __init__(self, field, coefficient, temperature): self.results = Results() self.results.temperature = temperature - - if coefficient is not None: - self.results.coefficient = coefficient + self.results.coefficient = coefficient self.assemble = Assemble( vector=self._vector, matrix=self._matrix, multiplier=-1.0 From 04eb524d8628d9376caabf6c08d40f8cb25ca454 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 23 Mar 2026 07:54:10 +0100 Subject: [PATCH 11/71] update tests for `thermal.SolidBodySurfaceHeatTransfer` --- tests/test_thermal.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/tests/test_thermal.py b/tests/test_thermal.py index 1d9e40861..f87af35e7 100644 --- a/tests/test_thermal.py +++ b/tests/test_thermal.py @@ -27,9 +27,9 @@ def test_thermal(): temperature = fem.Field(region, dim=1) field = fem.FieldContainer([temperature]) - region_convection = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) - temperature_convection = fem.Field(region_convection, dim=1) - field_convection = fem.FieldContainer([temperature_convection]) + region_heat_transfer = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) + temperature_heat_transfer = fem.Field(region_heat_transfer, dim=1) + field_heat_transfer = fem.FieldContainer([temperature_heat_transfer]) region_flux = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) temperature_flux = fem.Field(region_flux, dim=1) @@ -57,8 +57,8 @@ def test_thermal(): thermal_conductivity=1.0, # W/(m*K) ) - convection = fem.thermal.SolidBodyThermalConvection( - field=field_convection, + heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + field=field_heat_transfer, coefficient=1.0, # W/(m^2*K) temperature=10.0, # K ) @@ -70,8 +70,8 @@ def test_thermal(): solid.assemble.vector(field) solid.assemble.matrix(field) - convection.assemble.vector(field) - convection.assemble.matrix(field) + heat_transfer.assemble.vector(field) + heat_transfer.assemble.matrix(field) heat_flux.assemble.vector(field) heat_flux.assemble.matrix(field) @@ -80,11 +80,11 @@ def test_thermal(): ramp = { boundaries["right"]: 10 * table, time: 0.1 * table, - convection: 100 * table, + heat_transfer: 100 * table, heat_flux: 10 * table, } step = fem.Step( - items=[time, solid, convection, heat_flux], ramp=ramp, boundaries=boundaries + items=[time, solid, heat_transfer, heat_flux], ramp=ramp, boundaries=boundaries ) job = fem.Job(steps=[step]).evaluate( filename="result.xdmf", # result file for Paraview @@ -100,9 +100,9 @@ def test_thermal_axi(): temperature = fem.Field(region, dim=1) field = fem.FieldContainer([temperature]) - region_convection = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) - temperature_convection = fem.FieldAxisymmetric(region_convection, dim=1) - field_convection = fem.FieldContainer([temperature_convection]) + region_heat_transfer = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) + temperature_heat_transfer = fem.FieldAxisymmetric(region_heat_transfer, dim=1) + field_heat_transfer = fem.FieldContainer([temperature_heat_transfer]) region_flux = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) temperature_flux = fem.FieldAxisymmetric(region_flux, dim=1) @@ -130,8 +130,8 @@ def test_thermal_axi(): thermal_conductivity=1.0, # W/(m*K) ) - convection = fem.thermal.SolidBodyThermalConvection( - field=field_convection, + heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + field=field_heat_transfer, coefficient=1.0, # W/(m^2*K) temperature=10.0, # K ) @@ -143,8 +143,8 @@ def test_thermal_axi(): solid.assemble.vector(field) solid.assemble.matrix(field) - convection.assemble.vector(field) - convection.assemble.matrix(field) + heat_transfer.assemble.vector(field) + heat_transfer.assemble.matrix(field) heat_flux.assemble.vector(field) heat_flux.assemble.matrix(field) @@ -183,11 +183,11 @@ def test_thermal_axi(): ramp = { boundaries["right"]: 10 * table, time: 0.1 * table, - convection: 100 * table, + heat_transfer: 100 * table, heat_flux: 10 * table, } step = fem.Step( - items=[time, solid, convection, heat_flux], ramp=ramp, boundaries=boundaries + items=[time, solid, heat_transfer, heat_flux], ramp=ramp, boundaries=boundaries ) job = fem.Job(steps=[step]).evaluate( filename="result.xdmf", # result file for Paraview From fb0b177be7273f38d92b8a497bcb9a8632ad4af3 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Mon, 23 Mar 2026 07:55:11 +0100 Subject: [PATCH 12/71] Update CHANGELOG.md --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9275497ed..93a68bb39 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,7 @@ All notable changes to this project will be documented in this file. The format ### Added - Add Becker's logarithmic material model. -- Add `felupe.thermal` with `thermal.SolidBodyThermal`, `thermal.SolidBodyThermalConvection`, `thermal.SolidBodyThermalHeatFlux` and `thermal.TimeStep`. +- Add `felupe.thermal` with `thermal.SolidBodyThermal`, `thermal.SolidBodySurfaceHeatTransfer`, `thermal.SolidBodyThermalHeatFlux` and `thermal.TimeStep`. - Add `Mesh.cells_in(other_mesh, decimals=8)` to generate a cells-mask for cells, which are also included in another mesh. ### Changed From 0ea024d3c210a0f1eafaa9a859cda69e0fba2921 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Mon, 27 Apr 2026 17:19:43 +0200 Subject: [PATCH 13/71] Editorial changes to class documentation. --- .../_solidbody_surface_heat_transfer.py | 12 +++++----- .../thermal/_solidbody_surface_radiation.py | 23 ++++++++++++++----- 2 files changed, 23 insertions(+), 12 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_heat_transfer.py b/src/felupe/thermal/_solidbody_surface_heat_transfer.py index e46a85e93..70d9b61d2 100644 --- a/src/felupe/thermal/_solidbody_surface_heat_transfer.py +++ b/src/felupe/thermal/_solidbody_surface_heat_transfer.py @@ -30,17 +30,17 @@ class SolidBodySurfaceHeatTransfer: field : felupe.FieldContainer The field container with the temperature as first field. coefficient : float - The convection coefficient :math:`h` in W/(m^2 K). + The heat transfer coefficient :math:`h` in W/(m^2 K). temperature : float - The ambient temperature :math:`T_\infty` in °C. + The ambient temperature :math:`\theta_\infty` in °C. Notes ----- This class represents a boundary condition for a thermal solid body, which - is used to model heat transfer (convection, radiation) at the boundary of a - solid material. The coefficient is used to calculate the heat flux at the - boundary based on the difference between the temperature at the boundary - and the ambient temperature. + is used to model heat transfer (typically convection + radiation) at the + boundary of a solid material. The coefficient is used to calculate the heat + flux at the boundary based on the difference between the temperature at the + boundary and the ambient temperature. Examples -------- diff --git a/src/felupe/thermal/_solidbody_surface_radiation.py b/src/felupe/thermal/_solidbody_surface_radiation.py index ca946fe97..ad6546897 100644 --- a/src/felupe/thermal/_solidbody_surface_radiation.py +++ b/src/felupe/thermal/_solidbody_surface_radiation.py @@ -24,24 +24,35 @@ class SolidBodySurfaceRadiation: - r"""Radiative heat transfer on the surface of a thermal solid body. + r"""Long wave radiative heat transfer on the surface of a thermal solid body. Parameters ---------- field : felupe.FieldContainer - Field container with the temperature as first field. + Field container with the temperature in °C as first field. emissivity : float Emissivity :math:`\varepsilon` of the surface (dimensionless, :math:`0 \le \varepsilon \le 1`). temperature : float - The ambient temperature :math:`T_\infty` in °C. + The surrounding temperature :math:`\theta_{sur}` in °C. Notes ----- This class represents a boundary condition for a thermal solid body, which - is used to model radiative heat transfer at the boundary of a solid material. The - emissivity is used to calculate the heat flux at the boundary based on the - difference between the temperature at the boundary and the ambient temperature. + is used to model long wave radiative heat transfer between the boundary of a + solid material and a surrounding hemispherical body with an emissivity of + :math: `\epsilon_{sur}` = 1.0 and temperature :math:`\theta_{sur}` in °C. + + The hemispherical long wave surface emissivity :math: `\epsilon`is used to + define the corresponding solid surface property. + + The the heat flux at the boundary is calculated according to Eq. + :eq:`long-wave-radiation-flux` where 'T' denotes temperatures in K. + + .. math:: + :label: long-wave-radiation-flux + + q_r = \epsilon\,\sigma_B\,\left(T_s^4 - T_{sur}^4\right) .. note: From 909dde1453670649bee837d015717c8f473da294 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Mon, 27 Apr 2026 17:40:17 +0200 Subject: [PATCH 14/71] Add class 'SolidBodySurfaceConvection'. WIP, code save. --- .../thermal/_solidbody_surface_convection.py | 201 ++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 src/felupe/thermal/_solidbody_surface_convection.py diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py new file mode 100644 index 000000000..fc45c1658 --- /dev/null +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -0,0 +1,201 @@ +# -*- coding: utf-8 -*- +""" +This file is part of FElupe. + +FElupe is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +FElupe is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with FElupe. If not, see . +""" +import numpy as np +from scipy.constants import sigma +from scipy.sparse import csr_matrix + +from ..assembly import IntegralForm +from ..mechanics import Assemble, Results, UpdateItem + + +class SolidBodySurfaceConvection: + r"""Convective heat transfer on the surface of a thermal solid body. + + Parameters + ---------- + field : felupe.FieldContainer + Field container with the temperature in °C as first field. + convection coefficient : float | callable + Convection coefficient :math:`h_c` in W/(m^2 K). + temperature : float + The surrounding temperature :math:`\theta_{sur}` in °C. + + Notes + ----- + This class represents a boundary condition for a thermal solid body, which + is used to model convectinve heat transfer between the boundary of a + solid material and the adjacent ambient air with temperature + :math:`\theta_\infty` in °C. + + The the heat flux at the boundary is calculated according to Eq. + :eq:`convective-flux`. + + Eq. :eq: `example-horizontal-plate` gives an example for the detailed + calculation of the heat transfer coefficient for a warm horizontal plate + with heat flux upward. 'T' denotes temperatures in K. `Ra`is the Rayleigh + number, `Pr`is the Prandtl number (air) and `Nu` is the Nusselt number. + + .. math:: + :label: convective-flux + + q_c = h_c\,\left(\theta_s - \theta_\infty\right) + + .. math:: + :label: example-horizontal-plate + + \alpha = \frac{\lambda_{air}}{\rho_{air}cp_{air}} + + \T_m = 0.5\left(T_s + T_\infty\right) + + Ra = \frac{g \frac{1}{T_m}\left|\theta_s - \theta_\infty\right|\,L^3}{\alpha\mu} + + Nu = + + h_c = Nu\,\lambda_{air}\,L + + + Examples + -------- + .. pyvista-plot:: + + >>> import felupe as fem + >>> import numpy as np + >>> + >>> mesh = fem.Rectangle(n=11) + >>> region = fem.RegionQuad(mesh) + >>> temperature = fem.Field(region, dim=1, values=20.0) + >>> field = fem.FieldContainer([temperature]) + >>> + >>> region_radiation = fem.RegionQuadBoundary(mesh, mask=mesh.x == 1.0) + >>> temperature_radiation = fem.Field(region_radiation, dim=1) + >>> field_radiation = fem.FieldContainer([temperature_radiation]) + >>> + >>> boundaries = fem.BoundaryDict( + ... left=fem.Boundary(temperature, fx=0, value=20.0), + ... ) + >>> + >>> solid = fem.thermal.SolidBodyThermal( + ... field=field, + ... mass_density=1400.0, # kg / m^3 + ... specific_heat_capacity=1000.0, # J / (kg K) + ... time_step=720.0, # s + ... thermal_conductivity=1.0, # W / (m K) + ... ) + >>> radiation = fem.thermal.SolidBodySurfaceConvection( + ... field=field_radiation, + ... emissivity=0.8, + ... temperature=20.0, # °C + ... ) + >>> time = fem.thermal.TimeStep([solid]) + >>> table = fem.math.linsteps([0, 1], num=15) + >>> air_temperature = fem.math.linsteps([20, 40], num=15) # air temperature + >>> emissivity = fem.math.linsteps([0.6, 0.8], num=15) # a value between 0 ... 1 + >>> ramp = { + ... time: 18000 * table, # five hours + ... radiation["temperature"]: air_temperature, + ... radiation["emissivity"]: emissivity, + ... } + >>> step = fem.Step( + ... items=[time, solid, radiation], ramp=ramp, boundaries=boundaries + ... ) + >>> job = fem.Job(steps=[step]).evaluate() + >>> + >>> mesh.view( + ... point_data={"Temperature in °C": temperature.values} + ... ).plot("Temperature in °C").show() + + See Also + -------- + felupe.thermal.TimeStep : A time step item. + felupe.thermal.SolidBodyThermal : A thermal solid body for heat conduction. + + """ + + def __init__(self, field, emissivity, temperature): + self.field = field + self.time_step = None + + self.results = Results() + self.results.temperature = temperature # ambient temperature in °C + self.results.emissivity = emissivity + + self._sigma = sigma # Stefan-Boltzmann constant + + self.assemble = Assemble( + vector=self._vector, matrix=self._matrix, multiplier=-1.0 + ) + + def __getitem__(self, key): + return UpdateItem(self, key) + + def update(self, temperature): + self._update_temperature(temperature) + + def _update_temperature(self, temperature): + self.results.temperature = temperature + + def _update_emissivity(self, emissivity): + self.results.emissivity = emissivity + + def _vector(self, field=None, **kwargs): + if field is not None: + self.field = field + + if self.time_step is not None and self.time_step == 0: # inactive time step + return csr_matrix(([0.0], ([0], [0])), shape=(1, 1)) + + temperature = self.field.extract(grad=False)[0] + fun = [ + -self.results.emissivity + * self._sigma + * ((temperature + 273.15) ** 4 - (self.results.temperature + 273.15) ** 4) + ] + + self.results.force = IntegralForm( + fun=fun, v=self.field, dV=self.field.region.dV, grad_v=[False] + ).assemble(**kwargs) + + return self.results.force + + def _matrix(self, field=None, **kwargs): + if field is not None: + self.field = field + + if self.time_step is not None and self.time_step == 0: # inactive time step + return csr_matrix(([0.0], ([0], [0])), shape=(1, 1)) + + dim = self.field[0].dim + temperature = self.field.extract(grad=False)[0] + fun = [ + -self.results.emissivity + * self._sigma + * 4 + * (temperature + 273.15) ** 3 + * np.eye(dim).reshape(dim, dim, 1, 1) + ] + + self.results.stiffness = IntegralForm( + fun=fun, + v=self.field, + u=self.field, + dV=self.field.region.dV, + grad_v=[False], + grad_u=[False], + ).assemble(**kwargs) + + return self.results.stiffness From e47f380ea977c54b6bd67a9ee5375e6c6b116c3d Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Tue, 28 Apr 2026 09:20:32 +0200 Subject: [PATCH 15/71] Fix equation alignment. --- src/felupe/thermal/_solidbody_thermal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/felupe/thermal/_solidbody_thermal.py b/src/felupe/thermal/_solidbody_thermal.py index 0455d908b..94fef5ceb 100644 --- a/src/felupe/thermal/_solidbody_thermal.py +++ b/src/felupe/thermal/_solidbody_thermal.py @@ -68,7 +68,7 @@ class SolidBodyThermal(SolidBody): \boldsymbol{r} + \frac{\partial \boldsymbol{r}}{\partial \boldsymbol{T}} - \delta \boldsymbol {T} = \boldsymbol{0} + \delta \boldsymbol {T} &= \boldsymbol{0} \boldsymbol{K} \delta \boldsymbol{T} &= -\boldsymbol{r} From dac43c2595141a9ed25b47800478fee367c44207 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Tue, 28 Apr 2026 09:21:33 +0200 Subject: [PATCH 16/71] Update SolidBodySurfaceConvection code. Next step: test. --- .../thermal/_solidbody_surface_convection.py | 92 ++++++++++++------- 1 file changed, 60 insertions(+), 32 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index fc45c1658..e8673382c 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -30,7 +30,7 @@ class SolidBodySurfaceConvection: ---------- field : felupe.FieldContainer Field container with the temperature in °C as first field. - convection coefficient : float | callable + convection_coefficient : float | callable Convection coefficient :math:`h_c` in W/(m^2 K). temperature : float The surrounding temperature :math:`\theta_{sur}` in °C. @@ -47,8 +47,9 @@ class SolidBodySurfaceConvection: Eq. :eq: `example-horizontal-plate` gives an example for the detailed calculation of the heat transfer coefficient for a warm horizontal plate - with heat flux upward. 'T' denotes temperatures in K. `Ra`is the Rayleigh - number, `Pr`is the Prandtl number (air) and `Nu` is the Nusselt number. + with heat flux upward. `A` is the plate area, `P` is the plate perimeter + length. `T` denotes temperatures in K. `Ra` is the Rayleigh + number, `Pr` is the Prandtl number (air) and `Nu` is the Nusselt number. .. math:: :label: convective-flux @@ -58,15 +59,19 @@ class SolidBodySurfaceConvection: .. math:: :label: example-horizontal-plate - \alpha = \frac{\lambda_{air}}{\rho_{air}cp_{air}} + L &= \frac{A}{P} - \T_m = 0.5\left(T_s + T_\infty\right) + \alpha &= \frac{\lambda_\text{air}}{\rho_\text{air}cp_{air}} - Ra = \frac{g \frac{1}{T_m}\left|\theta_s - \theta_\infty\right|\,L^3}{\alpha\mu} + \T_m &= 0.5 \left(T_s + T_\infty\right) - Nu = + Ra &= \frac{g \frac{1}{T_m} \left|\theta_s - \theta_\infty\right| L^3}{\alpha\mu} - h_c = Nu\,\lambda_{air}\,L + Nu(10^4\leq Ra \leq 10^7) &= 0.54 Ra^{1/4} and + + Nu(10^7 < Ra \leq 10^11) &= 0.15 Ra^{1/3} + + h_c &= Nu \lambda_\text{air} L Examples @@ -76,17 +81,20 @@ class SolidBodySurfaceConvection: >>> import felupe as fem >>> import numpy as np >>> + >>> def hc_fun(ts, tamb): + >>> return(abs((ts-tamb)*2.0)) + >>> >>> mesh = fem.Rectangle(n=11) >>> region = fem.RegionQuad(mesh) - >>> temperature = fem.Field(region, dim=1, values=20.0) + >>> temperature = fem.Field(region, dim=1, values=30.0) >>> field = fem.FieldContainer([temperature]) >>> - >>> region_radiation = fem.RegionQuadBoundary(mesh, mask=mesh.x == 1.0) - >>> temperature_radiation = fem.Field(region_radiation, dim=1) - >>> field_radiation = fem.FieldContainer([temperature_radiation]) + >>> region_convection = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) + >>> temperature_convection = fem.Field(region_radiation, dim=1) + >>> field_convection = fem.FieldContainer([temperature_radiation]) >>> >>> boundaries = fem.BoundaryDict( - ... left=fem.Boundary(temperature, fx=0, value=20.0), + ... bottom=fem.Boundary(temperature, fy=0, value=30.0), ... ) >>> >>> solid = fem.thermal.SolidBodyThermal( @@ -96,22 +104,38 @@ class SolidBodySurfaceConvection: ... time_step=720.0, # s ... thermal_conductivity=1.0, # W / (m K) ... ) - >>> radiation = fem.thermal.SolidBodySurfaceConvection( - ... field=field_radiation, - ... emissivity=0.8, + >>> convection_constant = fem.thermal.SolidBodySurfaceConvection( + ... field=field_convection, + ... convection_coefficient=8.0, + ... temperature=20.0, # °C + ... ) + >>> convection_function = fem.thermal.SolidBodySurfaceConvection( + ... field=field_convection, + ... convection_coefficient=hc_fun, ... temperature=20.0, # °C ... ) >>> time = fem.thermal.TimeStep([solid]) >>> table = fem.math.linsteps([0, 1], num=15) - >>> air_temperature = fem.math.linsteps([20, 40], num=15) # air temperature - >>> emissivity = fem.math.linsteps([0.6, 0.8], num=15) # a value between 0 ... 1 + >>> air_temperature = fem.math.linsteps([15, 25], num=11) # air temperature + >>> ramp = { + ... time: 18000 * table, # five hours + ... convection_constant["temperature"]: air_temperature, + ... } + >>> step = fem.Step( + ... items=[time, solid, convection], ramp=ramp, boundaries=boundaries + ... ) + >>> job = fem.Job(steps=[step]).evaluate() + >>> + >>> mesh.view( + ... point_data={"Temperature in °C": temperature.values} + ... ).plot("Temperature in °C").show() + >>> >>> ramp = { ... time: 18000 * table, # five hours - ... radiation["temperature"]: air_temperature, - ... radiation["emissivity"]: emissivity, + ... convection_function["temperature"]: air_temperature, ... } >>> step = fem.Step( - ... items=[time, solid, radiation], ramp=ramp, boundaries=boundaries + ... items=[time, solid, convection_function], ramp=ramp, boundaries=boundaries ... ) >>> job = fem.Job(steps=[step]).evaluate() >>> @@ -126,13 +150,13 @@ class SolidBodySurfaceConvection: """ - def __init__(self, field, emissivity, temperature): + def __init__(self, field, convection_coefficient, temperature): self.field = field + self.convection_coefficient = convection_coefficient self.time_step = None self.results = Results() self.results.temperature = temperature # ambient temperature in °C - self.results.emissivity = emissivity self._sigma = sigma # Stefan-Boltzmann constant @@ -145,12 +169,20 @@ def __getitem__(self, key): def update(self, temperature): self._update_temperature(temperature) + self._update_convection_coefficient() # adapt hc using cur. temp. def _update_temperature(self, temperature): self.results.temperature = temperature - def _update_emissivity(self, emissivity): - self.results.emissivity = emissivity + def _update_convection_coefficient(self): + if callable(self.convection_coefficient): + self.results.convection_coefficient =\ + self.convection_coefficient( + self.results.temperature, + self.field.extract(grad=False)[0] + ) + else: + self.results.convection_coefficient = self.convection_coefficient def _vector(self, field=None, **kwargs): if field is not None: @@ -161,9 +193,8 @@ def _vector(self, field=None, **kwargs): temperature = self.field.extract(grad=False)[0] fun = [ - -self.results.emissivity - * self._sigma - * ((temperature + 273.15) ** 4 - (self.results.temperature + 273.15) ** 4) + -self.results.convection_coefficient + * (temperature - self.results.temperature) ] self.results.force = IntegralForm( @@ -182,10 +213,7 @@ def _matrix(self, field=None, **kwargs): dim = self.field[0].dim temperature = self.field.extract(grad=False)[0] fun = [ - -self.results.emissivity - * self._sigma - * 4 - * (temperature + 273.15) ** 3 + -self.results.convection_coefficient * np.eye(dim).reshape(dim, dim, 1, 1) ] From 3617ad4c29bc9cbd3c747d7965ef3720399f747a Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Tue, 28 Apr 2026 10:04:15 +0200 Subject: [PATCH 17/71] Documentation refinement. Comment out unused parameter. Test still pending. --- .../thermal/_solidbody_surface_convection.py | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index e8673382c..7d743261f 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -31,14 +31,18 @@ class SolidBodySurfaceConvection: field : felupe.FieldContainer Field container with the temperature in °C as first field. convection_coefficient : float | callable - Convection coefficient :math:`h_c` in W/(m^2 K). + Convection heat transfer coefficient :math: `h_c` in W/(m^2 K). A + callable requires the parameters 'surface temperature' and 'ambient + temperature'. + Additional parameters should be passed by **kwargs (and are not + adapted during simulation). temperature : float - The surrounding temperature :math:`\theta_{sur}` in °C. + The ambient air temperature :math:`\theta_\infty` in °C. Notes ----- This class represents a boundary condition for a thermal solid body, which - is used to model convectinve heat transfer between the boundary of a + is used to model convective heat transfer between the boundary of a solid material and the adjacent ambient air with temperature :math:`\theta_\infty` in °C. @@ -67,7 +71,7 @@ class SolidBodySurfaceConvection: Ra &= \frac{g \frac{1}{T_m} \left|\theta_s - \theta_\infty\right| L^3}{\alpha\mu} - Nu(10^4\leq Ra \leq 10^7) &= 0.54 Ra^{1/4} and + Nu(10^4\leq Ra \leq 10^7) &= 0.54 Ra^{1/4} \text{ and} Nu(10^7 < Ra \leq 10^11) &= 0.15 Ra^{1/3} @@ -122,7 +126,7 @@ class SolidBodySurfaceConvection: ... convection_constant["temperature"]: air_temperature, ... } >>> step = fem.Step( - ... items=[time, solid, convection], ramp=ramp, boundaries=boundaries + ... items=[time, solid, convection_constant], ramp=ramp, boundaries=boundaries ... ) >>> job = fem.Job(steps=[step]).evaluate() >>> @@ -143,6 +147,7 @@ class SolidBodySurfaceConvection: ... point_data={"Temperature in °C": temperature.values} ... ).plot("Temperature in °C").show() + See Also -------- felupe.thermal.TimeStep : A time step item. @@ -178,8 +183,8 @@ def _update_convection_coefficient(self): if callable(self.convection_coefficient): self.results.convection_coefficient =\ self.convection_coefficient( - self.results.temperature, - self.field.extract(grad=False)[0] + self.field.extract(grad=False)[0], # ts + self.results.temperature # tamb ) else: self.results.convection_coefficient = self.convection_coefficient @@ -211,7 +216,7 @@ def _matrix(self, field=None, **kwargs): return csr_matrix(([0.0], ([0], [0])), shape=(1, 1)) dim = self.field[0].dim - temperature = self.field.extract(grad=False)[0] + # temperature = self.field.extract(grad=False)[0] fun = [ -self.results.convection_coefficient * np.eye(dim).reshape(dim, dim, 1, 1) From 16f6c2ca4921e3072668b26e53eb4e523fcb2403 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Tue, 28 Apr 2026 14:09:26 +0200 Subject: [PATCH 18/71] Address various pytest issues (surface_convection). Now seems to run, however, newtonraphson() not converging (=> pytest error). --- src/felupe/thermal/__init__.py | 2 ++ src/felupe/thermal/_solidbody_surface_convection.py | 13 +++++++------ src/felupe/thermal/_solidbody_thermal.py | 3 ++- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/src/felupe/thermal/__init__.py b/src/felupe/thermal/__init__.py index 30eb6e1fd..9ceede790 100644 --- a/src/felupe/thermal/__init__.py +++ b/src/felupe/thermal/__init__.py @@ -1,5 +1,6 @@ from ._solidbody_heat_flux import SolidBodyHeatFlux from ._solidbody_surface_heat_transfer import SolidBodySurfaceHeatTransfer +from ._solidbody_surface_convection import SolidBodySurfaceConvection from ._solidbody_surface_radiation import SolidBodySurfaceRadiation from ._solidbody_thermal import SolidBodyThermal from ._time_step import TimeStep @@ -7,6 +8,7 @@ __all__ = [ "SolidBodyThermal", "SolidBodySurfaceHeatTransfer", + "SolidBodySurfaceConvection", "SolidBodySurfaceRadiation", "SolidBodyHeatFlux", "TimeStep", diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 7d743261f..e8352548e 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -86,7 +86,7 @@ class SolidBodySurfaceConvection: >>> import numpy as np >>> >>> def hc_fun(ts, tamb): - >>> return(abs((ts-tamb)*2.0)) + ... return((ts-tamb)*2.0) >>> >>> mesh = fem.Rectangle(n=11) >>> region = fem.RegionQuad(mesh) @@ -94,8 +94,8 @@ class SolidBodySurfaceConvection: >>> field = fem.FieldContainer([temperature]) >>> >>> region_convection = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) - >>> temperature_convection = fem.Field(region_radiation, dim=1) - >>> field_convection = fem.FieldContainer([temperature_radiation]) + >>> temperature_convection = fem.Field(region_convection, dim=1) + >>> field_convection = fem.FieldContainer([temperature_convection]) >>> >>> boundaries = fem.BoundaryDict( ... bottom=fem.Boundary(temperature, fy=0, value=30.0), @@ -110,7 +110,7 @@ class SolidBodySurfaceConvection: ... ) >>> convection_constant = fem.thermal.SolidBodySurfaceConvection( ... field=field_convection, - ... convection_coefficient=8.0, + ... convection_coefficient=8.0, # W/(m^2 K) ... temperature=20.0, # °C ... ) >>> convection_function = fem.thermal.SolidBodySurfaceConvection( @@ -119,8 +119,8 @@ class SolidBodySurfaceConvection: ... temperature=20.0, # °C ... ) >>> time = fem.thermal.TimeStep([solid]) - >>> table = fem.math.linsteps([0, 1], num=15) - >>> air_temperature = fem.math.linsteps([15, 25], num=11) # air temperature + >>> table = fem.math.linsteps([0, 1], num=25) + >>> air_temperature = fem.math.linsteps([15, 25], num=10) # air temperature >>> ramp = { ... time: 18000 * table, # five hours ... convection_constant["temperature"]: air_temperature, @@ -162,6 +162,7 @@ def __init__(self, field, convection_coefficient, temperature): self.results = Results() self.results.temperature = temperature # ambient temperature in °C + self.results.convection_coefficient = convection_coefficient self._sigma = sigma # Stefan-Boltzmann constant diff --git a/src/felupe/thermal/_solidbody_thermal.py b/src/felupe/thermal/_solidbody_thermal.py index 94fef5ceb..edef6d85f 100644 --- a/src/felupe/thermal/_solidbody_thermal.py +++ b/src/felupe/thermal/_solidbody_thermal.py @@ -132,7 +132,8 @@ class SolidBodyThermal(SolidBody): See Also -------- felupe.thermal.TimeStep : A time step item. - felupe.thermal.SolidBodySurfaceHeatTransfer : A surface heat transfer boundary condition. + felupe.thermal.SolidBodySurfaceHeatTransfer : A general surface heat transfer boundary condition. + felupe.thermal.SolidBodySurfaceConvection : A thermal convection boundary condition. felupe.thermal.SolidBodySurfaceRadiation : A thermal radiation boundary condition. felupe.thermal.SolidBodyHeatFlux : A thermal heat flux boundary condition. From 4437d56aacc1b54f14df0430e5144568404315d8 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Tue, 28 Apr 2026 16:43:22 +0200 Subject: [PATCH 19/71] Remove 'sigma'. Various further doc edits. --- src/felupe/thermal/_solidbody_surface_convection.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index e8352548e..95aabdc9a 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -16,7 +16,6 @@ along with FElupe. If not, see . """ import numpy as np -from scipy.constants import sigma from scipy.sparse import csr_matrix from ..assembly import IntegralForm @@ -88,7 +87,7 @@ class SolidBodySurfaceConvection: >>> def hc_fun(ts, tamb): ... return((ts-tamb)*2.0) >>> - >>> mesh = fem.Rectangle(n=11) + >>> mesh = fem.Rectangle(n=11) # rectangle w/ 10x10 cells >>> region = fem.RegionQuad(mesh) >>> temperature = fem.Field(region, dim=1, values=30.0) >>> field = fem.FieldContainer([temperature]) @@ -110,7 +109,7 @@ class SolidBodySurfaceConvection: ... ) >>> convection_constant = fem.thermal.SolidBodySurfaceConvection( ... field=field_convection, - ... convection_coefficient=8.0, # W/(m^2 K) + ... convection_coefficient=5.0, # W/(m^2 K) ... temperature=20.0, # °C ... ) >>> convection_function = fem.thermal.SolidBodySurfaceConvection( @@ -120,10 +119,10 @@ class SolidBodySurfaceConvection: ... ) >>> time = fem.thermal.TimeStep([solid]) >>> table = fem.math.linsteps([0, 1], num=25) - >>> air_temperature = fem.math.linsteps([15, 25], num=10) # air temperature + >>> air_temperature = fem.math.linsteps([15, 25], num=10) >>> ramp = { ... time: 18000 * table, # five hours - ... convection_constant["temperature"]: air_temperature, + ... convection_constant["temperature"]: air_temperature ... } >>> step = fem.Step( ... items=[time, solid, convection_constant], ramp=ramp, boundaries=boundaries @@ -164,8 +163,6 @@ def __init__(self, field, convection_coefficient, temperature): self.results.temperature = temperature # ambient temperature in °C self.results.convection_coefficient = convection_coefficient - self._sigma = sigma # Stefan-Boltzmann constant - self.assemble = Assemble( vector=self._vector, matrix=self._matrix, multiplier=-1.0 ) From 4a43cf5bc8e62bc70dedc521fa3e3a1a0d2f8b39 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Tue, 5 May 2026 18:03:43 +0200 Subject: [PATCH 20/71] Fix bug - all ramp lists seem to necessarily be of the same length. Duh. --- src/felupe/thermal/_solidbody_surface_convection.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 95aabdc9a..81dfbde24 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -87,12 +87,12 @@ class SolidBodySurfaceConvection: >>> def hc_fun(ts, tamb): ... return((ts-tamb)*2.0) >>> - >>> mesh = fem.Rectangle(n=11) # rectangle w/ 10x10 cells + >>> mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells >>> region = fem.RegionQuad(mesh) >>> temperature = fem.Field(region, dim=1, values=30.0) >>> field = fem.FieldContainer([temperature]) >>> - >>> region_convection = fem.RegionQuadBoundary(mesh, mask=mesh.y == 1.0) + >>> region_convection = fem.RegionQuadBoundary(mesh, mask=mesh.y == 0.25) >>> temperature_convection = fem.Field(region_convection, dim=1) >>> field_convection = fem.FieldContainer([temperature_convection]) >>> @@ -118,7 +118,7 @@ class SolidBodySurfaceConvection: ... temperature=20.0, # °C ... ) >>> time = fem.thermal.TimeStep([solid]) - >>> table = fem.math.linsteps([0, 1], num=25) + >>> table = fem.math.linsteps([0, 1], num=10) >>> air_temperature = fem.math.linsteps([15, 25], num=10) >>> ramp = { ... time: 18000 * table, # five hours @@ -128,11 +128,12 @@ class SolidBodySurfaceConvection: ... items=[time, solid, convection_constant], ramp=ramp, boundaries=boundaries ... ) >>> job = fem.Job(steps=[step]).evaluate() - >>> + >>> ... >>> mesh.view( ... point_data={"Temperature in °C": temperature.values} ... ).plot("Temperature in °C").show() >>> + >>> time = fem.thermal.TimeStep([solid]) >>> ramp = { ... time: 18000 * table, # five hours ... convection_function["temperature"]: air_temperature, From 8701f4eea6226b9e2602edb9dae93183d6711a9d Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 09:35:02 +0200 Subject: [PATCH 21/71] Adapt documentation. Adapt example to actually calculate h_c for a horizontal plate w/ heat flow upwards. Not working yet, the function call does not seem to work at runtime. --- .../thermal/_solidbody_surface_convection.py | 67 ++++++++++++------- 1 file changed, 42 insertions(+), 25 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 81dfbde24..56e3eb7d6 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -70,22 +70,35 @@ class SolidBodySurfaceConvection: Ra &= \frac{g \frac{1}{T_m} \left|\theta_s - \theta_\infty\right| L^3}{\alpha\mu} - Nu(10^4\leq Ra \leq 10^7) &= 0.54 Ra^{1/4} \text{ and} + Nu(10^4\leq Ra \leq 10^7) &= 0.54 Ra^{1/4} \text{(Pr > 0.7) and} Nu(10^7 < Ra \leq 10^11) &= 0.15 Ra^{1/3} - h_c &= Nu \lambda_\text{air} L + h_c &= \frac{Nu \lambda_\text{air}}{L} Examples -------- .. pyvista-plot:: + >>> import math >>> import felupe as fem >>> import numpy as np >>> >>> def hc_fun(ts, tamb): - ... return((ts-tamb)*2.0) + ... l = 0.25 # slab 1 x 1 m^2 + ... alpha = 2.25E-05 # m^2/s, air at 300 K + ... lam_air = 0.0263 # W/(m K), air at 300 K + ... pr = 0.707 # air at 300 K + ... t_m = 0.5*(ts + tamb) + 273.15 # K + ... ra = (9.81 / t_m * abs(ts - tamb) * l**3)/alpha/1.59E-7 + ... if (10**4 <= ra <= 10**7) and pr > 0.7: + ... nu = 0.54 * math.pow(ra, 0.25) + ... elif (10**7 < ra <= 10**11): + ... nu = 0.15 * math.pow(ra, 0.33) + ... else: + ... nu = 1 + ... return(nu*lam_air/l) >>> >>> mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells >>> region = fem.RegionQuad(mesh) @@ -107,11 +120,6 @@ class SolidBodySurfaceConvection: ... time_step=720.0, # s ... thermal_conductivity=1.0, # W / (m K) ... ) - >>> convection_constant = fem.thermal.SolidBodySurfaceConvection( - ... field=field_convection, - ... convection_coefficient=5.0, # W/(m^2 K) - ... temperature=20.0, # °C - ... ) >>> convection_function = fem.thermal.SolidBodySurfaceConvection( ... field=field_convection, ... convection_coefficient=hc_fun, @@ -122,20 +130,6 @@ class SolidBodySurfaceConvection: >>> air_temperature = fem.math.linsteps([15, 25], num=10) >>> ramp = { ... time: 18000 * table, # five hours - ... convection_constant["temperature"]: air_temperature - ... } - >>> step = fem.Step( - ... items=[time, solid, convection_constant], ramp=ramp, boundaries=boundaries - ... ) - >>> job = fem.Job(steps=[step]).evaluate() - >>> ... - >>> mesh.view( - ... point_data={"Temperature in °C": temperature.values} - ... ).plot("Temperature in °C").show() - >>> - >>> time = fem.thermal.TimeStep([solid]) - >>> ramp = { - ... time: 18000 * table, # five hours ... convection_function["temperature"]: air_temperature, ... } >>> step = fem.Step( @@ -144,8 +138,8 @@ class SolidBodySurfaceConvection: >>> job = fem.Job(steps=[step]).evaluate() >>> >>> mesh.view( - ... point_data={"Temperature in °C": temperature.values} - ... ).plot("Temperature in °C").show() + ... point_data={"Temperature 2 in °C": temperature.values} + ... ).plot("Temperature 2 in °C").show() See Also @@ -162,7 +156,14 @@ def __init__(self, field, convection_coefficient, temperature): self.results = Results() self.results.temperature = temperature # ambient temperature in °C - self.results.convection_coefficient = convection_coefficient + if callable(convection_coefficient): + self.results.convection_coefficient =\ + convection_coefficient( + self.field.extract(grad=False)[0], # ts + temperature # tamb + ) + else: + self.results.convection_coefficient = convection_coefficient self.assemble = Assemble( vector=self._vector, matrix=self._matrix, multiplier=-1.0 @@ -231,3 +232,19 @@ def _matrix(self, field=None, **kwargs): ).assemble(**kwargs) return self.results.stiffness + +# import math + +# def hc_fun(ts, tamb): +# l = 0.25 # slab 1 x 1 m^2 +# alpha = 2.25E-05 # m^2/s, air at 300 K +# pr = 0.707 # air at 300 K +# t_m = 0.5*(ts + tamb) + 273.15 # K +# ra = (9.81 / t_m * abs(ts - tamb) * l**3)/alpha/1.59E-5 +# if (10**4 <= ra <= 10**7) and pr > 0.7: +# nu = 0.54 * math.pow(ra, 0.25) +# elif (10**7 < ra <= 10**11): +# nu = 0.15 * math.pow(ra, 0.33) +# else: +# nu = 1 +# return(nu*0.0263/l) From 3a748991ca914eac1992cc8dac7fedd1f80c5745 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 11:06:41 +0200 Subject: [PATCH 22/71] Variant of hc_fun inclusion. Also does not work. --- .../thermal/_solidbody_surface_convection.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 56e3eb7d6..2bbdab8c3 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -128,9 +128,11 @@ class SolidBodySurfaceConvection: >>> time = fem.thermal.TimeStep([solid]) >>> table = fem.math.linsteps([0, 1], num=10) >>> air_temperature = fem.math.linsteps([15, 25], num=10) + >>> hc_const = fem.math.linsteps([4, 4], num=10) >>> ramp = { ... time: 18000 * table, # five hours ... convection_function["temperature"]: air_temperature, + ... convection_function["convection_coefficient"]: hc_const, ... } >>> step = fem.Step( ... items=[time, solid, convection_function], ramp=ramp, boundaries=boundaries @@ -151,11 +153,12 @@ class SolidBodySurfaceConvection: def __init__(self, field, convection_coefficient, temperature): self.field = field - self.convection_coefficient = convection_coefficient + self.convection_coefficient = convection_coefficient # value or callable self.time_step = None self.results = Results() self.results.temperature = temperature # ambient temperature in °C + if callable(convection_coefficient): self.results.convection_coefficient =\ convection_coefficient( @@ -174,20 +177,20 @@ def __getitem__(self, key): def update(self, temperature): self._update_temperature(temperature) - self._update_convection_coefficient() # adapt hc using cur. temp. + # self._update_convection_coefficient() # adapt hc using cur. temp. def _update_temperature(self, temperature): self.results.temperature = temperature - def _update_convection_coefficient(self): - if callable(self.convection_coefficient): + def _update_convection_coefficient(self, convection_coefficient): + if callable(convection_coefficient): self.results.convection_coefficient =\ - self.convection_coefficient( + convection_coefficient( self.field.extract(grad=False)[0], # ts self.results.temperature # tamb ) else: - self.results.convection_coefficient = self.convection_coefficient + self.results.convection_coefficient = convection_coefficient def _vector(self, field=None, **kwargs): if field is not None: From a763092d8dcb49b641cfa3cdf9aced9aab1628a6 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 11:52:39 +0200 Subject: [PATCH 23/71] Refine example documentation. --- .../thermal/_solidbody_surface_convection.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 2bbdab8c3..d1f6146ab 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -49,10 +49,12 @@ class SolidBodySurfaceConvection: :eq:`convective-flux`. Eq. :eq: `example-horizontal-plate` gives an example for the detailed - calculation of the heat transfer coefficient for a warm horizontal plate - with heat flux upward. `A` is the plate area, `P` is the plate perimeter - length. `T` denotes temperatures in K. `Ra` is the Rayleigh - number, `Pr` is the Prandtl number (air) and `Nu` is the Nusselt number. + calculation of the convective heat transfer coefficient for a warm + horizontal plate with heat flux upward. `A` is the plate area, `P` is the + plate perimeter length. `T` denotes temperatures in K, :math:`T_m` is + the 'film temperature' for which fluid properties are evaluated. `Ra` is + the Rayleigh number, `Pr` is the Prandtl number (air) and `Nu` is the + Nusselt number. .. math:: :label: convective-flux @@ -64,9 +66,9 @@ class SolidBodySurfaceConvection: L &= \frac{A}{P} - \alpha &= \frac{\lambda_\text{air}}{\rho_\text{air}cp_{air}} + T_m &= 0.5 \left(T_s + T_\infty\right) - \T_m &= 0.5 \left(T_s + T_\infty\right) + \alpha\left(T_m\right) &= \frac{\lambda_\text{air}}{\rho_\text{air}cp_{air}} Ra &= \frac{g \frac{1}{T_m} \left|\theta_s - \theta_\infty\right| L^3}{\alpha\mu} From 5d7dad6374081330e3735ae7f8b599e4f509aea2 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 11:58:56 +0200 Subject: [PATCH 24/71] Change "Ra over the limit" calculation of Nu to constant "Ra limit". --- src/felupe/thermal/_solidbody_surface_convection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index d1f6146ab..a3c8d5bb9 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -99,7 +99,7 @@ class SolidBodySurfaceConvection: ... elif (10**7 < ra <= 10**11): ... nu = 0.15 * math.pow(ra, 0.33) ... else: - ... nu = 1 + ... nu = 0.15 * math.pow(10**11, 0.33) ... return(nu*lam_air/l) >>> >>> mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells From 166812e9d56330776d803e71c760675c51d1c99a Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 12:01:59 +0200 Subject: [PATCH 25/71] Adapt Ra limit values to engineering notation. --- src/felupe/thermal/_solidbody_surface_convection.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index a3c8d5bb9..e6a5c0cda 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -94,12 +94,12 @@ class SolidBodySurfaceConvection: ... pr = 0.707 # air at 300 K ... t_m = 0.5*(ts + tamb) + 273.15 # K ... ra = (9.81 / t_m * abs(ts - tamb) * l**3)/alpha/1.59E-7 - ... if (10**4 <= ra <= 10**7) and pr > 0.7: + ... if (1E04 <= ra <= 1E07) and pr > 0.7: ... nu = 0.54 * math.pow(ra, 0.25) - ... elif (10**7 < ra <= 10**11): + ... elif (1E07 < ra <= 1E11): ... nu = 0.15 * math.pow(ra, 0.33) ... else: - ... nu = 0.15 * math.pow(10**11, 0.33) + ... nu = 0.15 * math.pow(1E11, 0.33) ... return(nu*lam_air/l) >>> >>> mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells From 962a68ad76bcd2f3a558505260ebb862f6fdf9fb Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 12:51:44 +0200 Subject: [PATCH 26/71] Add convection class to documentation (API). Update example description after first successful render. --- docs/felupe/thermal.rst | 6 ++++++ src/felupe/thermal/_solidbody_surface_convection.py | 4 ++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/docs/felupe/thermal.rst b/docs/felupe/thermal.rst index 4bbec539b..1a7a9842b 100644 --- a/docs/felupe/thermal.rst +++ b/docs/felupe/thermal.rst @@ -16,6 +16,7 @@ Thermal thermal.SolidBodyThermal thermal.SolidBodySurfaceHeatTransfer + thermal.SolidBodySurfaceConvection thermal.SolidBodySurfaceRadiation thermal.SolidBodyHeatFlux @@ -36,6 +37,11 @@ Thermal :undoc-members: :show-inheritance: +.. autoclass:: felupe.thermal.SolidBodySurfaceConvection + :members: + :undoc-members: + :show-inheritance: + .. autoclass:: felupe.thermal.SolidBodySurfaceRadiation :members: :undoc-members: diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index e6a5c0cda..1ab8c3914 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -68,9 +68,9 @@ class SolidBodySurfaceConvection: T_m &= 0.5 \left(T_s + T_\infty\right) - \alpha\left(T_m\right) &= \frac{\lambda_\text{air}}{\rho_\text{air}cp_{air}} + \alpha\left(T_m\right) &= \frac{\lambda_\text{air}}{\rho_\text{air}\,cp_{air}} - Ra &= \frac{g \frac{1}{T_m} \left|\theta_s - \theta_\infty\right| L^3}{\alpha\mu} + Ra &= \frac{g \frac{1}{T_m} \left|\theta_s - \theta_\infty\right| L^3}{\alpha\nu} Nu(10^4\leq Ra \leq 10^7) &= 0.54 Ra^{1/4} \text{(Pr > 0.7) and} From 0dec89a069b9f7c9df7aaae8c3f67ba81868c6c6 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 12:53:11 +0200 Subject: [PATCH 27/71] Update .gitignore. --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 80af874f5..7bc3083a0 100644 --- a/.gitignore +++ b/.gitignore @@ -122,6 +122,7 @@ examples/result.xdmf docs/examples/ docs/tutorial/* !docs/tutorial/examples +docs/mesh.png sg_execution_times.rst result.h5 result.xdmf From 06360caab7a809ae0e9a17efb4ce4c04a62a29f7 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 13:02:03 +0200 Subject: [PATCH 28/71] Further update example description after renewed render. --- src/felupe/thermal/_solidbody_surface_convection.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 1ab8c3914..59471d418 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -48,7 +48,7 @@ class SolidBodySurfaceConvection: The the heat flux at the boundary is calculated according to Eq. :eq:`convective-flux`. - Eq. :eq: `example-horizontal-plate` gives an example for the detailed + Eq. :eq:`example-horizontal-plate` gives an example for the detailed calculation of the convective heat transfer coefficient for a warm horizontal plate with heat flux upward. `A` is the plate area, `P` is the plate perimeter length. `T` denotes temperatures in K, :math:`T_m` is @@ -68,15 +68,15 @@ class SolidBodySurfaceConvection: T_m &= 0.5 \left(T_s + T_\infty\right) - \alpha\left(T_m\right) &= \frac{\lambda_\text{air}}{\rho_\text{air}\,cp_{air}} + \alpha\left(T_m\right) &= \frac{\lambda_\text{air}}{\rho_\text{air}\,c_{p,air}} - Ra &= \frac{g \frac{1}{T_m} \left|\theta_s - \theta_\infty\right| L^3}{\alpha\nu} + Ra &= \frac{g \frac{1}{T_m} \left|\theta_s - \theta_\infty\right| L^3}{\alpha\,\nu} - Nu(10^4\leq Ra \leq 10^7) &= 0.54 Ra^{1/4} \text{(Pr > 0.7) and} + Nu(10^4\leq Ra \leq 10^7) &= 0.54\,Ra^{1/4} \text{ (Pr > 0.7) and} - Nu(10^7 < Ra \leq 10^11) &= 0.15 Ra^{1/3} + Nu(10^7 < Ra \leq 10^11) &= 0.15\,Ra^{1/3} - h_c &= \frac{Nu \lambda_\text{air}}{L} + h_c &= \frac{Nu\,\lambda_\text{air}}{L} Examples From 593dab923d975e2990c43af71a6d225d768e7b92 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 13:18:49 +0200 Subject: [PATCH 29/71] Yet another update of example doc. --- src/felupe/thermal/_solidbody_surface_convection.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 59471d418..2059eeceb 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -29,12 +29,10 @@ class SolidBodySurfaceConvection: ---------- field : felupe.FieldContainer Field container with the temperature in °C as first field. - convection_coefficient : float | callable - Convection heat transfer coefficient :math: `h_c` in W/(m^2 K). A + convection_coefficient : float or callable + Convection heat transfer coefficient :math:`h_c` in W/(m^2 K). A callable requires the parameters 'surface temperature' and 'ambient temperature'. - Additional parameters should be passed by **kwargs (and are not - adapted during simulation). temperature : float The ambient air temperature :math:`\theta_\infty` in °C. From 08d343d4d4d727c21575cae7fc8d6dc124c832d2 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 13:50:09 +0200 Subject: [PATCH 30/71] In this version, function call seems to work. Testing required. --- .../thermal/_solidbody_surface_convection.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 2059eeceb..2a6e7f36d 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -122,22 +122,20 @@ class SolidBodySurfaceConvection: ... ) >>> convection_function = fem.thermal.SolidBodySurfaceConvection( ... field=field_convection, - ... convection_coefficient=hc_fun, + ... convection_coefficient=hc_fun(30, 20), ... temperature=20.0, # °C ... ) >>> time = fem.thermal.TimeStep([solid]) >>> table = fem.math.linsteps([0, 1], num=10) >>> air_temperature = fem.math.linsteps([15, 25], num=10) - >>> hc_const = fem.math.linsteps([4, 4], num=10) >>> ramp = { ... time: 18000 * table, # five hours ... convection_function["temperature"]: air_temperature, - ... convection_function["convection_coefficient"]: hc_const, ... } >>> step = fem.Step( ... items=[time, solid, convection_function], ramp=ramp, boundaries=boundaries ... ) - >>> job = fem.Job(steps=[step]).evaluate() + >>> job = fem.Job(steps=[step]).evaluate(verbose=False) >>> >>> mesh.view( ... point_data={"Temperature 2 in °C": temperature.values} @@ -177,20 +175,20 @@ def __getitem__(self, key): def update(self, temperature): self._update_temperature(temperature) - # self._update_convection_coefficient() # adapt hc using cur. temp. + self._update_convection_coefficient() # adapt hc using cur. temp. def _update_temperature(self, temperature): self.results.temperature = temperature - def _update_convection_coefficient(self, convection_coefficient): - if callable(convection_coefficient): + def _update_convection_coefficient(self): + if callable(self.convection_coefficient): self.results.convection_coefficient =\ - convection_coefficient( + self.convection_coefficient( self.field.extract(grad=False)[0], # ts self.results.temperature # tamb ) else: - self.results.convection_coefficient = convection_coefficient + self.results.convection_coefficient = self.convection_coefficient def _vector(self, field=None, **kwargs): if field is not None: From c4eea2a7890b0ad477ee1444050319746473f393 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 14:28:33 +0200 Subject: [PATCH 31/71] Refactor the example. The function call now really seems to be working(!). It would be good to be able to add a h_c.vs.time graph to check ... (callback?). --- .../thermal/_solidbody_surface_convection.py | 59 ++++++++++++++++++- 1 file changed, 57 insertions(+), 2 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 2a6e7f36d..49093790f 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -79,11 +79,14 @@ class SolidBodySurfaceConvection: Examples -------- + The examples here show how to use the :class:`~felupe.SolidBodySurfaceConvection` + both with a constant value for :math:`h_c` and a function + :math:`h_c=f(\theta_s, \theta_\infty)` according to + .. pyvista-plot:: + :context: >>> import math - >>> import felupe as fem - >>> import numpy as np >>> >>> def hc_fun(ts, tamb): ... l = 0.25 # slab 1 x 1 m^2 @@ -99,6 +102,14 @@ class SolidBodySurfaceConvection: ... else: ... nu = 0.15 * math.pow(1E11, 0.33) ... return(nu*lam_air/l) + + Set up the model (a horizontal slab with dimensions 1 x 1 m^2, 0.25 m thick). + + .. pyvista-plot:: + :context: + + >>> import felupe as fem + >>> import numpy as np >>> >>> mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells >>> region = fem.RegionQuad(mesh) @@ -120,6 +131,43 @@ class SolidBodySurfaceConvection: ... time_step=720.0, # s ... thermal_conductivity=1.0, # W / (m K) ... ) + + We will start with the example using a constant value, this is basically + identical in functionality to :class:`~felupe.SolidBodySurfaceHeatTransfer` + when the value entered for `coefficient` corresponds to the convective + part, only. + + .. pyvista-plot:: + :context: + :force_static: + + >>> convection_constant = fem.thermal.SolidBodySurfaceConvection( + ... field=field_convection, + ... convection_coefficient=5.0, + ... temperature=20.0, # °C + ... ) + >>> time = fem.thermal.TimeStep([solid]) + >>> table = fem.math.linsteps([0, 1], num=10) + >>> air_temperature = fem.math.linsteps([15, 25], num=10) + >>> ramp = { + ... time: 18000 * table, # five hours + ... convection_constant["temperature"]: air_temperature, + ... } + >>> step = fem.Step( + ... items=[time, solid, convection_constant], ramp=ramp, boundaries=boundaries + ... ) + >>> job = fem.Job(steps=[step]).evaluate(verbose=False) + >>> + >>> mesh.view( + ... point_data={"Temperature in °C": temperature.values} + ... ).plot("Temperature in °C").show() + + And now set up convection to use the function for :math:`h_c` defined above + using the same air temperature boundary conditions ... + + .. pyvista-plot:: + :context: + >>> convection_function = fem.thermal.SolidBodySurfaceConvection( ... field=field_convection, ... convection_coefficient=hc_fun(30, 20), @@ -132,6 +180,13 @@ class SolidBodySurfaceConvection: ... time: 18000 * table, # five hours ... convection_function["temperature"]: air_temperature, ... } + + ... and run. + + .. pyvista-plot:: + :context: + :force_static: + >>> step = fem.Step( ... items=[time, solid, convection_function], ramp=ramp, boundaries=boundaries ... ) From d4f4d38013e0151cd64dff37bb969d104f5632d4 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 6 May 2026 12:53:11 +0200 Subject: [PATCH 32/71] Update .gitignore. (cherry picked from commit 0dec89a069b9f7c9df7aaae8c3f67ba81868c6c6) --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 80af874f5..7bc3083a0 100644 --- a/.gitignore +++ b/.gitignore @@ -122,6 +122,7 @@ examples/result.xdmf docs/examples/ docs/tutorial/* !docs/tutorial/examples +docs/mesh.png sg_execution_times.rst result.h5 result.xdmf From b04c6e55c92ce765fa28463a45a747877112b6d9 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Thu, 7 May 2026 08:20:39 +0200 Subject: [PATCH 33/71] Add reference for horizontal plate convection coefficient. --- src/felupe/thermal/_solidbody_surface_convection.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 49093790f..5d721dcfe 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -46,7 +46,7 @@ class SolidBodySurfaceConvection: The the heat flux at the boundary is calculated according to Eq. :eq:`convective-flux`. - Eq. :eq:`example-horizontal-plate` gives an example for the detailed + Eq. :eq:`example-horizontal-plate` [1]_ gives an example for the detailed calculation of the convective heat transfer coefficient for a warm horizontal plate with heat flux upward. `A` is the plate area, `P` is the plate perimeter length. `T` denotes temperatures in K, :math:`T_m` is @@ -196,6 +196,11 @@ class SolidBodySurfaceConvection: ... point_data={"Temperature 2 in °C": temperature.values} ... ).plot("Temperature 2 in °C").show() + References + ---------- + .. [1] F. P. Incropera, D. P. DeWitt, and et. al., Fundamentals of Heat + and Mass Transfer, 6th Edition. John Wiley & Sons, 2007; + ISBN 0-471-45728-0. See Also -------- From d8a642756d10f08209d0cf13d5be9d266053fda2 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Thu, 7 May 2026 16:50:26 +0200 Subject: [PATCH 34/71] Add tester function code for _convection example. Can be deleted after testing complete. --- src/felupe/thermal/_sbsc_tester.py | 169 +++++++++++++++++++++++++++++ 1 file changed, 169 insertions(+) create mode 100644 src/felupe/thermal/_sbsc_tester.py diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py new file mode 100644 index 000000000..7c50f100c --- /dev/null +++ b/src/felupe/thermal/_sbsc_tester.py @@ -0,0 +1,169 @@ +# -*- coding: utf-8 -*- +""" +Tester version for _solid_body_surface_convection.py example case. +""" +import math + +from pyfluids import HumidAir, InputHumidAir +from scipy.constants import g + +import felupe as fem + + +T0 = 273.15 +P0 = 101325 + +def pyfluids_units(): + check = HumidAir().factory() + if str(check.units_system) == 'SIWithCelsiusAndPercents': + dt_ = 0 # use °C + rh_ = 1 # use % + else: + dt_ = 273.15 # use K + rh_ = 100 # use absolute value + return dt_, rh_ + + +def rayleigh(ts_c, ti_c, length_, rh=10): + """ + Calculate dimensionless Rayleigh number Ra. + """ + dtk, rhf = pyfluids_units() + tm_c = (ts_c + ti_c)/2 + + # Humid air properties at p0 and Tm (indoors). + air = HumidAir().with_state( + InputHumidAir.pressure(P0), + InputHumidAir.temperature(tm_c + dtk), + InputHumidAir.relative_humidity(rh/rhf), + ) + rho = air.density + cp = air.specific_heat + uv = air.kinematic_viscosity # m^2/s + k = air.conductivity # W/(m K) + alpha = k/(rho*cp) # m^2/s thermal diffusivity + beta = 1/(tm_c + T0) + + # Eqn. 9.25, page 571. + ra = g*beta*abs(ts_c - ti_c)*length_*length_*length_/alpha/uv + + return ra + +def nusselt_horizontal(ra, pr, hflux='z+'): + """ + Calculate dimensionless Nusselt number Nu for horizontal plates for various + cases of heat flux direction. + """ + if hflux == 'z+': # warm plate, top face or cold plate, bottom face + if ra < 1E04: + nu = 0.54*math.pow(1E04,0.25) + elif (1E04 <= ra <= 1E07) and (pr >= 0.7): + nu = 0.54*math.pow(ra,0.25) + elif 1E07 < ra <= 1E11: + nu = 0.15*math.pow(ra,0.33333) + else: + nu = 0.15*math.pow(1E11,0.33333) + else: # warm plate, bottom face or cold plate, top face + if ra < 1E04: + nu = 0.52*math.pow(1E04,0.2) + elif 1E04 <= ra <= 1E09 and pr >= 0.7: + nu = 0.52*math.pow(ra,0.2) + else: + nu = 0.52*math.pow(1E09,0.2) + return nu + + +def hc_fun(ts, tamb): + """ + Calculate convection coefficient for horizontal plate. + """ + dtk, rhf = pyfluids_units() + tm_c = (ts + tamb)/2 + rh = 50 + + air = HumidAir().with_state( + InputHumidAir.pressure(P0), + InputHumidAir.temperature(tm_c + dtk), + InputHumidAir.relative_humidity(rh/rhf), + ) + + l = 1.25 # slab 1 x 1 m^2 + ra = rayleigh(ts, tamb, l) + # alpha = 2.25E-05 # m^2/s, air at 300 K + # lam_air = 0.0263 # W/(m K), air at 300 K + # pr = 0.707 # air at 300 K + # t_m = 0.5*(ts + tamb) + 273.15 # K + # ra = (9.81 / t_m * abs(ts - tamb) * math.pow(l, 3.0))/alpha/1.59E-7 + if ts > tamb: + nu = nusselt_horizontal(ra, air.prandtl, hflux='z+') + else: + nu = nusselt_horizontal(ra, air.prandtl, hflux='z-') + + return(nu*air.conductivity/l) + +mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells +region = fem.RegionQuad(mesh) +temperature = fem.Field(region, dim=1, values=30.0) +field = fem.FieldContainer([temperature]) + +region_convection = fem.RegionQuadBoundary(mesh, mask=mesh.y == 0.25) +temperature_convection = fem.Field(region_convection, dim=1) +field_convection = fem.FieldContainer([temperature_convection]) + +boundaries = fem.BoundaryDict( + bottom=fem.Boundary(temperature, fy=0, value=30.0), +) + +solid = fem.thermal.SolidBodyThermal( + field=field, + mass_density=1400.0, # kg / m^3 + specific_heat_capacity=1000.0, # J / (kg K) + time_step=720.0, # s + thermal_conductivity=1.0, # W / (m K) +) + +convection_function = fem.thermal.SolidBodySurfaceConvection( + field=field_convection, + convection_coefficient=hc_fun(30, 20), + temperature=20.0, # °C +) + +def callback(stepnumber, substepnumber, substep, flux_data): + """Save mean surface heat flux at top (convective) boundary. + """ + tamb = list(ramp.values())[1][substepnumber] + ts = list(convection_function.field.extract(grad=False)[0])[0][0:1][0].mean() + heat_flux = solid.heat_flux_boundary + qc = heat_flux(region=region_convection) + flux_data["top.W.m-2"].append(qc) + flux_data["tamb.degC"].append(tamb) + flux_data["ts_top.degC"].append(ts) + flux_data["hc_top.W.m-2.K-1"].append(convection_function.results.convection_coefficient) + # flux_data["hc_top.W.m-2.K-1"].append(qc/(abs(ts-tamb))) + # flux_data["top_t"].append(heat_flux(region=internal_region)) + # solid.results.statevars.data.tolist() # current temperature (mesh-points) + # region_convection.mask.data.tolist() # surface points True(!) + +n_steps = 20 +time = fem.thermal.TimeStep([solid]) +table = fem.math.linsteps([0, 1], num=n_steps) +air_temperature = fem.math.linsteps([15, 25], num=n_steps) + +ramp = { + time: 18000 * table, # five hours + convection_function["temperature"]: air_temperature, +} + +step = fem.Step( + items=[time, solid, convection_function], ramp=ramp, boundaries=boundaries +) + +flux_data = {"tamb.degC": [], "ts_top.degC": [], "hc_top.W.m-2.K-1": [], "top.W.m-2": []} + +job = fem.Job(steps=[step], callback=callback, flux_data=flux_data).evaluate( + verbose=False +) + +# mesh.view( +# point_data={"Temperature 2 in °C": temperature.values} +# ).plot("Temperature 2 in °C").show() From 25136f8300d536a1415037516bd44ccddd47bf63 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Thu, 7 May 2026 17:26:04 +0200 Subject: [PATCH 35/71] Slight clean-up. Adapt naming. --- src/felupe/thermal/_sbsc_tester.py | 35 +++++++++++++++--------------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 7c50f100c..31523c420 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -87,13 +87,8 @@ def hc_fun(ts, tamb): InputHumidAir.relative_humidity(rh/rhf), ) - l = 1.25 # slab 1 x 1 m^2 + l = 1.25 # assume slab size 5 x 5 m^2 for h_c ra = rayleigh(ts, tamb, l) - # alpha = 2.25E-05 # m^2/s, air at 300 K - # lam_air = 0.0263 # W/(m K), air at 300 K - # pr = 0.707 # air at 300 K - # t_m = 0.5*(ts + tamb) + 273.15 # K - # ra = (9.81 / t_m * abs(ts - tamb) * math.pow(l, 3.0))/alpha/1.59E-7 if ts > tamb: nu = nusselt_horizontal(ra, air.prandtl, hflux='z+') else: @@ -101,6 +96,7 @@ def hc_fun(ts, tamb): return(nu*air.conductivity/l) + mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells region = fem.RegionQuad(mesh) temperature = fem.Field(region, dim=1, values=30.0) @@ -128,21 +124,21 @@ def hc_fun(ts, tamb): temperature=20.0, # °C ) -def callback(stepnumber, substepnumber, substep, flux_data): - """Save mean surface heat flux at top (convective) boundary. +def callback(stepnumber, substepnumber, substep, tstep_data): + """Save time step data at top (convective) boundary. """ tamb = list(ramp.values())[1][substepnumber] ts = list(convection_function.field.extract(grad=False)[0])[0][0:1][0].mean() heat_flux = solid.heat_flux_boundary qc = heat_flux(region=region_convection) - flux_data["top.W.m-2"].append(qc) - flux_data["tamb.degC"].append(tamb) - flux_data["ts_top.degC"].append(ts) - flux_data["hc_top.W.m-2.K-1"].append(convection_function.results.convection_coefficient) - # flux_data["hc_top.W.m-2.K-1"].append(qc/(abs(ts-tamb))) - # flux_data["top_t"].append(heat_flux(region=internal_region)) - # solid.results.statevars.data.tolist() # current temperature (mesh-points) - # region_convection.mask.data.tolist() # surface points True(!) + tstep_data["tstep.s"].append(list(ramp.values())[0][substepnumber]) + tstep_data["qc_top.W.m-2"].append(qc) + tstep_data["tamb.degC"].append(tamb) + tstep_data["ts_top.degC"].append(ts) + tstep_data["hc_top.W.m-2.K-1"].append(convection_function.results.convection_coefficient) + tstep_data["hc_fun_top.W.m-2.K-1"].append(hc_fun(ts, tamb)) + tstep_data["hc_calc_top.W.m-2.K-1"].append(qc/(abs(ts-tamb))) + n_steps = 20 time = fem.thermal.TimeStep([solid]) @@ -158,9 +154,12 @@ def callback(stepnumber, substepnumber, substep, flux_data): items=[time, solid, convection_function], ramp=ramp, boundaries=boundaries ) -flux_data = {"tamb.degC": [], "ts_top.degC": [], "hc_top.W.m-2.K-1": [], "top.W.m-2": []} +tstep_data = {"tstep.s": [], "tamb.degC": [], "ts_top.degC": [], + "hc_top.W.m-2.K-1": [], "qc_top.W.m-2": [], + "hc_fun_top.W.m-2.K-1": [], + "hc_calc_top.W.m-2.K-1": []} -job = fem.Job(steps=[step], callback=callback, flux_data=flux_data).evaluate( +job = fem.Job(steps=[step], callback=callback, tstep_data=tstep_data).evaluate( verbose=False ) From 9d96d3df9f9d000bc950f0c24034cc0d3d67f719 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Fri, 8 May 2026 08:03:49 +0200 Subject: [PATCH 36/71] Undo change to .gitignore on main. --- .gitignore | 1 - 1 file changed, 1 deletion(-) diff --git a/.gitignore b/.gitignore index 7bc3083a0..80af874f5 100644 --- a/.gitignore +++ b/.gitignore @@ -122,7 +122,6 @@ examples/result.xdmf docs/examples/ docs/tutorial/* !docs/tutorial/examples -docs/mesh.png sg_execution_times.rst result.h5 result.xdmf From a1236bd0b74e724ef6b27f2bbd7b1b8142676e30 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Fri, 8 May 2026 16:06:28 +0200 Subject: [PATCH 37/71] Fix h_c calculation from q_c/Dt. Remove obsolete hc_fun(). --- src/felupe/thermal/_sbsc_tester.py | 2 +- .../thermal/_solidbody_surface_convection.py | 16 ---------------- 2 files changed, 1 insertion(+), 17 deletions(-) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 31523c420..6d2b84a42 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -137,7 +137,7 @@ def callback(stepnumber, substepnumber, substep, tstep_data): tstep_data["ts_top.degC"].append(ts) tstep_data["hc_top.W.m-2.K-1"].append(convection_function.results.convection_coefficient) tstep_data["hc_fun_top.W.m-2.K-1"].append(hc_fun(ts, tamb)) - tstep_data["hc_calc_top.W.m-2.K-1"].append(qc/(abs(ts-tamb))) + tstep_data["hc_calc_top.W.m-2.K-1"].append(abs(qc/(ts-tamb))) n_steps = 20 diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 5d721dcfe..1e583ba1e 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -293,19 +293,3 @@ def _matrix(self, field=None, **kwargs): ).assemble(**kwargs) return self.results.stiffness - -# import math - -# def hc_fun(ts, tamb): -# l = 0.25 # slab 1 x 1 m^2 -# alpha = 2.25E-05 # m^2/s, air at 300 K -# pr = 0.707 # air at 300 K -# t_m = 0.5*(ts + tamb) + 273.15 # K -# ra = (9.81 / t_m * abs(ts - tamb) * l**3)/alpha/1.59E-5 -# if (10**4 <= ra <= 10**7) and pr > 0.7: -# nu = 0.54 * math.pow(ra, 0.25) -# elif (10**7 < ra <= 10**11): -# nu = 0.15 * math.pow(ra, 0.33) -# else: -# nu = 1 -# return(nu*0.0263/l) From c03fa0b1d983213dc19a10fdfbf1b83bedcb2935 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 17 May 2026 17:13:52 +0200 Subject: [PATCH 38/71] Test parameter updating. No success. --- src/felupe/thermal/_sbsc_tester.py | 1 + src/felupe/thermal/_solidbody_surface_convection.py | 7 ++++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 6d2b84a42..2be066555 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -148,6 +148,7 @@ def callback(stepnumber, substepnumber, substep, tstep_data): ramp = { time: 18000 * table, # five hours convection_function["temperature"]: air_temperature, + convection_function["convection_coefficient"]: air_temperature, } step = fem.Step( diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 1e583ba1e..3fcecfcb2 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -235,17 +235,18 @@ def __getitem__(self, key): def update(self, temperature): self._update_temperature(temperature) - self._update_convection_coefficient() # adapt hc using cur. temp. + self._update_convection_coefficient(temperature) # adapt hc using cur. temp. def _update_temperature(self, temperature): self.results.temperature = temperature - def _update_convection_coefficient(self): + def _update_convection_coefficient(self, temperature): if callable(self.convection_coefficient): self.results.convection_coefficient =\ self.convection_coefficient( self.field.extract(grad=False)[0], # ts - self.results.temperature # tamb + # self.results.temperature # tamb + temperature # tamb ) else: self.results.convection_coefficient = self.convection_coefficient From b3fe146a2df64c6220b931e52279c04978316a03 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Mon, 18 May 2026 21:05:50 +0200 Subject: [PATCH 39/71] Tester now basically works. Embedded example not checked, make html gives unclear error. Code save. --- src/felupe/thermal/_sbsc_tester.py | 14 +++++++++----- .../thermal/_solidbody_surface_convection.py | 15 +++++++-------- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 2be066555..a694bb8e4 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -3,6 +3,7 @@ Tester version for _solid_body_surface_convection.py example case. """ import math +import numpy as np from pyfluids import HumidAir, InputHumidAir from scipy.constants import g @@ -24,6 +25,7 @@ def pyfluids_units(): return dt_, rh_ +@np.vectorize def rayleigh(ts_c, ti_c, length_, rh=10): """ Calculate dimensionless Rayleigh number Ra. @@ -49,6 +51,7 @@ def rayleigh(ts_c, ti_c, length_, rh=10): return ra +@np.vectorize def nusselt_horizontal(ra, pr, hflux='z+'): """ Calculate dimensionless Nusselt number Nu for horizontal plates for various @@ -73,6 +76,7 @@ def nusselt_horizontal(ra, pr, hflux='z+'): return nu +@np.vectorize def hc_fun(ts, tamb): """ Calculate convection coefficient for horizontal plate. @@ -120,7 +124,7 @@ def hc_fun(ts, tamb): convection_function = fem.thermal.SolidBodySurfaceConvection( field=field_convection, - convection_coefficient=hc_fun(30, 20), + convection_coefficient=hc_fun, #(30, 20), temperature=20.0, # °C ) @@ -135,8 +139,9 @@ def callback(stepnumber, substepnumber, substep, tstep_data): tstep_data["qc_top.W.m-2"].append(qc) tstep_data["tamb.degC"].append(tamb) tstep_data["ts_top.degC"].append(ts) - tstep_data["hc_top.W.m-2.K-1"].append(convection_function.results.convection_coefficient) - tstep_data["hc_fun_top.W.m-2.K-1"].append(hc_fun(ts, tamb)) + tstep_data["hc_top.W.m-2.K-1"].append( + convection_function.results.convection_coefficient.mean()) + tstep_data["hc_fun_top.W.m-2.K-1"].append(hc_fun(ts, tamb).mean()) tstep_data["hc_calc_top.W.m-2.K-1"].append(abs(qc/(ts-tamb))) @@ -147,8 +152,7 @@ def callback(stepnumber, substepnumber, substep, tstep_data): ramp = { time: 18000 * table, # five hours - convection_function["temperature"]: air_temperature, - convection_function["convection_coefficient"]: air_temperature, + convection_function: air_temperature, } step = fem.Step( diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 3fcecfcb2..35b939864 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -151,7 +151,7 @@ class SolidBodySurfaceConvection: >>> air_temperature = fem.math.linsteps([15, 25], num=10) >>> ramp = { ... time: 18000 * table, # five hours - ... convection_constant["temperature"]: air_temperature, + ... convection_constant: air_temperature, ... } >>> step = fem.Step( ... items=[time, solid, convection_constant], ramp=ramp, boundaries=boundaries @@ -170,7 +170,7 @@ class SolidBodySurfaceConvection: >>> convection_function = fem.thermal.SolidBodySurfaceConvection( ... field=field_convection, - ... convection_coefficient=hc_fun(30, 20), + ... convection_coefficient=hc_fun, ... temperature=20.0, # °C ... ) >>> time = fem.thermal.TimeStep([solid]) @@ -178,9 +178,9 @@ class SolidBodySurfaceConvection: >>> air_temperature = fem.math.linsteps([15, 25], num=10) >>> ramp = { ... time: 18000 * table, # five hours - ... convection_function["temperature"]: air_temperature, + ... convection_function: air_temperature, ... } - + ... and run. .. pyvista-plot:: @@ -235,18 +235,17 @@ def __getitem__(self, key): def update(self, temperature): self._update_temperature(temperature) - self._update_convection_coefficient(temperature) # adapt hc using cur. temp. + self._update_convection_coefficient() # adapt hc using cur. temp. def _update_temperature(self, temperature): self.results.temperature = temperature - def _update_convection_coefficient(self, temperature): + def _update_convection_coefficient(self): if callable(self.convection_coefficient): self.results.convection_coefficient =\ self.convection_coefficient( self.field.extract(grad=False)[0], # ts - # self.results.temperature # tamb - temperature # tamb + self.results.temperature # tamb ) else: self.results.convection_coefficient = self.convection_coefficient From 427ff03eb5aab7b45554ef0b36b3524af87bcfb1 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Tue, 19 May 2026 13:43:02 +0200 Subject: [PATCH 40/71] Pytest now returns 'passed'. --- .../thermal/_solidbody_surface_convection.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 35b939864..14c97fbcc 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -85,10 +85,12 @@ class SolidBodySurfaceConvection: .. pyvista-plot:: :context: + :force_static: >>> import math + >>> import numpy as np >>> - >>> def hc_fun(ts, tamb): + >>> def _hc_fun(ts, tamb): ... l = 0.25 # slab 1 x 1 m^2 ... alpha = 2.25E-05 # m^2/s, air at 300 K ... lam_air = 0.0263 # W/(m K), air at 300 K @@ -102,11 +104,13 @@ class SolidBodySurfaceConvection: ... else: ... nu = 0.15 * math.pow(1E11, 0.33) ... return(nu*lam_air/l) + >>> hc_fun = np.vectorize(_hc_fun) Set up the model (a horizontal slab with dimensions 1 x 1 m^2, 0.25 m thick). - + .. pyvista-plot:: :context: + :force_static: >>> import felupe as fem >>> import numpy as np @@ -160,13 +164,14 @@ class SolidBodySurfaceConvection: >>> >>> mesh.view( ... point_data={"Temperature in °C": temperature.values} - ... ).plot("Temperature in °C").show() + ... ).plot("Temperature in °C", off_screen=True).show() And now set up convection to use the function for :math:`h_c` defined above using the same air temperature boundary conditions ... .. pyvista-plot:: :context: + :force_static: >>> convection_function = fem.thermal.SolidBodySurfaceConvection( ... field=field_convection, @@ -193,8 +198,8 @@ class SolidBodySurfaceConvection: >>> job = fem.Job(steps=[step]).evaluate(verbose=False) >>> >>> mesh.view( - ... point_data={"Temperature 2 in °C": temperature.values} - ... ).plot("Temperature 2 in °C").show() + ... point_data={"Temperature hc_fun() in °C": temperature.values} + ... ).plot("Temperature hc_fun() in °C", off_screen=True).show() References ---------- From b1e399ca5c1c564081c866ce11977a0312f3bfe1 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 20 May 2026 13:14:11 +0200 Subject: [PATCH 41/71] Add code outline for line-graph of h_c, t_s, t_i. --- src/felupe/thermal/_sbsc_tester.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index a694bb8e4..702de78a6 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -171,3 +171,12 @@ def callback(stepnumber, substepnumber, substep, tstep_data): # mesh.view( # point_data={"Temperature 2 in °C": temperature.values} # ).plot("Temperature 2 in °C").show() + +# Plot h_c / surface temp. / air temp. vs. time. +# https://felupe.readthedocs.io/en/latest/examples/ex21_nonlinear-truss-analysis.html +# fig, ax = plt.subplots() +# ax.plot(*X[:, [0, -1]].T, ".-", label="Point 3") +# ax.plot(*X[:, [3, -1]].T, ".-", label="Point 4") +# ax.set_xlabel("Displacement X") +# ax.set_ylabel("LPF") +# ax.legend() From 2edd4be00262357e261a33ac86bb7ab9d044bf93 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 20 May 2026 17:19:30 +0200 Subject: [PATCH 42/71] Switch plot code to matplotlib 'parasite simple' approach. Works as expected, to be extended. --- src/felupe/thermal/_sbsc_tester.py | 35 ++++++++++++++++++++++++------ 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 702de78a6..28e5b3530 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -173,10 +173,31 @@ def callback(stepnumber, substepnumber, substep, tstep_data): # ).plot("Temperature 2 in °C").show() # Plot h_c / surface temp. / air temp. vs. time. -# https://felupe.readthedocs.io/en/latest/examples/ex21_nonlinear-truss-analysis.html -# fig, ax = plt.subplots() -# ax.plot(*X[:, [0, -1]].T, ".-", label="Point 3") -# ax.plot(*X[:, [3, -1]].T, ".-", label="Point 4") -# ax.set_xlabel("Displacement X") -# ax.set_ylabel("LPF") -# ax.legend() +# https://matplotlib.org/stable/gallery/axes_grid1/parasite_simple.html +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import host_subplot + +host = host_subplot(111) +par = host.twinx() + +host.set_xlabel("Time (s)") +host.set_ylabel("Convection coefficient (W/(m2 K))") +par.set_ylabel("Temperature (°C)") + +p1, = host.plot(tstep_data["tstep.s"], + tstep_data["hc_top.W.m-2.K-1"], + label="hc") +p2, = par.plot(tstep_data["tstep.s"], + tstep_data["tamb.degC"], + label="t_amb") +p3, = par.plot(tstep_data["tstep.s"], + tstep_data["ts_top.degC"], + label="ts_top") +p3.set_color(p2.get_color()) + +host.legend(labelcolor="linecolor") + +host.yaxis.label.set_color(p1.get_color()) +par.yaxis.label.set_color(p2.get_color()) + +plt.savefig("_test.png") From 17b0a611d48628a4591782edda7e43a933b8fffc Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 20 May 2026 17:24:39 +0200 Subject: [PATCH 43/71] Add alternative hc calculation results. --- src/felupe/thermal/_sbsc_tester.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 28e5b3530..3db4f0ec1 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -184,20 +184,21 @@ def callback(stepnumber, substepnumber, substep, tstep_data): host.set_ylabel("Convection coefficient (W/(m2 K))") par.set_ylabel("Temperature (°C)") -p1, = host.plot(tstep_data["tstep.s"], - tstep_data["hc_top.W.m-2.K-1"], +p1, = host.plot(tstep_data["tstep.s"], tstep_data["hc_top.W.m-2.K-1"], label="hc") -p2, = par.plot(tstep_data["tstep.s"], - tstep_data["tamb.degC"], +p2, = host.plot(tstep_data["tstep.s"], tstep_data["hc_fun_top.W.m-2.K-1"], + label="hc_fun(ts - tamb)") +p3, = host.plot(tstep_data["tstep.s"], tstep_data["hc_calc_top.W.m-2.K-1"], + label="hc from q") +p4, = par.plot(tstep_data["tstep.s"], tstep_data["tamb.degC"], label="t_amb") -p3, = par.plot(tstep_data["tstep.s"], - tstep_data["ts_top.degC"], +p5, = par.plot(tstep_data["tstep.s"], tstep_data["ts_top.degC"], label="ts_top") -p3.set_color(p2.get_color()) +# p3.set_color(p2.get_color()) host.legend(labelcolor="linecolor") -host.yaxis.label.set_color(p1.get_color()) -par.yaxis.label.set_color(p2.get_color()) +# host.yaxis.label.set_color(p1.get_color()) +# par.yaxis.label.set_color(p2.get_color()) plt.savefig("_test.png") From bacac5608e0a391b494c8049f6d0c2b3e0d39db9 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 31 May 2026 11:44:53 +0200 Subject: [PATCH 44/71] Add more complex thermal example. --- examples/ex23_solid_body_thermal-tabs.py | 181 +++++++++++++++++++++++ 1 file changed, 181 insertions(+) create mode 100644 examples/ex23_solid_body_thermal-tabs.py diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py new file mode 100644 index 000000000..0a90d47a0 --- /dev/null +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -0,0 +1,181 @@ +r""" +Thermal Analysis +---------------- + +.. topic:: Analysis of a thermally activated slab setup. + + * use :class:`~felupe.thermal.SolidBodyThermal`, + :class:`~felupe.thermal.SolidBodyHeatFlux`, + :class:`~felupe.thermal.SolidBodySurfaceRadiation` and + :class:`~felupe.thermal.SolidBodySurfaceConvection` + + * evaluate the surface heat flux at internal and external boundaries + with a job :class:`~felupe.Plugin` + + * view the temperature field + + +This example describes a thermally activated concrete slab using a simplified +model and geometry. The model is two-dimensional. The system is set up with two +:class:`solids `. The temperature boundary +conditions include the floor temperature (constant), the ceiling temperature +(constant) and the room air temperature with a :math:`\pm 2` K sinusoidal +variation around its average value with a period of 24 h. + +The heat injection via the pipe layer is constant at 295 W/m2 and directly +injected at the internal concrete surfaces (no pipe material is modelled). + +Surface heat transfer is modelled separately for convection and radiation. +""" +import matplotlib.pyplot as plt +import numpy as np + +import felupe as fem + +# %% +# Define material properties as lists for (reinforced) concrete, plastic and +# insulation. This includes mass density, specific heat capacity and thermal +# conductivity. +density = [2100, 700, 20] # kg/m^3 +specific_heat = [1000, 1000, 1450] # J/(kg K) +thermal_conductivity = [2.1, 0.3, 0.035] # W/(m K) + +# %% +# Set up one mesh per material. If a material consists of multiple areas, these +# are collected in a :class:`mesh container ` and are +# merged into one mesh per material. These meshes per material are then added +# to a mesh container for the construction. + +# << copy ex22 approach >> + +container.plot( + colors=["lightgrey", "sepia"], + labels=["Concrete", "Insulation"], + show_edges=False, +).show() + +# %% +# A top-level temperature field is defined on the whole construction with an initial +# temperature value of 10 °C, and separate fields are defined for each material. The +# surface heat transfer coefficients and ambient temperatures are defined for the +# internal and external boundaries. Thermal solid bodies are created for each material. +regions = [fem.RegionQuad(m) for m in container] +fields = [fem.Field(r, dim=1).as_container() for r in regions] + +# top level temperature field +mesh = container.stack() +region = fem.RegionQuad(mesh) +temperature = fem.Field(region, dim=1, values=10.0) # initial temperature 10 °C +field = fem.FieldContainer([temperature]) + +external_region = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.min()) +external_temperature = fem.Field(external_region, dim=1) +external_field = fem.FieldContainer([external_temperature]) + +internal_region = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.max()) +internal_temperature = fem.Field(internal_region, dim=1) +internal_field = fem.FieldContainer([internal_temperature]) + +external_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + field=external_field, + coefficient=25.0, # W/(m^2 K) + temperature=0.0, # °C +) +internal_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + field=internal_field, + coefficient=7.69, # W/(m^2 K) + temperature=20.0, # °C +) + +materials = [] +for mfield, rho, cp, k in zip(fields, density, specific_heat, thermal_conductivity): + materials.append( + fem.thermal.SolidBodyThermal( + field=mfield, + mass_density=rho, + specific_heat_capacity=cp, + thermal_conductivity=k, + ) + ) + + +# %% +# A callback-function records the mean surface heat flux at the internal and external +# boundaries after each completed time step. The mean surface heat flux is calculated +# by the :meth:`~felupe.thermal.SolidBodyThermal.heat_flux_boundary` method of the +# thermal solid body, which returns the integrated surface heat flux for a given +# boundary region and time step. The mean surface heat flux is stored in the +# ``flux_data`` dictionary, which is passed to the callback function as an argument. +def callback(stepnumber, substepnumber, substep, flux_data): + """Save mean surface heat flux at internal and external boundaries.""" + + heat_flux = materials[0].heat_flux_boundary + flux_data["external"].append(heat_flux(region=external_region)) + flux_data["internal"].append(heat_flux(region=internal_region)) + + +time_steps = fem.math.linsteps([0, 24 * 3600], num=int(24 * 3600 / 720))[1:] + +t_ext = 0 + 1 * np.sin(2 * np.pi * time_steps / 86400) +t_int = 20 + 1 * np.sin(2 * np.pi * time_steps / 86400) + + +# %% +# The time step item is created with the thermal solid bodies. It must be located as the +# first item in the step to properly update the time step in the materials. The internal +# and external heat transfer item values are defined in the ramp, which specifies how +# their values change over time. Finally, a job is created with the step and the +# callback function, and evaluated with the top-level temperature field. A result file +# is created for visualization in Paraview, and the temperature field is saved as point- +# data in the result file. +time = fem.thermal.TimeStep( + [*materials, external_heat_transfer, internal_heat_transfer] +) +ramp = { + time: time_steps, + internal_heat_transfer: t_int, + external_heat_transfer: t_ext, +} +step = fem.Step( + items=[time, *materials, internal_heat_transfer, external_heat_transfer], + ramp=ramp, +) + +flux_data = {"external": [], "internal": []} + +job = fem.Job(steps=[step], callback=callback, flux_data=flux_data).evaluate( + x0=field, + filename="result.xdmf", # create a result file for Paraview + point_data={"Temperature": lambda field, substep: temperature.values}, + point_data_default=False, + cell_data_default=False, +) + +# %% +# Internal and external surface heat flux values are plotted over time. +# +# .. note:: +# +# The heat flux is **positive** when **heat leaves the construction** (here, on the +# external surface), and **negative** when **heat enters the construction** (here, on +# the internal surface). +# +fig, ax = plt.subplots() +ax.plot(time_steps / 3600, flux_data["external"], color="C3", label="external") +ax.plot(time_steps / 3600, flux_data["internal"], color="C0", label="internal") + +tmin, tmax = ax.get_xlim() +ax.plot([tmin, tmax], np.zeros(2), "black", lw=0.5) + +text_kwargs = dict(transform=ax.transAxes, ha="center", va="center") +ax.text(0.5, 0.97, "heat leaves construction", **text_kwargs) +ax.text(0.5, 0.03, "heat enters construction", **text_kwargs) + +ax.legend() +ax.set(xlim=(tmin, tmax), xlabel="time in h", ylabel=r"surface heat flux in W/m$^2$") + + +# %% +# A view on the temperature field at the end of the simulation period visualizes the +# temperature distribution. +field.plot("Field", scalar_bar_vertical=True).show() From 9f52c0fa216115608bd2248cf51bd509ca50f940 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 31 May 2026 14:39:47 +0200 Subject: [PATCH 45/71] Generate simplified geometry. Looks as expected. --- examples/ex23_solid_body_thermal-tabs.py | 34 ++++++++++++++++++++++-- 1 file changed, 32 insertions(+), 2 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index 0a90d47a0..62f9f7b07 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -45,8 +45,38 @@ # are collected in a :class:`mesh container ` and are # merged into one mesh per material. These meshes per material are then added # to a mesh container for the construction. - -# << copy ex22 approach >> +concrete_1a = fem.Rectangle(a=(0.02, 0.0), b=(0.20, 0.22), n=(11, 16)) # left / right +concrete_1b = fem.Rectangle(a=(0.20, 0.0), b=(0.22, 0.10), n=(5, 8)) # pipe bottom / top +concrete_1 = fem.MeshContainer( + [ + concrete_1a, # left + concrete_1b, # pipe 1, bottom + concrete_1b.translate(0.12, axis=1), # pipe1, top + ], + merge=True, +).stack() + +concrete = fem.MeshContainer( + [ + concrete_1, # left + concrete_1.translate(0.2, axis=0), # + concrete_1.translate(0.4, axis=0), # + concrete_1.translate(0.6, axis=0), # + concrete_1a.translate(0.8, axis=0), # right + ], + merge=True, +).stack() + +insulation_1 = fem.Rectangle(a=(0.0, 0.0), b=(0.02, 0.22), n=(4, 16)) # left / right +insulation = fem.MeshContainer( + [ + insulation_1, + insulation_1.translate(1.0, axis=0), + ], + merge=True, +).stack() + +container = fem.MeshContainer([concrete, insulation], merge=True) container.plot( colors=["lightgrey", "sepia"], From 00a814b91f73ddeda2638241b17771564714fe89 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 31 May 2026 17:37:23 +0200 Subject: [PATCH 46/71] Add basic fem model (regions and fields). --- examples/ex23_solid_body_thermal-tabs.py | 63 +++++++++++++++--------- 1 file changed, 41 insertions(+), 22 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index 62f9f7b07..e2be8f0ea 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -95,28 +95,9 @@ # top level temperature field mesh = container.stack() region = fem.RegionQuad(mesh) -temperature = fem.Field(region, dim=1, values=10.0) # initial temperature 10 °C +temperature = fem.Field(region, dim=1, values=20.0) # initial temperature 10 °C field = fem.FieldContainer([temperature]) -external_region = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.min()) -external_temperature = fem.Field(external_region, dim=1) -external_field = fem.FieldContainer([external_temperature]) - -internal_region = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.max()) -internal_temperature = fem.Field(internal_region, dim=1) -internal_field = fem.FieldContainer([internal_temperature]) - -external_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( - field=external_field, - coefficient=25.0, # W/(m^2 K) - temperature=0.0, # °C -) -internal_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( - field=internal_field, - coefficient=7.69, # W/(m^2 K) - temperature=20.0, # °C -) - materials = [] for mfield, rho, cp, k in zip(fields, density, specific_heat, thermal_conductivity): materials.append( @@ -128,6 +109,44 @@ ) ) +# %% +# The surface heat transfer coefficients and ambient temperatures are defined +# for the top and bottom surfaces. +bottom_region = fem.RegionQuadBoundary(mesh, mask=mesh.y == mesh.y.min()) +bottom_temperature = fem.Field(external_region, dim=1) +bottom_field = fem.FieldContainer([external_temperature]) + +top_region = fem.RegionQuadBoundary(mesh, mask=mesh.y == mesh.y.max()) +top_temperature = fem.Field(internal_region, dim=1) +top_field = fem.FieldContainer([internal_temperature]) + +top_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + field=external_field, + coefficient=7.69, # W/(m^2 K) + temperature=20.0, # °C +) +bottom_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + field=internal_field, + coefficient=5.0, # W/(m^2 K) + temperature=20.0, # °C +) + +# %% +# Heat flux on pipe walls is defined. +center_points = np.array([[0.2, 0.1], [0.4, 0.], [0.6, 0.1], [0.8, 0.1]]) + +pipe_region = [] +pipe_field = [] +pipe_flux = [] +for idx, p in enumerate(center_points): + mask = (mesh.points[:, None, :] == p[None, :, :]).all(axis=2).any(axis=1) + pipe_region.append(fem.RegionQuadBoundary(mesh, mask=mask)) + pipe_field.append(fem.FieldContainer([fem.Field(pipe_region[idx], dim=1)])) + pipe_flux.append(fem.thermal.SolidBodyHeatFlux( + field=pipe_field[idx], + heat_flux=-294.6, # W / m^2, 74/(4*3.14*2*ri) + )) + # %% # A callback-function records the mean surface heat flux at the internal and external @@ -140,8 +159,8 @@ def callback(stepnumber, substepnumber, substep, flux_data): """Save mean surface heat flux at internal and external boundaries.""" heat_flux = materials[0].heat_flux_boundary - flux_data["external"].append(heat_flux(region=external_region)) - flux_data["internal"].append(heat_flux(region=internal_region)) + flux_data["top"].append(heat_flux(region=top_region)) + flux_data["bottom"].append(heat_flux(region=bottom_region)) time_steps = fem.math.linsteps([0, 24 * 3600], num=int(24 * 3600 / 720))[1:] From 534dbcd087b827c40092ff24e4d92dae518787be Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 31 May 2026 19:06:29 +0200 Subject: [PATCH 47/71] Add basic code to set heat flux at 'pipe' surfaces. Actual center points used do not seem to be correct. --- examples/ex23_solid_body_thermal-tabs.py | 29 ++++++++++++------------ 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index e2be8f0ea..61f4ff282 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -78,11 +78,11 @@ container = fem.MeshContainer([concrete, insulation], merge=True) -container.plot( - colors=["lightgrey", "sepia"], - labels=["Concrete", "Insulation"], - show_edges=False, -).show() +# container.plot( +# colors=["lightgrey", "sepia"], +# labels=["Concrete", "Insulation"], +# show_edges=False, +# ).show() # %% # A top-level temperature field is defined on the whole construction with an initial @@ -113,41 +113,40 @@ # The surface heat transfer coefficients and ambient temperatures are defined # for the top and bottom surfaces. bottom_region = fem.RegionQuadBoundary(mesh, mask=mesh.y == mesh.y.min()) -bottom_temperature = fem.Field(external_region, dim=1) -bottom_field = fem.FieldContainer([external_temperature]) +bottom_temperature = fem.Field(bottom_region, dim=1) +bottom_field = fem.FieldContainer([bottom_temperature]) top_region = fem.RegionQuadBoundary(mesh, mask=mesh.y == mesh.y.max()) -top_temperature = fem.Field(internal_region, dim=1) -top_field = fem.FieldContainer([internal_temperature]) +top_temperature = fem.Field(top_region, dim=1) +top_field = fem.FieldContainer([top_temperature]) top_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( - field=external_field, + field=top_field, coefficient=7.69, # W/(m^2 K) temperature=20.0, # °C ) bottom_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( - field=internal_field, + field=bottom_field, coefficient=5.0, # W/(m^2 K) temperature=20.0, # °C ) # %% # Heat flux on pipe walls is defined. -center_points = np.array([[0.2, 0.1], [0.4, 0.], [0.6, 0.1], [0.8, 0.1]]) +center_points = np.asarray([[0.2, 0.1], [0.4, 0.1], [0.6, 0.1], [0.8, 0.1]]) pipe_region = [] pipe_field = [] pipe_flux = [] for idx, p in enumerate(center_points): - mask = (mesh.points[:, None, :] == p[None, :, :]).all(axis=2).any(axis=1) + mask = np.isclose(mesh.points[:, None, :], p[:], atol=0.015).all(axis=2).any(axis=1) pipe_region.append(fem.RegionQuadBoundary(mesh, mask=mask)) pipe_field.append(fem.FieldContainer([fem.Field(pipe_region[idx], dim=1)])) pipe_flux.append(fem.thermal.SolidBodyHeatFlux( field=pipe_field[idx], - heat_flux=-294.6, # W / m^2, 74/(4*3.14*2*ri) + heat_flux=-231.25, # W / m^2, 74/(4*4*0.02) )) - # %% # A callback-function records the mean surface heat flux at the internal and external # boundaries after each completed time step. The mean surface heat flux is calculated From 81db8aada7e49a939f4f79a7042ae56d8cd3bf51 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 31 May 2026 19:07:26 +0200 Subject: [PATCH 48/71] Code save for tester line plot creation (different approach w/o "esoteric" packages). Not working, yet. --- src/felupe/thermal/_sbsc_tester.py | 49 +++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 15 deletions(-) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 3db4f0ec1..07f7a8c40 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -174,29 +174,48 @@ def callback(stepnumber, substepnumber, substep, tstep_data): # Plot h_c / surface temp. / air temp. vs. time. # https://matplotlib.org/stable/gallery/axes_grid1/parasite_simple.html +# or +# https://stackoverflow.com/questions/9103166/multiple-axis-in-matplotlib-with-different-scales import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import host_subplot -host = host_subplot(111) -par = host.twinx() +# host = host_subplot(111) +# par = host.twinx() +# par2 = host.twinx() -host.set_xlabel("Time (s)") -host.set_ylabel("Convection coefficient (W/(m2 K))") -par.set_ylabel("Temperature (°C)") -p1, = host.plot(tstep_data["tstep.s"], tstep_data["hc_top.W.m-2.K-1"], +fig, ax = plt.subplots() +fig.subplots_adjust(right=0.75) + +twin1 = ax.twinx() +twin2 = ax.twinx() + +# twin2.spines.right.set_position(("axes", 1.2)) + +ax.set_xlabel("Time (s)") +ax.set_ylabel("Convection coefficient (W/(m2 K))") +twin1.set_ylabel("Temperature (°C)") +twin2.set_ylabel("Heat flux (W/m2)") + +p1 = ax.plot(tstep_data["tstep.s"], tstep_data["hc_top.W.m-2.K-1"], label="hc") -p2, = host.plot(tstep_data["tstep.s"], tstep_data["hc_fun_top.W.m-2.K-1"], - label="hc_fun(ts - tamb)") -p3, = host.plot(tstep_data["tstep.s"], tstep_data["hc_calc_top.W.m-2.K-1"], - label="hc from q") -p4, = par.plot(tstep_data["tstep.s"], tstep_data["tamb.degC"], - label="t_amb") -p5, = par.plot(tstep_data["tstep.s"], tstep_data["ts_top.degC"], - label="ts_top") +p2 = ax.plot(tstep_data["tstep.s"], tstep_data["hc_fun_top.W.m-2.K-1"], + label="hc_fun(ts_top - t_amb)") +p3 = twin1.plot(tstep_data["tstep.s"], tstep_data["tamb.degC"], + label="t_amb") +p4 = twin1.plot(tstep_data["tstep.s"], tstep_data["ts_top.degC"], + label="ts_top") +# p3, = host.plot(tstep_data["tstep.s"], tstep_data["hc_calc_top.W.m-2.K-1"], +# label="hc from q") # p3.set_color(p2.get_color()) +p5 = twin2.plot(tstep_data["tstep.s"], tstep_data["qc_top.W.m-2"], + label="qc_top") + +ax.legend(handles=p1+p2+p3+p4+p5, loc='best') + +ax.legend(labelcolor="linecolor") -host.legend(labelcolor="linecolor") +twin2.spines['right'].set_position(('outward', 45)) # host.yaxis.label.set_color(p1.get_color()) # par.yaxis.label.set_color(p2.get_color()) From 072fdac3e8e7e91ca8d4820e1a144cb6dbacb808 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Mon, 1 Jun 2026 18:54:33 +0200 Subject: [PATCH 49/71] Complete modelling. Basically works, however, results show issues (probably with mesh). Refactor mesh generation - not quite working yet, code save. --- examples/ex23_solid_body_thermal-tabs.py | 48 +++++++++++++----------- 1 file changed, 26 insertions(+), 22 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index 61f4ff282..fe6e6bc5c 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -45,13 +45,13 @@ # are collected in a :class:`mesh container ` and are # merged into one mesh per material. These meshes per material are then added # to a mesh container for the construction. -concrete_1a = fem.Rectangle(a=(0.02, 0.0), b=(0.20, 0.22), n=(11, 16)) # left / right -concrete_1b = fem.Rectangle(a=(0.20, 0.0), b=(0.22, 0.10), n=(5, 8)) # pipe bottom / top +concrete_1a = fem.Rectangle(a=(0.0, 0.0), b=(0.18, 0.22), n=(19, 23)) # left / right +concrete_1b = fem.Rectangle(a=(0.0, 0.0), b=(0.02, 0.10), n=(3, 11)) # pipe bottom / top concrete_1 = fem.MeshContainer( [ - concrete_1a, # left - concrete_1b, # pipe 1, bottom - concrete_1b.translate(0.12, axis=1), # pipe1, top + concrete_1a.translate(0.02, axis=0), # left + concrete_1b.translate(0.18, axis=0), # pipe 1, bottom + concrete_1b.translate(0.18, axis=0).translate(0.12, axis=1), # pipe1, top ], merge=True, ).stack() @@ -67,7 +67,7 @@ merge=True, ).stack() -insulation_1 = fem.Rectangle(a=(0.0, 0.0), b=(0.02, 0.22), n=(4, 16)) # left / right +insulation_1 = fem.Rectangle(a=(0.0, 0.0), b=(0.02, 0.22), n=(3, 23)) # left / right insulation = fem.MeshContainer( [ insulation_1, @@ -78,11 +78,11 @@ container = fem.MeshContainer([concrete, insulation], merge=True) -# container.plot( -# colors=["lightgrey", "sepia"], -# labels=["Concrete", "Insulation"], -# show_edges=False, -# ).show() +container.plot( + colors=["lightgrey", "sepia"], + labels=["Concrete", "Insulation"], + show_edges=False, +).show() # %% # A top-level temperature field is defined on the whole construction with an initial @@ -133,13 +133,13 @@ # %% # Heat flux on pipe walls is defined. -center_points = np.asarray([[0.2, 0.1], [0.4, 0.1], [0.6, 0.1], [0.8, 0.1]]) +center_points = np.asarray([[0.21, 0.11], [0.41, 0.11], [0.61, 0.11], [0.81, 0.11]]) pipe_region = [] pipe_field = [] pipe_flux = [] for idx, p in enumerate(center_points): - mask = np.isclose(mesh.points[:, None, :], p[:], atol=0.015).all(axis=2).any(axis=1) + mask = np.isclose(mesh.points[:, None, :], p[:], rtol=0.11, atol=0.015).all(axis=2).any(axis=1) pipe_region.append(fem.RegionQuadBoundary(mesh, mask=mask)) pipe_field.append(fem.FieldContainer([fem.Field(pipe_region[idx], dim=1)])) pipe_flux.append(fem.thermal.SolidBodyHeatFlux( @@ -160,12 +160,16 @@ def callback(stepnumber, substepnumber, substep, flux_data): heat_flux = materials[0].heat_flux_boundary flux_data["top"].append(heat_flux(region=top_region)) flux_data["bottom"].append(heat_flux(region=bottom_region)) + pflux = 0 + for p_ in pipe_region: + pflux += heat_flux(region=p_) + flux_data["pipes"].append(pflux) time_steps = fem.math.linsteps([0, 24 * 3600], num=int(24 * 3600 / 720))[1:] -t_ext = 0 + 1 * np.sin(2 * np.pi * time_steps / 86400) -t_int = 20 + 1 * np.sin(2 * np.pi * time_steps / 86400) +t_ext = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) +t_int = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) # %% @@ -177,19 +181,19 @@ def callback(stepnumber, substepnumber, substep, flux_data): # is created for visualization in Paraview, and the temperature field is saved as point- # data in the result file. time = fem.thermal.TimeStep( - [*materials, external_heat_transfer, internal_heat_transfer] + [*materials, top_heat_transfer, bottom_heat_transfer] ) ramp = { time: time_steps, - internal_heat_transfer: t_int, - external_heat_transfer: t_ext, + top_heat_transfer: t_int, + bottom_heat_transfer: t_ext, } step = fem.Step( - items=[time, *materials, internal_heat_transfer, external_heat_transfer], + items=[time, *materials, top_heat_transfer, bottom_heat_transfer] + pipe_flux, ramp=ramp, ) -flux_data = {"external": [], "internal": []} +flux_data = {"top": [], "bottom": [], "pipes": []} job = fem.Job(steps=[step], callback=callback, flux_data=flux_data).evaluate( x0=field, @@ -209,8 +213,8 @@ def callback(stepnumber, substepnumber, substep, flux_data): # the internal surface). # fig, ax = plt.subplots() -ax.plot(time_steps / 3600, flux_data["external"], color="C3", label="external") -ax.plot(time_steps / 3600, flux_data["internal"], color="C0", label="internal") +ax.plot(time_steps / 3600, flux_data["top"], color="C3", label="top") +ax.plot(time_steps / 3600, flux_data["bottom"], color="C0", label="bottom") tmin, tmax = ax.get_xlim() ax.plot([tmin, tmax], np.zeros(2), "black", lw=0.5) From d178f5150a1103a7240ad970e6156c06b69d803e Mon Sep 17 00:00:00 2001 From: AGeissler Date: Mon, 1 Jun 2026 19:45:30 +0200 Subject: [PATCH 50/71] Fix of geometry creation w/ refactored approach. Issue (probably): boundary mesh identification for 'pipes'. --- examples/ex23_solid_body_thermal-tabs.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index fe6e6bc5c..8b08133fa 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -50,8 +50,8 @@ concrete_1 = fem.MeshContainer( [ concrete_1a.translate(0.02, axis=0), # left - concrete_1b.translate(0.18, axis=0), # pipe 1, bottom - concrete_1b.translate(0.18, axis=0).translate(0.12, axis=1), # pipe1, top + concrete_1b.translate(0.20, axis=0), # pipe 1, bottom + concrete_1b.translate(0.20, axis=0).translate(0.12, axis=1), # pipe1, top ], merge=True, ).stack() @@ -62,7 +62,7 @@ concrete_1.translate(0.2, axis=0), # concrete_1.translate(0.4, axis=0), # concrete_1.translate(0.6, axis=0), # - concrete_1a.translate(0.8, axis=0), # right + concrete_1a.translate(0.82, axis=0), # right ], merge=True, ).stack() @@ -78,11 +78,11 @@ container = fem.MeshContainer([concrete, insulation], merge=True) -container.plot( - colors=["lightgrey", "sepia"], - labels=["Concrete", "Insulation"], - show_edges=False, -).show() +# container.plot( +# colors=["lightgrey", "sepia"], +# labels=["Concrete", "Insulation"], +# show_edges=False, +# ).show() # %% # A top-level temperature field is defined on the whole construction with an initial @@ -139,7 +139,7 @@ pipe_field = [] pipe_flux = [] for idx, p in enumerate(center_points): - mask = np.isclose(mesh.points[:, None, :], p[:], rtol=0.11, atol=0.015).all(axis=2).any(axis=1) + mask = np.isclose(mesh.points[:, None, :], p[:], rtol=0.048, atol=0.0101).all(axis=2).any(axis=1) pipe_region.append(fem.RegionQuadBoundary(mesh, mask=mask)) pipe_field.append(fem.FieldContainer([fem.Field(pipe_region[idx], dim=1)])) pipe_flux.append(fem.thermal.SolidBodyHeatFlux( From 00296fad5617bb09a5e941db9eb888f043d7b48b Mon Sep 17 00:00:00 2001 From: AGeissler Date: Wed, 3 Jun 2026 07:44:32 +0200 Subject: [PATCH 51/71] Add explicit 'square edge point' calculation. Remove unused material properties. Masking still not working correctly. --- examples/ex23_solid_body_thermal-tabs.py | 52 ++++++++++++++++++++---- 1 file changed, 45 insertions(+), 7 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index 8b08133fa..ddf6b910e 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -33,12 +33,11 @@ import felupe as fem # %% -# Define material properties as lists for (reinforced) concrete, plastic and -# insulation. This includes mass density, specific heat capacity and thermal -# conductivity. -density = [2100, 700, 20] # kg/m^3 -specific_heat = [1000, 1000, 1450] # J/(kg K) -thermal_conductivity = [2.1, 0.3, 0.035] # W/(m K) +# Define material properties as lists for (reinforced) concrete and insulation. +# This includes mass density, specific heat capacity and thermal conductivity. +density = [2100, 20] # kg/m^3 +specific_heat = [1000, 1450] # J/(kg K) +thermal_conductivity = [2.1, 0.035] # W/(m K) # %% # Set up one mesh per material. If a material consists of multiple areas, these @@ -135,11 +134,50 @@ # Heat flux on pipe walls is defined. center_points = np.asarray([[0.21, 0.11], [0.41, 0.11], [0.61, 0.11], [0.81, 0.11]]) +# Calculate Square Coordinates +def calculate_square_coordinates(center_x, center_y, size, sections): + step = size / sections + coordinates = [] + + for i in range(sections + 1): + for j in range(sections + 1): + x = round(center_x - (size / 2) + (i * step), 6) + y = round(center_y - (size / 2) + (j * step), 6) + coordinates.append([x, y]) + + return coordinates + +# Calculate Square Edge Points +def calculate_square_edge_points(center_x, center_y, size, sections): + points = [] + step = size / sections + + for i in range(sections + 1): + # Bottom edge + points.append((center_x - size / 2 + i * step, center_y + size / 2)) + # Right edge + points.append((center_x + size / 2, center_y - size / 2 + i * step)) + # Top edge + points.append((center_x + size / 2 - i * step, center_y - size / 2)) + # Left edge + points.append((center_x - size / 2, center_y + size / 2 - i * step)) + + return list(set(points)) + + +coords = [] +for p in center_points: + coords = coords + calculate_square_edge_points(p[0], p[1], 0.02, 2) + +coords = np.asarray(coords) + pipe_region = [] pipe_field = [] pipe_flux = [] -for idx, p in enumerate(center_points): +for idx, p in enumerate(coords): mask = np.isclose(mesh.points[:, None, :], p[:], rtol=0.048, atol=0.0101).all(axis=2).any(axis=1) + # mask = np.isclose(mesh.points[:, None, :], p[:]).all(axis=2).any(axis=1) + # mask = (mesh.points[:, None, :] == p[:]).all(axis=2).any(axis=1) pipe_region.append(fem.RegionQuadBoundary(mesh, mask=mask)) pipe_field.append(fem.FieldContainer([fem.Field(pipe_region[idx], dim=1)])) pipe_flux.append(fem.thermal.SolidBodyHeatFlux( From f36ff14508e1bdadd1881bd43447b18c99b1d00e Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 3 Jun 2026 14:27:49 +0200 Subject: [PATCH 52/71] Code save prior to clean-up. --- examples/ex23_solid_body_thermal-tabs.py | 102 +++++++++++++---------- 1 file changed, 60 insertions(+), 42 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index ddf6b910e..7371519b1 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -53,6 +53,7 @@ concrete_1b.translate(0.20, axis=0).translate(0.12, axis=1), # pipe1, top ], merge=True, + decimals=6, ).stack() concrete = fem.MeshContainer( @@ -64,6 +65,7 @@ concrete_1a.translate(0.82, axis=0), # right ], merge=True, + decimals=6, ).stack() insulation_1 = fem.Rectangle(a=(0.0, 0.0), b=(0.02, 0.22), n=(3, 23)) # left / right @@ -73,9 +75,10 @@ insulation_1.translate(1.0, axis=0), ], merge=True, + decimals=6, ).stack() -container = fem.MeshContainer([concrete, insulation], merge=True) +container = fem.MeshContainer([concrete, insulation], merge=True, decimals=6) # container.plot( # colors=["lightgrey", "sepia"], @@ -134,50 +137,66 @@ # Heat flux on pipe walls is defined. center_points = np.asarray([[0.21, 0.11], [0.41, 0.11], [0.61, 0.11], [0.81, 0.11]]) -# Calculate Square Coordinates -def calculate_square_coordinates(center_x, center_y, size, sections): - step = size / sections - coordinates = [] - - for i in range(sections + 1): - for j in range(sections + 1): - x = round(center_x - (size / 2) + (i * step), 6) - y = round(center_y - (size / 2) + (j * step), 6) - coordinates.append([x, y]) - - return coordinates - -# Calculate Square Edge Points -def calculate_square_edge_points(center_x, center_y, size, sections): - points = [] - step = size / sections - - for i in range(sections + 1): - # Bottom edge - points.append((center_x - size / 2 + i * step, center_y + size / 2)) - # Right edge - points.append((center_x + size / 2, center_y - size / 2 + i * step)) - # Top edge - points.append((center_x + size / 2 - i * step, center_y - size / 2)) - # Left edge - points.append((center_x - size / 2, center_y + size / 2 - i * step)) - - return list(set(points)) - - -coords = [] -for p in center_points: - coords = coords + calculate_square_edge_points(p[0], p[1], 0.02, 2) - -coords = np.asarray(coords) +# Calculate Square Coordinates (code by MS Copilot, three iterations). +# def calculate_square_edge_points(centers, size, sections, tol=1e-12): +# centers = np.atleast_2d(centers) # (N, 2) +# cx = centers[:, 0][:, None] # (N, 1) +# cy = centers[:, 1][:, None] # (N, 1) + +# half = size / 2 +# t = np.linspace(0, size, sections + 1) # (sections+1,) + +# # --- Build edges (broadcasting ensures matching shapes) --- + +# # Bottom edge: left → right +# bottom_x = cx - half + t +# bottom_y = np.full_like(bottom_x, cy - half) +# bottom = np.stack((bottom_x, bottom_y), axis=2) + +# # Right edge: bottom → top (skip first) +# right_y = cy - half + t[1:] +# right_x = np.full_like(right_y, cx + half) +# right = np.stack((right_x, right_y), axis=2) + +# # Top edge: right → left (skip first) +# top_x = cx + half - t[1:] +# top_y = np.full_like(top_x, cy + half) +# top = np.stack((top_x, top_y), axis=2) + +# # Left edge: top → bottom (skip first) +# left_y = cy + half - t[1:] +# left_x = np.full_like(left_y, cx - half) +# left = np.stack((left_x, left_y), axis=2) + +# # --- Flatten each edge to 2‑D --- +# bottom = bottom.reshape(-1, 2) +# right = right.reshape(-1, 2) +# top = top.reshape(-1, 2) +# left = left.reshape(-1, 2) + +# # --- Combine --- +# pts = np.vstack((bottom, right, top, left)) + +# # --- Deduplicate with tolerance --- +# pts_rounded = np.round(pts / tol) * tol +# _, idx = np.unique(pts_rounded, axis=0, return_index=True) + +# return pts[idx] + +# coords = calculate_square_edge_points(center_points, 0.02, 2) pipe_region = [] pipe_field = [] pipe_flux = [] -for idx, p in enumerate(coords): - mask = np.isclose(mesh.points[:, None, :], p[:], rtol=0.048, atol=0.0101).all(axis=2).any(axis=1) - # mask = np.isclose(mesh.points[:, None, :], p[:]).all(axis=2).any(axis=1) - # mask = (mesh.points[:, None, :] == p[:]).all(axis=2).any(axis=1) + +# tolerance = 0.001 # 1 mm if your units are meters +# mask = ( +# np.abs(mesh.points[:, None, :] - coords[None, :, :]) < tolerance +# ).all(axis=2).any(axis=1) + +# for idx, p in enumerate(coords): +for idx, p in enumerate(center_points): + mask = np.isclose(mesh.points[:, None, :], p[:], rtol=0.05, atol=0.0101).all(axis=2).any(axis=1) pipe_region.append(fem.RegionQuadBoundary(mesh, mask=mask)) pipe_field.append(fem.FieldContainer([fem.Field(pipe_region[idx], dim=1)])) pipe_flux.append(fem.thermal.SolidBodyHeatFlux( @@ -203,7 +222,6 @@ def callback(stepnumber, substepnumber, substep, flux_data): pflux += heat_flux(region=p_) flux_data["pipes"].append(pflux) - time_steps = fem.math.linsteps([0, 24 * 3600], num=int(24 * 3600 / 720))[1:] t_ext = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) From 05ef285b28bb00ad2e23a59587fcef01e28ba74f Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 3 Jun 2026 14:51:38 +0200 Subject: [PATCH 53/71] Basically works with simple heat transfer. Extension with side surface heat transfer does not work, reason unclear! --- examples/ex23_solid_body_thermal-tabs.py | 90 +++++++----------------- 1 file changed, 25 insertions(+), 65 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index 7371519b1..1fc4d961a 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -113,7 +113,15 @@ # %% # The surface heat transfer coefficients and ambient temperatures are defined -# for the top and bottom surfaces. +# for the side, top and bottom surfaces. +side_region1 = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.min()) +side_temperature1 = fem.Field(side_region1, dim=1) + +side_region2 = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.max()) +side_temperature2 = fem.Field(side_region2, dim=1) + +side_field = fem.FieldContainer([side_temperature1, side_temperature2]) + bottom_region = fem.RegionQuadBoundary(mesh, mask=mesh.y == mesh.y.min()) bottom_temperature = fem.Field(bottom_region, dim=1) bottom_field = fem.FieldContainer([bottom_temperature]) @@ -122,9 +130,14 @@ top_temperature = fem.Field(top_region, dim=1) top_field = fem.FieldContainer([top_temperature]) +side_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + field=side_field, + coefficient=7.69, # W/(m^2 K) + temperature=20.0, # °C +) top_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( field=top_field, - coefficient=7.69, # W/(m^2 K) + coefficient=10.0, # W/(m^2 K) temperature=20.0, # °C ) bottom_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( @@ -137,64 +150,9 @@ # Heat flux on pipe walls is defined. center_points = np.asarray([[0.21, 0.11], [0.41, 0.11], [0.61, 0.11], [0.81, 0.11]]) -# Calculate Square Coordinates (code by MS Copilot, three iterations). -# def calculate_square_edge_points(centers, size, sections, tol=1e-12): -# centers = np.atleast_2d(centers) # (N, 2) -# cx = centers[:, 0][:, None] # (N, 1) -# cy = centers[:, 1][:, None] # (N, 1) - -# half = size / 2 -# t = np.linspace(0, size, sections + 1) # (sections+1,) - -# # --- Build edges (broadcasting ensures matching shapes) --- - -# # Bottom edge: left → right -# bottom_x = cx - half + t -# bottom_y = np.full_like(bottom_x, cy - half) -# bottom = np.stack((bottom_x, bottom_y), axis=2) - -# # Right edge: bottom → top (skip first) -# right_y = cy - half + t[1:] -# right_x = np.full_like(right_y, cx + half) -# right = np.stack((right_x, right_y), axis=2) - -# # Top edge: right → left (skip first) -# top_x = cx + half - t[1:] -# top_y = np.full_like(top_x, cy + half) -# top = np.stack((top_x, top_y), axis=2) - -# # Left edge: top → bottom (skip first) -# left_y = cy + half - t[1:] -# left_x = np.full_like(left_y, cx - half) -# left = np.stack((left_x, left_y), axis=2) - -# # --- Flatten each edge to 2‑D --- -# bottom = bottom.reshape(-1, 2) -# right = right.reshape(-1, 2) -# top = top.reshape(-1, 2) -# left = left.reshape(-1, 2) - -# # --- Combine --- -# pts = np.vstack((bottom, right, top, left)) - -# # --- Deduplicate with tolerance --- -# pts_rounded = np.round(pts / tol) * tol -# _, idx = np.unique(pts_rounded, axis=0, return_index=True) - -# return pts[idx] - -# coords = calculate_square_edge_points(center_points, 0.02, 2) - pipe_region = [] pipe_field = [] pipe_flux = [] - -# tolerance = 0.001 # 1 mm if your units are meters -# mask = ( -# np.abs(mesh.points[:, None, :] - coords[None, :, :]) < tolerance -# ).all(axis=2).any(axis=1) - -# for idx, p in enumerate(coords): for idx, p in enumerate(center_points): mask = np.isclose(mesh.points[:, None, :], p[:], rtol=0.05, atol=0.0101).all(axis=2).any(axis=1) pipe_region.append(fem.RegionQuadBoundary(mesh, mask=mask)) @@ -224,8 +182,8 @@ def callback(stepnumber, substepnumber, substep, flux_data): time_steps = fem.math.linsteps([0, 24 * 3600], num=int(24 * 3600 / 720))[1:] -t_ext = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) -t_int = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) +t_air = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) +# t_int = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) # %% @@ -236,16 +194,18 @@ def callback(stepnumber, substepnumber, substep, flux_data): # callback function, and evaluated with the top-level temperature field. A result file # is created for visualization in Paraview, and the temperature field is saved as point- # data in the result file. -time = fem.thermal.TimeStep( - [*materials, top_heat_transfer, bottom_heat_transfer] -) +model_list = [*materials, top_heat_transfer, bottom_heat_transfer] +# model_list = [*materials, side_heat_transfer, top_heat_transfer, bottom_heat_transfer] + +time = fem.thermal.TimeStep(model_list) ramp = { time: time_steps, - top_heat_transfer: t_int, - bottom_heat_transfer: t_ext, + side_heat_transfer: t_air, + top_heat_transfer: t_air, + bottom_heat_transfer: t_air, } step = fem.Step( - items=[time, *materials, top_heat_transfer, bottom_heat_transfer] + pipe_flux, + items=[time] + model_list + pipe_flux, ramp=ramp, ) From 6468c7bcb8d34abe970c7d32be07a71b5651ef9f Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 3 Jun 2026 16:16:27 +0200 Subject: [PATCH 54/71] Defining the side heat transfer separately works. Now base "simple heat transfer" can be viewed as O.K. --- examples/ex23_solid_body_thermal-tabs.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index 1fc4d961a..eab31a161 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -116,11 +116,11 @@ # for the side, top and bottom surfaces. side_region1 = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.min()) side_temperature1 = fem.Field(side_region1, dim=1) +side_field1 = fem.FieldContainer([side_temperature1]) side_region2 = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.max()) side_temperature2 = fem.Field(side_region2, dim=1) - -side_field = fem.FieldContainer([side_temperature1, side_temperature2]) +side_field2 = fem.FieldContainer([side_temperature2]) bottom_region = fem.RegionQuadBoundary(mesh, mask=mesh.y == mesh.y.min()) bottom_temperature = fem.Field(bottom_region, dim=1) @@ -130,8 +130,13 @@ top_temperature = fem.Field(top_region, dim=1) top_field = fem.FieldContainer([top_temperature]) -side_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( - field=side_field, +side1_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + field=side_field1, + coefficient=7.69, # W/(m^2 K) + temperature=20.0, # °C +) +side2_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + field=side_field2, coefficient=7.69, # W/(m^2 K) temperature=20.0, # °C ) @@ -194,13 +199,15 @@ def callback(stepnumber, substepnumber, substep, flux_data): # callback function, and evaluated with the top-level temperature field. A result file # is created for visualization in Paraview, and the temperature field is saved as point- # data in the result file. -model_list = [*materials, top_heat_transfer, bottom_heat_transfer] -# model_list = [*materials, side_heat_transfer, top_heat_transfer, bottom_heat_transfer] +# model_list = [*materials, top_heat_transfer, bottom_heat_transfer] +model_list = [*materials, side1_heat_transfer, side2_heat_transfer, + top_heat_transfer, bottom_heat_transfer] time = fem.thermal.TimeStep(model_list) ramp = { time: time_steps, - side_heat_transfer: t_air, + side1_heat_transfer: t_air, + side2_heat_transfer: t_air, top_heat_transfer: t_air, bottom_heat_transfer: t_air, } From 1392fddefc319a034562badda36f8ccf4c65225e Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Wed, 3 Jun 2026 16:51:15 +0200 Subject: [PATCH 55/71] Initial steps to add class 'FreeConvection' as a kind of material (not sure, if best location ...). Not working yet, "not found". --- examples/ex23_solid_body_thermal-tabs.py | 2 +- src/felupe/__init__.py | 2 + src/felupe/constitution/__init__.py | 2 + .../constitution/heat_transfer/__init__.py | 11 ++ .../heat_transfer/_free_convection.py | 117 ++++++++++++++++++ src/felupe/thermal/_sbsc_tester.py | 101 +-------------- 6 files changed, 135 insertions(+), 100 deletions(-) create mode 100644 src/felupe/constitution/heat_transfer/__init__.py create mode 100644 src/felupe/constitution/heat_transfer/_free_convection.py diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index eab31a161..88f9583fb 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -188,7 +188,7 @@ def callback(stepnumber, substepnumber, substep, flux_data): time_steps = fem.math.linsteps([0, 24 * 3600], num=int(24 * 3600 / 720))[1:] t_air = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) -# t_int = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) +# t_sur = 20 + 0.5 * np.sin(2 * np.pi * time_steps / 86400) # %% diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index 1a6597edb..477cb49a4 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -42,6 +42,7 @@ linear_elastic, linear_elastic_plastic_isotropic_hardening, linear_elastic_viscoelastic, + FreeConvection, ) from .dof import Boundary, BoundaryDict from .element import ArbitraryOrderLagrange as ArbitraryOrderLagrangeElement @@ -202,6 +203,7 @@ "VolumeChange", "linear_elastic", "linear_elastic_viscoelastic", + "FreeConvection", "linear_elastic_plastic_isotropic_hardening", "Boundary", "BoundaryDict", diff --git a/src/felupe/constitution/__init__.py b/src/felupe/constitution/__init__.py index 5f6e1af55..a685c559d 100644 --- a/src/felupe/constitution/__init__.py +++ b/src/felupe/constitution/__init__.py @@ -23,6 +23,7 @@ linear_elastic_plastic_isotropic_hardening, linear_elastic_viscoelastic, ) +from .heat_transfer import FreeConvection __all__ = [ "NeoHooke", @@ -55,6 +56,7 @@ "constitutive_material", "CompositeMaterial", "Volumetric", + "FreeConvection", ] try: from .tensortrax import Hyperelastic diff --git a/src/felupe/constitution/heat_transfer/__init__.py b/src/felupe/constitution/heat_transfer/__init__.py new file mode 100644 index 000000000..17b0509f3 --- /dev/null +++ b/src/felupe/constitution/heat_transfer/__init__.py @@ -0,0 +1,11 @@ +""" +constitution.heat_transfer +============================ +This module contains ... +""" + +from ._free_convection import FreeConvection + +__all__ = [ + "FreeConvection", +] diff --git a/src/felupe/constitution/heat_transfer/_free_convection.py b/src/felupe/constitution/heat_transfer/_free_convection.py new file mode 100644 index 000000000..0bee90e54 --- /dev/null +++ b/src/felupe/constitution/heat_transfer/_free_convection.py @@ -0,0 +1,117 @@ +# -*- coding: utf-8 -*- +""" +This file is part of FElupe. + +FElupe is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +FElupe is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with FElupe. If not, see . +""" + +import numpy as np + +from pyfluids import HumidAir, InputHumidAir +from scipy.constants import g + + +class FreeConvection: + r"Free convection heat transfer formulation for flat plates." + + def __init__(self): + self.T0 = 273.15 + self.P0 = 101325 + + + def _pyfluids_units(): + check = HumidAir().factory() + if str(check.units_system) == 'SIWithCelsiusAndPercents': + dt_ = 0 # use °C + rh_ = 1 # use % + else: + dt_ = 273.15 # use K + rh_ = 100 # use absolute value + return dt_, rh_ + + + @np.vectorize + def _rayleigh(ts_c, ti_c, length_, rh=10): + """ + Calculate dimensionless Rayleigh number Ra. + """ + dtk, rhf = self._pyfluids_units() + tm_c = (ts_c + ti_c)/2 + + # Humid air properties at p0 and Tm (indoors). + air = HumidAir().with_state( + InputHumidAir.pressure(self.P0), + InputHumidAir.temperature(tm_c + dtk), + InputHumidAir.relative_humidity(rh/rhf), + ) + rho = air.density + cp = air.specific_heat + uv = air.kinematic_viscosity # m^2/s + k = air.conductivity # W/(m K) + alpha = k/(rho*cp) # m^2/s thermal diffusivity + beta = 1/(tm_c + self.T0) + + # Eqn. 9.25, page 571. + ra = g*beta*abs(ts_c - ti_c)*length_*length_*length_/alpha/uv + + return ra + + @np.vectorize + def _nusselt_horizontal(ra, pr, hflux='z+'): + """ + Calculate dimensionless Nusselt number Nu for horizontal plates for various + cases of heat flux direction. + """ + if hflux == 'z+': # warm plate, top face or cold plate, bottom face + if ra < 1E04: + nu = 0.54*math.pow(1E04,0.25) + elif (1E04 <= ra <= 1E07) and (pr >= 0.7): + nu = 0.54*math.pow(ra,0.25) + elif 1E07 < ra <= 1E11: + nu = 0.15*math.pow(ra,0.33333) + else: + nu = 0.15*math.pow(1E11,0.33333) + else: # warm plate, bottom face or cold plate, top face + if ra < 1E04: + nu = 0.52*math.pow(1E04,0.2) + elif 1E04 <= ra <= 1E09 and pr >= 0.7: + nu = 0.52*math.pow(ra,0.2) + else: + nu = 0.52*math.pow(1E09,0.2) + return nu + + + @np.vectorize + def hc_fun(ts, tamb): + """ + Calculate convection coefficient for horizontal plate. + """ + dtk, rhf = self._pyfluids_units() + tm_c = (ts + tamb)/2 + rh = 50 + + air = HumidAir().with_state( + InputHumidAir.pressure(self.P0), + InputHumidAir.temperature(tm_c + dtk), + InputHumidAir.relative_humidity(rh/rhf), + ) + + l = 1.25 # assume slab size 5 x 5 m^2 for h_c + ra = rayleigh(ts, tamb, l) + if ts > tamb: + nu = nusselt_horizontal(ra, air.prandtl, hflux='z+') + else: + nu = nusselt_horizontal(ra, air.prandtl, hflux='z-') + + return(nu*air.conductivity/l) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 07f7a8c40..6b5cabd7d 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -2,105 +2,9 @@ """ Tester version for _solid_body_surface_convection.py example case. """ -import math -import numpy as np - -from pyfluids import HumidAir, InputHumidAir -from scipy.constants import g import felupe as fem - -T0 = 273.15 -P0 = 101325 - -def pyfluids_units(): - check = HumidAir().factory() - if str(check.units_system) == 'SIWithCelsiusAndPercents': - dt_ = 0 # use °C - rh_ = 1 # use % - else: - dt_ = 273.15 # use K - rh_ = 100 # use absolute value - return dt_, rh_ - - -@np.vectorize -def rayleigh(ts_c, ti_c, length_, rh=10): - """ - Calculate dimensionless Rayleigh number Ra. - """ - dtk, rhf = pyfluids_units() - tm_c = (ts_c + ti_c)/2 - - # Humid air properties at p0 and Tm (indoors). - air = HumidAir().with_state( - InputHumidAir.pressure(P0), - InputHumidAir.temperature(tm_c + dtk), - InputHumidAir.relative_humidity(rh/rhf), - ) - rho = air.density - cp = air.specific_heat - uv = air.kinematic_viscosity # m^2/s - k = air.conductivity # W/(m K) - alpha = k/(rho*cp) # m^2/s thermal diffusivity - beta = 1/(tm_c + T0) - - # Eqn. 9.25, page 571. - ra = g*beta*abs(ts_c - ti_c)*length_*length_*length_/alpha/uv - - return ra - -@np.vectorize -def nusselt_horizontal(ra, pr, hflux='z+'): - """ - Calculate dimensionless Nusselt number Nu for horizontal plates for various - cases of heat flux direction. - """ - if hflux == 'z+': # warm plate, top face or cold plate, bottom face - if ra < 1E04: - nu = 0.54*math.pow(1E04,0.25) - elif (1E04 <= ra <= 1E07) and (pr >= 0.7): - nu = 0.54*math.pow(ra,0.25) - elif 1E07 < ra <= 1E11: - nu = 0.15*math.pow(ra,0.33333) - else: - nu = 0.15*math.pow(1E11,0.33333) - else: # warm plate, bottom face or cold plate, top face - if ra < 1E04: - nu = 0.52*math.pow(1E04,0.2) - elif 1E04 <= ra <= 1E09 and pr >= 0.7: - nu = 0.52*math.pow(ra,0.2) - else: - nu = 0.52*math.pow(1E09,0.2) - return nu - - -@np.vectorize -def hc_fun(ts, tamb): - """ - Calculate convection coefficient for horizontal plate. - """ - dtk, rhf = pyfluids_units() - tm_c = (ts + tamb)/2 - rh = 50 - - air = HumidAir().with_state( - InputHumidAir.pressure(P0), - InputHumidAir.temperature(tm_c + dtk), - InputHumidAir.relative_humidity(rh/rhf), - ) - - l = 1.25 # assume slab size 5 x 5 m^2 for h_c - ra = rayleigh(ts, tamb, l) - if ts > tamb: - nu = nusselt_horizontal(ra, air.prandtl, hflux='z+') - else: - nu = nusselt_horizontal(ra, air.prandtl, hflux='z-') - - return(nu*air.conductivity/l) - - mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells region = fem.RegionQuad(mesh) temperature = fem.Field(region, dim=1, values=30.0) @@ -124,7 +28,7 @@ def hc_fun(ts, tamb): convection_function = fem.thermal.SolidBodySurfaceConvection( field=field_convection, - convection_coefficient=hc_fun, #(30, 20), + convection_coefficient=fem.FreeConvection.hc_fun, #(30, 20), temperature=20.0, # °C ) @@ -141,7 +45,7 @@ def callback(stepnumber, substepnumber, substep, tstep_data): tstep_data["ts_top.degC"].append(ts) tstep_data["hc_top.W.m-2.K-1"].append( convection_function.results.convection_coefficient.mean()) - tstep_data["hc_fun_top.W.m-2.K-1"].append(hc_fun(ts, tamb).mean()) + tstep_data["hc_fun_top.W.m-2.K-1"].append(fem.FreeConvection.hc_fun(ts, tamb).mean()) tstep_data["hc_calc_top.W.m-2.K-1"].append(abs(qc/(ts-tamb))) @@ -177,7 +81,6 @@ def callback(stepnumber, substepnumber, substep, tstep_data): # or # https://stackoverflow.com/questions/9103166/multiple-axis-in-matplotlib-with-different-scales import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import host_subplot # host = host_subplot(111) # par = host.twinx() From 1fac1c0a6018c93eb006d6c6b6c94eca5f73f5ee Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Thu, 4 Jun 2026 12:17:28 +0200 Subject: [PATCH 56/71] Address various typos and minor code issues. Code save. --- .../heat_transfer/_free_convection.py | 37 ++++++++++--------- src/felupe/thermal/_sbsc_tester.py | 6 ++- 2 files changed, 25 insertions(+), 18 deletions(-) diff --git a/src/felupe/constitution/heat_transfer/_free_convection.py b/src/felupe/constitution/heat_transfer/_free_convection.py index 0bee90e54..858d80a04 100644 --- a/src/felupe/constitution/heat_transfer/_free_convection.py +++ b/src/felupe/constitution/heat_transfer/_free_convection.py @@ -16,6 +16,7 @@ along with FElupe. If not, see . """ +import math import numpy as np from pyfluids import HumidAir, InputHumidAir @@ -25,12 +26,18 @@ class FreeConvection: r"Free convection heat transfer formulation for flat plates." - def __init__(self): + def __init__(self, plate_width, plate_length, rh=50): self.T0 = 273.15 self.P0 = 101325 + self.plate_width = plate_width + self.plate_length = plate_length + self.rh = rh + # Characteristic length for horizontal plate. + self.length = self.plate_width*self.plate_length/\ + (2*(self.plate_width+self.plate_length)) - def _pyfluids_units(): + def _pyfluids_units(self): check = HumidAir().factory() if str(check.units_system) == 'SIWithCelsiusAndPercents': dt_ = 0 # use °C @@ -41,8 +48,8 @@ def _pyfluids_units(): return dt_, rh_ - @np.vectorize - def _rayleigh(ts_c, ti_c, length_, rh=10): + # @np.vectorize + def _rayleigh(self, ts_c, ti_c, length_, rh=10): """ Calculate dimensionless Rayleigh number Ra. """ @@ -67,8 +74,7 @@ def _rayleigh(ts_c, ti_c, length_, rh=10): return ra - @np.vectorize - def _nusselt_horizontal(ra, pr, hflux='z+'): + def _nusselt_horizontal(self, ra, pr, hflux='z+'): """ Calculate dimensionless Nusselt number Nu for horizontal plates for various cases of heat flux direction. @@ -92,26 +98,23 @@ def _nusselt_horizontal(ra, pr, hflux='z+'): return nu - @np.vectorize - def hc_fun(ts, tamb): + def hc_fun(self, ts, tamb): """ Calculate convection coefficient for horizontal plate. """ dtk, rhf = self._pyfluids_units() tm_c = (ts + tamb)/2 - rh = 50 air = HumidAir().with_state( InputHumidAir.pressure(self.P0), InputHumidAir.temperature(tm_c + dtk), - InputHumidAir.relative_humidity(rh/rhf), + InputHumidAir.relative_humidity(self.rh/rhf), ) - - l = 1.25 # assume slab size 5 x 5 m^2 for h_c - ra = rayleigh(ts, tamb, l) + + ra = self._rayleigh(ts_c=ts, ti_c=tamb, length_=self.length) if ts > tamb: - nu = nusselt_horizontal(ra, air.prandtl, hflux='z+') + nu = self._nusselt_horizontal(ra, air.prandtl, hflux='z+') else: - nu = nusselt_horizontal(ra, air.prandtl, hflux='z-') - - return(nu*air.conductivity/l) + nu = self._nusselt_horizontal(ra, air.prandtl, hflux='z-') + + return(nu*air.conductivity/self.length) diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 6b5cabd7d..384c1cd51 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -3,8 +3,12 @@ Tester version for _solid_body_surface_convection.py example case. """ +import numpy as np import felupe as fem +hc_fun = np.vectorize(fem.FreeConvection(5, 5).hc_fun) + + mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells region = fem.RegionQuad(mesh) temperature = fem.Field(region, dim=1, values=30.0) @@ -28,7 +32,7 @@ convection_function = fem.thermal.SolidBodySurfaceConvection( field=field_convection, - convection_coefficient=fem.FreeConvection.hc_fun, #(30, 20), + convection_coefficient=hc_fun, temperature=20.0, # °C ) From 2cb06250f8a2599df6e9fda1520945ca3dea4e80 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Thu, 4 Jun 2026 12:28:52 +0200 Subject: [PATCH 57/71] Free convection function via class FreeConvection now works in tester code as expected. --- src/felupe/constitution/heat_transfer/_free_convection.py | 2 -- src/felupe/thermal/_sbsc_tester.py | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/felupe/constitution/heat_transfer/_free_convection.py b/src/felupe/constitution/heat_transfer/_free_convection.py index 858d80a04..0a762c687 100644 --- a/src/felupe/constitution/heat_transfer/_free_convection.py +++ b/src/felupe/constitution/heat_transfer/_free_convection.py @@ -17,7 +17,6 @@ """ import math -import numpy as np from pyfluids import HumidAir, InputHumidAir from scipy.constants import g @@ -48,7 +47,6 @@ def _pyfluids_units(self): return dt_, rh_ - # @np.vectorize def _rayleigh(self, ts_c, ti_c, length_, rh=10): """ Calculate dimensionless Rayleigh number Ra. diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index 384c1cd51..fe3ccd4ae 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -49,7 +49,7 @@ def callback(stepnumber, substepnumber, substep, tstep_data): tstep_data["ts_top.degC"].append(ts) tstep_data["hc_top.W.m-2.K-1"].append( convection_function.results.convection_coefficient.mean()) - tstep_data["hc_fun_top.W.m-2.K-1"].append(fem.FreeConvection.hc_fun(ts, tamb).mean()) + tstep_data["hc_fun_top.W.m-2.K-1"].append(hc_fun(ts, tamb).mean()) tstep_data["hc_calc_top.W.m-2.K-1"].append(abs(qc/(ts-tamb))) From 52f55fef1703c38f12b8cf769dd157d3be981894 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Thu, 4 Jun 2026 15:55:45 +0200 Subject: [PATCH 58/71] Adapt surface radiation class definition to calculate radiation coefficient as result. Seems to work. --- .../thermal/_solidbody_surface_radiation.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/felupe/thermal/_solidbody_surface_radiation.py b/src/felupe/thermal/_solidbody_surface_radiation.py index ad6546897..97338b69d 100644 --- a/src/felupe/thermal/_solidbody_surface_radiation.py +++ b/src/felupe/thermal/_solidbody_surface_radiation.py @@ -123,6 +123,8 @@ def __init__(self, field, emissivity, temperature): self.results = Results() self.results.temperature = temperature # ambient temperature in °C self.results.emissivity = emissivity + self.results.radiation_coefficient =\ + 4 * emissivity * sigma * (temperature + 273.15) ** 3 self._sigma = sigma # Stefan-Boltzmann constant @@ -171,11 +173,16 @@ def _matrix(self, field=None, **kwargs): dim = self.field[0].dim temperature = self.field.extract(grad=False)[0] + self.results.radiation_coefficient = 4 * self.results.emissivity\ + * self._sigma\ + * ((temperature + self.results.temperature)/2 + 273.15) ** 3 + fun = [ - -self.results.emissivity - * self._sigma - * 4 - * (temperature + 273.15) ** 3 + -self.results.radiation_coefficient + # -self.results.emissivity + # * self._sigma + # * 4 + # * (temperature + 273.15) ** 3 * np.eye(dim).reshape(dim, dim, 1, 1) ] From 66680cb32d2953ef619b0ab5ccc6cd1d27e99741 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Thu, 4 Jun 2026 15:56:19 +0200 Subject: [PATCH 59/71] Extend example 23 to use detailed convection and radiation transfer models. Seems to work. --- examples/ex23_solid_body_thermal-tabs.py | 84 +++++++++++++++++------- 1 file changed, 60 insertions(+), 24 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index 88f9583fb..b657aaf81 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -18,11 +18,11 @@ This example describes a thermally activated concrete slab using a simplified model and geometry. The model is two-dimensional. The system is set up with two :class:`solids `. The temperature boundary -conditions include the floor temperature (constant), the ceiling temperature -(constant) and the room air temperature with a :math:`\pm 2` K sinusoidal -variation around its average value with a period of 24 h. +conditions include the floor temperature, the ceiling temperature and the room +air temperatures, each with a :math:`\pm \Delta\theta` K sinusoidal variation +around its average value with a period of 24 h. -The heat injection via the pipe layer is constant at 295 W/m2 and directly +The heat injection via the pipe layer is constant at 231 W/m2 and directly injected at the internal concrete surfaces (no pipe material is modelled). Surface heat transfer is modelled separately for convection and radiation. @@ -140,14 +140,29 @@ coefficient=7.69, # W/(m^2 K) temperature=20.0, # °C ) -top_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + +hc_fun = np.vectorize(fem.FreeConvection(5, 5).hc_fun) + +top_convection = fem.thermal.SolidBodySurfaceConvection( + field=top_field, + convection_coefficient=hc_fun, # W/(m^2 K) + temperature=20.0, # °C +) +bottom_convection = fem.thermal.SolidBodySurfaceConvection( + field=bottom_field, + convection_coefficient=hc_fun, # W/(m^2 K) + temperature=20.0, # °C +) + +top_radiation = fem.thermal.SolidBodySurfaceRadiation( field=top_field, - coefficient=10.0, # W/(m^2 K) + emissivity=0.9, temperature=20.0, # °C ) -bottom_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( + +bottom_radiation = fem.thermal.SolidBodySurfaceRadiation( field=bottom_field, - coefficient=5.0, # W/(m^2 K) + emissivity=0.9, temperature=20.0, # °C ) @@ -173,22 +188,36 @@ # by the :meth:`~felupe.thermal.SolidBodyThermal.heat_flux_boundary` method of the # thermal solid body, which returns the integrated surface heat flux for a given # boundary region and time step. The mean surface heat flux is stored in the -# ``flux_data`` dictionary, which is passed to the callback function as an argument. -def callback(stepnumber, substepnumber, substep, flux_data): +# ``tstep_data`` dictionary, which is passed to the callback function as an argument. +def callback(stepnumber, substepnumber, substep, tstep_data): """Save mean surface heat flux at internal and external boundaries.""" heat_flux = materials[0].heat_flux_boundary - flux_data["top"].append(heat_flux(region=top_region)) - flux_data["bottom"].append(heat_flux(region=bottom_region)) + tstep_data["top"].append(heat_flux(region=top_region)) + tstep_data["bottom"].append(heat_flux(region=bottom_region)) + + tstep_data["hc_top.W.m-2.K-1"].append( + top_convection.results.convection_coefficient.mean()) + tstep_data["hr_top.W.m-2.K-1"].append( + top_radiation.results.radiation_coefficient.mean()) + + tstep_data["hc_bottom.W.m-2.K-1"].append( + bottom_convection.results.convection_coefficient.mean()) + tstep_data["hr_bottom.W.m-2.K-1"].append( + bottom_radiation.results.radiation_coefficient.mean()) + pflux = 0 for p_ in pipe_region: pflux += heat_flux(region=p_) - flux_data["pipes"].append(pflux) + tstep_data["pipes"].append(pflux) -time_steps = fem.math.linsteps([0, 24 * 3600], num=int(24 * 3600 / 720))[1:] +N_DAYS = 2 +time_steps = fem.math.linsteps([0, N_DAYS * 24 * 3600], + num=int(N_DAYS * 24 * 3600 / 720))[1:] t_air = 20 + 2 * np.sin(2 * np.pi * time_steps / 86400) -# t_sur = 20 + 0.5 * np.sin(2 * np.pi * time_steps / 86400) +t_ceil = 20 + 0.5 * np.sin(2 * np.pi * time_steps / 86400) +t_floor = 18 + 0.5 * np.sin(2 * np.pi * time_steps / 86400) # %% @@ -201,24 +230,29 @@ def callback(stepnumber, substepnumber, substep, flux_data): # data in the result file. # model_list = [*materials, top_heat_transfer, bottom_heat_transfer] model_list = [*materials, side1_heat_transfer, side2_heat_transfer, - top_heat_transfer, bottom_heat_transfer] + top_convection, bottom_convection, top_radiation, bottom_radiation] time = fem.thermal.TimeStep(model_list) ramp = { time: time_steps, side1_heat_transfer: t_air, side2_heat_transfer: t_air, - top_heat_transfer: t_air, - bottom_heat_transfer: t_air, + top_convection: t_air, + bottom_convection: t_air, + top_radiation: t_ceil, + bottom_radiation: t_floor, } step = fem.Step( items=[time] + model_list + pipe_flux, ramp=ramp, ) -flux_data = {"top": [], "bottom": [], "pipes": []} +tstep_data = {"top": [], "bottom": [], + "hc_top.W.m-2.K-1": [], "hr_top.W.m-2.K-1": [], + "hc_bottom.W.m-2.K-1": [], "hr_bottom.W.m-2.K-1": [], + "pipes": []} -job = fem.Job(steps=[step], callback=callback, flux_data=flux_data).evaluate( +job = fem.Job(steps=[step], callback=callback, tstep_data=tstep_data).evaluate( x0=field, filename="result.xdmf", # create a result file for Paraview point_data={"Temperature": lambda field, substep: temperature.values}, @@ -227,7 +261,7 @@ def callback(stepnumber, substepnumber, substep, flux_data): ) # %% -# Internal and external surface heat flux values are plotted over time. +# Top and bottom surface heat flux values are plotted over time. # # .. note:: # @@ -236,8 +270,8 @@ def callback(stepnumber, substepnumber, substep, flux_data): # the internal surface). # fig, ax = plt.subplots() -ax.plot(time_steps / 3600, flux_data["top"], color="C3", label="top") -ax.plot(time_steps / 3600, flux_data["bottom"], color="C0", label="bottom") +ax.plot(time_steps / 3600, tstep_data["top"], color="C3", label="top") +ax.plot(time_steps / 3600, tstep_data["bottom"], color="C0", label="bottom") tmin, tmax = ax.get_xlim() ax.plot([tmin, tmax], np.zeros(2), "black", lw=0.5) @@ -249,8 +283,10 @@ def callback(stepnumber, substepnumber, substep, flux_data): ax.legend() ax.set(xlim=(tmin, tmax), xlabel="time in h", ylabel=r"surface heat flux in W/m$^2$") +plt.savefig("_ex23.png") + # %% # A view on the temperature field at the end of the simulation period visualizes the # temperature distribution. -field.plot("Field", scalar_bar_vertical=True).show() +# field.plot("Field", scalar_bar_vertical=True).show() From 53085c26a76ba122fc3f4dd43017c4929bea5ab6 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Thu, 4 Jun 2026 15:56:43 +0200 Subject: [PATCH 60/71] Start writing documentation for FreeConvection class. --- .../heat_transfer/_free_convection.py | 30 ++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/src/felupe/constitution/heat_transfer/_free_convection.py b/src/felupe/constitution/heat_transfer/_free_convection.py index 0a762c687..177b83ae9 100644 --- a/src/felupe/constitution/heat_transfer/_free_convection.py +++ b/src/felupe/constitution/heat_transfer/_free_convection.py @@ -23,7 +23,35 @@ class FreeConvection: - r"Free convection heat transfer formulation for flat plates." + r"""Free convection heat transfer formulation for flat plates. + + Parameters + ---------- + plate_width : float + Horizontal plate width in (m). + plate_length : float + Horizontal plate length in (m). + rh : float (optional, default 50 %) + Relative humidity of air in (%). + + Notes + ----- + This class represents ... + + References + ---------- + .. [1] F. P. Incropera, D. P. DeWitt, and et. al., Fundamentals of Heat + and Mass Transfer, 6th Edition. John Wiley & Sons, 2007; + ISBN 0-471-45728-0. + + See Also + -------- + felupe.thermal.SolidBodyThermal : A thermal solid body for heat conduction. + felupe.thermal.SolidBodySurfaceConvection : Detailed surface convection + heat transfer. + + + """ def __init__(self, plate_width, plate_length, rh=50): self.T0 = 273.15 From 572b030dec36ca8e6367d694fa1a059ee55836c8 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Thu, 4 Jun 2026 17:09:07 +0200 Subject: [PATCH 61/71] Add 2nd day with heat extraction via hydronic system. Seems to work. --- examples/ex23_solid_body_thermal-tabs.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index b657aaf81..f981f0840 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -174,6 +174,7 @@ pipe_field = [] pipe_flux = [] for idx, p in enumerate(center_points): + # Inelegant, but seems to work: mask = np.isclose(mesh.points[:, None, :], p[:], rtol=0.05, atol=0.0101).all(axis=2).any(axis=1) pipe_region.append(fem.RegionQuadBoundary(mesh, mask=mask)) pipe_field.append(fem.FieldContainer([fem.Field(pipe_region[idx], dim=1)])) @@ -219,6 +220,11 @@ def callback(stepnumber, substepnumber, substep, tstep_data): t_ceil = 20 + 0.5 * np.sin(2 * np.pi * time_steps / 86400) t_floor = 18 + 0.5 * np.sin(2 * np.pi * time_steps / 86400) +pipe_heat_flux = np.concatenate( + (fem.math.linsteps([-231.25, -231.25], num=int(len(time_steps)/2)-1), + fem.math.linsteps([231.25, 231.25], num=int(len(time_steps)/2)-1)) +) + # %% # The time step item is created with the thermal solid bodies. It must be located as the @@ -241,6 +247,10 @@ def callback(stepnumber, substepnumber, substep, tstep_data): bottom_convection: t_air, top_radiation: t_ceil, bottom_radiation: t_floor, + pipe_flux[0]: pipe_heat_flux, + pipe_flux[1]: pipe_heat_flux, + pipe_flux[2]: pipe_heat_flux, + pipe_flux[3]: pipe_heat_flux, } step = fem.Step( items=[time] + model_list + pipe_flux, From 1641db6f219e3421bbe843a2686af2fc1891fe88 Mon Sep 17 00:00:00 2001 From: Achim Geissler Date: Fri, 5 Jun 2026 16:39:46 +0200 Subject: [PATCH 62/71] Code save. --- examples/ex23_solid_body_thermal-tabs.py | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index f981f0840..324c3d48b 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -279,6 +279,7 @@ def callback(stepnumber, substepnumber, substep, tstep_data): # external surface), and **negative** when **heat enters the construction** (here, on # the internal surface). # +# "engine mount, ex07" => further graphic examples. fig, ax = plt.subplots() ax.plot(time_steps / 3600, tstep_data["top"], color="C3", label="top") ax.plot(time_steps / 3600, tstep_data["bottom"], color="C0", label="bottom") From 28bff45741ccd77bd2eb4901d70e2393690b95a3 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Tue, 9 Jun 2026 22:07:58 +0200 Subject: [PATCH 63/71] Add heat transfer and pipe flux figure to ex23. Adapt code documentation (WIP). --- examples/ex23_solid_body_thermal-tabs.py | 75 +++++++++++++++++++----- src/felupe/thermal/_sbsc_tester.py | 27 ++------- 2 files changed, 66 insertions(+), 36 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index 324c3d48b..98039e0a6 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -112,8 +112,7 @@ ) # %% -# The surface heat transfer coefficients and ambient temperatures are defined -# for the side, top and bottom surfaces. +# The surface heat transfer is defined for the side, top and bottom surfaces. side_region1 = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.min()) side_temperature1 = fem.Field(side_region1, dim=1) side_field1 = fem.FieldContainer([side_temperature1]) @@ -130,6 +129,7 @@ top_temperature = fem.Field(top_region, dim=1) top_field = fem.FieldContainer([top_temperature]) +# For the sides, combined transfer coefficients are used. side1_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer( field=side_field1, coefficient=7.69, # W/(m^2 K) @@ -141,6 +141,12 @@ temperature=20.0, # °C ) +# For the top and bottom surfaces, the detailed calculation approaches defined +# in :class:`~felupe.thermal.SolidBodySurfaceConvection` and +# :class:`~felupe.thermal.SolidBodySurfaceRadiation` are used for convection +# and radiation, respectively. For convection, the convection coefficient +# function defined in :class:`~felupe.constitution.heat_transfer.FreeConvection` +# is used. hc_fun = np.vectorize(fem.FreeConvection(5, 5).hc_fun) top_convection = fem.thermal.SolidBodySurfaceConvection( @@ -184,11 +190,19 @@ )) # %% -# A callback-function records the mean surface heat flux at the internal and external -# boundaries after each completed time step. The mean surface heat flux is calculated -# by the :meth:`~felupe.thermal.SolidBodyThermal.heat_flux_boundary` method of the +# A callback-function records the mean surface heat flux at the top and bottom +# boundaries, the top and bottom convection coefficients as well as the top and +# bottom radiation coefficients after each completed time step. +# The mean surface heat flux is calculated by the +# :meth:`~felupe.thermal.SolidBodyThermal.heat_flux_boundary` method of the # thermal solid body, which returns the integrated surface heat flux for a given -# boundary region and time step. The mean surface heat flux is stored in the +# boundary region and time step. +# +# The convection coefficient values are calculated by ... +# +# The radiation coefficient values are calculated by ... +# +# All values are stored in the # ``tstep_data`` dictionary, which is passed to the callback function as an argument. def callback(stepnumber, substepnumber, substep, tstep_data): """Save mean surface heat flux at internal and external boundaries.""" @@ -234,7 +248,6 @@ def callback(stepnumber, substepnumber, substep, tstep_data): # callback function, and evaluated with the top-level temperature field. A result file # is created for visualization in Paraview, and the temperature field is saved as point- # data in the result file. -# model_list = [*materials, top_heat_transfer, bottom_heat_transfer] model_list = [*materials, side1_heat_transfer, side2_heat_transfer, top_convection, bottom_convection, top_radiation, bottom_radiation] @@ -275,11 +288,10 @@ def callback(stepnumber, substepnumber, substep, tstep_data): # # .. note:: # -# The heat flux is **positive** when **heat leaves the construction** (here, on the -# external surface), and **negative** when **heat enters the construction** (here, on -# the internal surface). -# -# "engine mount, ex07" => further graphic examples. +# The heat flux is **positive** when **heat leaves the construction** (here, +# on both top and bottom surfaces in 'heating mode', and **negative** when +# **heat enters the construction** (here, on both the top and bottom +# surfaces in 'cooling mode'. fig, ax = plt.subplots() ax.plot(time_steps / 3600, tstep_data["top"], color="C3", label="top") ax.plot(time_steps / 3600, tstep_data["bottom"], color="C0", label="bottom") @@ -294,10 +306,43 @@ def callback(stepnumber, substepnumber, substep, tstep_data): ax.legend() ax.set(xlim=(tmin, tmax), xlabel="time in h", ylabel=r"surface heat flux in W/m$^2$") -plt.savefig("_ex23.png") +plt.savefig("_ex23a.png") + +# %% +# Top and bottom surface heat transfer coefficients and pipe heat flux are +# plotted over time. +fig, ax = plt.subplots() +fig.subplots_adjust(right=0.75) + +twin1 = ax.twinx() +twin2 = ax.twinx() + +ax.set_xlabel("Time (s)") +ax.set_ylabel("Convection coefficient (W/(m2 K))") +twin1.set_ylabel("Temperature (°C)") +twin2.set_ylabel("Pipe heat flux (W/m2)") + +time_steps_h = time_steps / 3600 + +p1 = ax.plot(time_steps_h, tstep_data["hc_top.W.m-2.K-1"], + label="hc_top", color='lightblue') +p2 = ax.plot(time_steps_h, tstep_data["hc_bottom.W.m-2.K-1"], + label="hc_bottom", color='darkblue') +p3 = twin1.plot(time_steps_h, tstep_data["hr_top.W.m-2.K-1"], + label="hr_top", color='blue') +p4 = twin1.plot(time_steps_h, tstep_data["hr_bottom.W.m-2.K-1"], + label="hr_bottom", color='red') +p5 = twin2.plot(time_steps_h, tstep_data["pipes"], + label="pipe_flux", color='magenta') + +ax.legend(handles=p1+p2+p3+p4+p5, labelcolor="linecolor") + +twin2.spines['right'].set_position(('outward', 45)) + +plt.savefig("_ex23b.png") # %% -# A view on the temperature field at the end of the simulation period visualizes the -# temperature distribution. +# A view on the temperature field at the end of the simulation period visualizes +# the temperature distribution. # field.plot("Field", scalar_bar_vertical=True).show() diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py index fe3ccd4ae..588d48f45 100644 --- a/src/felupe/thermal/_sbsc_tester.py +++ b/src/felupe/thermal/_sbsc_tester.py @@ -86,45 +86,30 @@ def callback(stepnumber, substepnumber, substep, tstep_data): # https://stackoverflow.com/questions/9103166/multiple-axis-in-matplotlib-with-different-scales import matplotlib.pyplot as plt -# host = host_subplot(111) -# par = host.twinx() -# par2 = host.twinx() - - fig, ax = plt.subplots() fig.subplots_adjust(right=0.75) twin1 = ax.twinx() twin2 = ax.twinx() -# twin2.spines.right.set_position(("axes", 1.2)) - ax.set_xlabel("Time (s)") ax.set_ylabel("Convection coefficient (W/(m2 K))") twin1.set_ylabel("Temperature (°C)") twin2.set_ylabel("Heat flux (W/m2)") p1 = ax.plot(tstep_data["tstep.s"], tstep_data["hc_top.W.m-2.K-1"], - label="hc") + label="hc", color='lightblue') p2 = ax.plot(tstep_data["tstep.s"], tstep_data["hc_fun_top.W.m-2.K-1"], - label="hc_fun(ts_top - t_amb)") + label="hc_fun(ts_top - t_amb)", color='darkblue') p3 = twin1.plot(tstep_data["tstep.s"], tstep_data["tamb.degC"], - label="t_amb") + label="t_amb", color='blue') p4 = twin1.plot(tstep_data["tstep.s"], tstep_data["ts_top.degC"], - label="ts_top") -# p3, = host.plot(tstep_data["tstep.s"], tstep_data["hc_calc_top.W.m-2.K-1"], -# label="hc from q") -# p3.set_color(p2.get_color()) + label="ts_top", color='red') p5 = twin2.plot(tstep_data["tstep.s"], tstep_data["qc_top.W.m-2"], - label="qc_top") + label="qc_top", color='magenta') -ax.legend(handles=p1+p2+p3+p4+p5, loc='best') - -ax.legend(labelcolor="linecolor") +ax.legend(handles=p1+p2+p3+p4+p5, labelcolor="linecolor") twin2.spines['right'].set_position(('outward', 45)) -# host.yaxis.label.set_color(p1.get_color()) -# par.yaxis.label.set_color(p2.get_color()) - plt.savefig("_test.png") From f352f0155f65adb8258bd35a8836a96c9050530d Mon Sep 17 00:00:00 2001 From: AGeissler Date: Thu, 11 Jun 2026 11:52:20 +0200 Subject: [PATCH 64/71] Add absolute pressure as variable. Identify issue with too similar values for top/bottom of plate (missing switch!). --- .../heat_transfer/_free_convection.py | 28 +++++++++++-------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/felupe/constitution/heat_transfer/_free_convection.py b/src/felupe/constitution/heat_transfer/_free_convection.py index 177b83ae9..c226f8873 100644 --- a/src/felupe/constitution/heat_transfer/_free_convection.py +++ b/src/felupe/constitution/heat_transfer/_free_convection.py @@ -25,12 +25,16 @@ class FreeConvection: r"""Free convection heat transfer formulation for flat plates. + <> + Parameters ---------- plate_width : float Horizontal plate width in (m). plate_length : float Horizontal plate length in (m). + p_abs: float (optional, default 101325) + Absolute (total) air pressure in (Pa). rh : float (optional, default 50 %) Relative humidity of air in (%). @@ -53,12 +57,12 @@ class FreeConvection: """ - def __init__(self, plate_width, plate_length, rh=50): - self.T0 = 273.15 - self.P0 = 101325 + def __init__(self, plate_width, plate_length, p_abs=101325, rh=50): + self.T0 = 273.15 # 0 °C in Kelvin, for °C <=> K conversion self.plate_width = plate_width self.plate_length = plate_length self.rh = rh + self.p_abs = p_abs # Characteristic length for horizontal plate. self.length = self.plate_width*self.plate_length/\ @@ -84,7 +88,7 @@ def _rayleigh(self, ts_c, ti_c, length_, rh=10): # Humid air properties at p0 and Tm (indoors). air = HumidAir().with_state( - InputHumidAir.pressure(self.P0), + InputHumidAir.pressure(self.p_abs), InputHumidAir.temperature(tm_c + dtk), InputHumidAir.relative_humidity(rh/rhf), ) @@ -107,20 +111,20 @@ def _nusselt_horizontal(self, ra, pr, hflux='z+'): """ if hflux == 'z+': # warm plate, top face or cold plate, bottom face if ra < 1E04: - nu = 0.54*math.pow(1E04,0.25) + nu = 0.54*math.pow(1E04, 0.25) elif (1E04 <= ra <= 1E07) and (pr >= 0.7): - nu = 0.54*math.pow(ra,0.25) + nu = 0.54*math.pow(ra, 0.25) elif 1E07 < ra <= 1E11: - nu = 0.15*math.pow(ra,0.33333) + nu = 0.15*math.pow(ra, 0.33333) else: - nu = 0.15*math.pow(1E11,0.33333) + nu = 0.15*math.pow(1E11, 0.33333) else: # warm plate, bottom face or cold plate, top face if ra < 1E04: - nu = 0.52*math.pow(1E04,0.2) + nu = 0.52*math.pow(1E04, 0.2) elif 1E04 <= ra <= 1E09 and pr >= 0.7: - nu = 0.52*math.pow(ra,0.2) + nu = 0.52*math.pow(ra, 0.2) else: - nu = 0.52*math.pow(1E09,0.2) + nu = 0.52*math.pow(1E09, 0.2) return nu @@ -132,7 +136,7 @@ def hc_fun(self, ts, tamb): tm_c = (ts + tamb)/2 air = HumidAir().with_state( - InputHumidAir.pressure(self.P0), + InputHumidAir.pressure(self.p_abs), InputHumidAir.temperature(tm_c + dtk), InputHumidAir.relative_humidity(self.rh/rhf), ) From fb1cc38ea6ab243d1f30fd83466fd387911652dc Mon Sep 17 00:00:00 2001 From: AGeissler Date: Thu, 11 Jun 2026 16:35:14 +0200 Subject: [PATCH 65/71] Add 'side of plate' toggle. Extend documentation of class (not finished). --- .../heat_transfer/_free_convection.py | 53 ++++++++++++++++--- 1 file changed, 45 insertions(+), 8 deletions(-) diff --git a/src/felupe/constitution/heat_transfer/_free_convection.py b/src/felupe/constitution/heat_transfer/_free_convection.py index c226f8873..4607441ff 100644 --- a/src/felupe/constitution/heat_transfer/_free_convection.py +++ b/src/felupe/constitution/heat_transfer/_free_convection.py @@ -25,14 +25,14 @@ class FreeConvection: r"""Free convection heat transfer formulation for flat plates. - <> - Parameters ---------- plate_width : float Horizontal plate width in (m). plate_length : float Horizontal plate length in (m). + plate_side : string + Face of plate considered, 'top' or 'bottom'. p_abs: float (optional, default 101325) Absolute (total) air pressure in (Pa). rh : float (optional, default 50 %) @@ -40,7 +40,36 @@ class FreeConvection: Notes ----- - This class represents ... + This class provides a convection heat transfer coefficient for horizontal + plates based on detailed empirical approaches from [1]_. + + The dimensionless Rayleigh number Ra is a function of gravity g, inverse + mean temperature (film temperature) :math:`\beta` (see :eq:`film-temperature`), + characteristic length L, thermal diffusivity :math:`\alpha` and kinematic + viscosity of the air :math:`\nu_v` according to :eq:`rayleigh`. + + .. math:: + :label: rayleigh + + Ra\,=\,\frac{g\,\beta\,|\theta_s - \theta_i|\,L^3}{\alpha\nu_v} + + where + + .. math:: + :label: film-temperature + + \beta\,=\,\frac{2}{T_s + T_i}\,\text{ in K}^{-1}. + + The characteristic length :math:`L` for horizontal plates is defined as + (eqn 9.29 from [1]_) + + .. math:: + :label: l-horiz-plate + + L\,=\,\frac{A_s}{P} + + where :math:`A_s` is the plate surface (one side) and :math:`P` is the + plate perimeter. References ---------- @@ -57,10 +86,12 @@ class FreeConvection: """ - def __init__(self, plate_width, plate_length, p_abs=101325, rh=50): + def __init__(self, plate_width, plate_length, plate_side, + p_abs=101325, rh=50): self.T0 = 273.15 # 0 °C in Kelvin, for °C <=> K conversion self.plate_width = plate_width self.plate_length = plate_length + self.plate_side = plate_side self.rh = rh self.p_abs = p_abs @@ -142,9 +173,15 @@ def hc_fun(self, ts, tamb): ) ra = self._rayleigh(ts_c=ts, ti_c=tamb, length_=self.length) - if ts > tamb: - nu = self._nusselt_horizontal(ra, air.prandtl, hflux='z+') - else: - nu = self._nusselt_horizontal(ra, air.prandtl, hflux='z-') + if self.plate_side == 'top': + if ts > tamb: + nu = self._nusselt_horizontal(ra, air.prandtl, hflux='z+') + else: + nu = self._nusselt_horizontal(ra, air.prandtl, hflux='z-') + else: # self.plate_side == 'bottom' (or actually any other string) + if ts < tamb: + nu = self._nusselt_horizontal(ra, air.prandtl, hflux='z+') + else: + nu = self._nusselt_horizontal(ra, air.prandtl, hflux='z-') return(nu*air.conductivity/self.length) From a3690d1cac164fcf4aa58d386a22308cfb7f9eb9 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Thu, 11 Jun 2026 16:36:14 +0200 Subject: [PATCH 66/71] Update ex23 to use correct plate sides for convection. Results for h_c are now different for top and bottom sides, as expected. --- examples/ex23_solid_body_thermal-tabs.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index 98039e0a6..be847ffa7 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -147,16 +147,17 @@ # and radiation, respectively. For convection, the convection coefficient # function defined in :class:`~felupe.constitution.heat_transfer.FreeConvection` # is used. -hc_fun = np.vectorize(fem.FreeConvection(5, 5).hc_fun) +hc_top = np.vectorize(fem.FreeConvection(5, 5, 'top').hc_fun) +hc_bottom = np.vectorize(fem.FreeConvection(5, 5, 'bottom').hc_fun) top_convection = fem.thermal.SolidBodySurfaceConvection( field=top_field, - convection_coefficient=hc_fun, # W/(m^2 K) + convection_coefficient=hc_top, # W/(m^2 K) temperature=20.0, # °C ) bottom_convection = fem.thermal.SolidBodySurfaceConvection( field=bottom_field, - convection_coefficient=hc_fun, # W/(m^2 K) + convection_coefficient=hc_bottom, # W/(m^2 K) temperature=20.0, # °C ) From 5c4161be56d93e723728e8579c2ede65ff6dc997 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Thu, 11 Jun 2026 16:40:11 +0200 Subject: [PATCH 67/71] Add reference to class 'free_convection'. --- src/felupe/thermal/_solidbody_surface_convection.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/felupe/thermal/_solidbody_surface_convection.py b/src/felupe/thermal/_solidbody_surface_convection.py index 14c97fbcc..ad2a630c3 100644 --- a/src/felupe/thermal/_solidbody_surface_convection.py +++ b/src/felupe/thermal/_solidbody_surface_convection.py @@ -106,6 +106,11 @@ class SolidBodySurfaceConvection: ... return(nu*lam_air/l) >>> hc_fun = np.vectorize(_hc_fun) + using constant air properties for brevity. + In :class:`~felupe.constitution.heat_transfer.free_convection` a more + detailed calculation of convection coefficients for horizontal plates is + given which use air properties based on current temperatures. + Set up the model (a horizontal slab with dimensions 1 x 1 m^2, 0.25 m thick). .. pyvista-plot:: From f79726a3e1884991947df300a3665b33bdbbed9f Mon Sep 17 00:00:00 2001 From: AGeissler Date: Thu, 11 Jun 2026 16:48:34 +0200 Subject: [PATCH 68/71] Further extend documentation of class free_convection. WIP. --- .../constitution/heat_transfer/_free_convection.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/felupe/constitution/heat_transfer/_free_convection.py b/src/felupe/constitution/heat_transfer/_free_convection.py index 4607441ff..35907b204 100644 --- a/src/felupe/constitution/heat_transfer/_free_convection.py +++ b/src/felupe/constitution/heat_transfer/_free_convection.py @@ -42,6 +42,13 @@ class FreeConvection: ----- This class provides a convection heat transfer coefficient for horizontal plates based on detailed empirical approaches from [1]_. + + The convection coefficient is calculated according to :eq:`hc`. + + .. math:: + :label: hc + + h_c\,=\,\text{Nu}\,\lambda_\text{air}\frac{1}{L} The dimensionless Rayleigh number Ra is a function of gravity g, inverse mean temperature (film temperature) :math:`\beta` (see :eq:`film-temperature`), @@ -61,7 +68,7 @@ class FreeConvection: \beta\,=\,\frac{2}{T_s + T_i}\,\text{ in K}^{-1}. The characteristic length :math:`L` for horizontal plates is defined as - (eqn 9.29 from [1]_) + given in :eq:`l-horiz-plate` (eqn 9.29 from [1]_). .. math:: :label: l-horiz-plate From 9f2c682590526afc5e342ac651ca891ef1a5e546 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sat, 13 Jun 2026 13:09:55 +0200 Subject: [PATCH 69/71] Clean up ex23. --- examples/ex23_solid_body_thermal-tabs.py | 77 +++++++++++------------- 1 file changed, 36 insertions(+), 41 deletions(-) diff --git a/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py index be847ffa7..187dbc1f1 100644 --- a/examples/ex23_solid_body_thermal-tabs.py +++ b/examples/ex23_solid_body_thermal-tabs.py @@ -9,10 +9,14 @@ :class:`~felupe.thermal.SolidBodySurfaceRadiation` and :class:`~felupe.thermal.SolidBodySurfaceConvection` - * evaluate the surface heat flux at internal and external boundaries + * calculate detailed convection transfer coefficient using + :class:`~felupe.constitution.heat_transfer.FreeConvection` + + * evaluate the surface heat flux at top and bottom boundaries with a job :class:`~felupe.Plugin` - * view the temperature field + * view top/bottom surface heat flux, convective and radiative transfer + coefficients and the temperature field This example describes a thermally activated concrete slab using a simplified @@ -25,7 +29,8 @@ The heat injection via the pipe layer is constant at 231 W/m2 and directly injected at the internal concrete surfaces (no pipe material is modelled). -Surface heat transfer is modelled separately for convection and radiation. +Surface heat transfer at the top and bottom surfaces is modelled separately for +convection and radiation. """ import matplotlib.pyplot as plt import numpy as np @@ -33,14 +38,14 @@ import felupe as fem # %% -# Define material properties as lists for (reinforced) concrete and insulation. +# Material properties are defined as lists for (reinforced) concrete and insulation. # This includes mass density, specific heat capacity and thermal conductivity. density = [2100, 20] # kg/m^3 specific_heat = [1000, 1450] # J/(kg K) thermal_conductivity = [2.1, 0.035] # W/(m K) # %% -# Set up one mesh per material. If a material consists of multiple areas, these +# One mesh per material is set up. If a material consists of multiple areas, these # are collected in a :class:`mesh container ` and are # merged into one mesh per material. These meshes per material are then added # to a mesh container for the construction. @@ -80,24 +85,21 @@ container = fem.MeshContainer([concrete, insulation], merge=True, decimals=6) -# container.plot( -# colors=["lightgrey", "sepia"], -# labels=["Concrete", "Insulation"], -# show_edges=False, -# ).show() +container.plot( + colors=["lightgrey", "sepia"], + labels=["Concrete", "Insulation"], + show_edges=False, +).show() # %% -# A top-level temperature field is defined on the whole construction with an initial -# temperature value of 10 °C, and separate fields are defined for each material. The -# surface heat transfer coefficients and ambient temperatures are defined for the -# internal and external boundaries. Thermal solid bodies are created for each material. +# A top-level temperature field is defined on the whole construction with an +# initial temperature value of 20 °C, and separate fields are defined for each +# material. Thermal solid bodies are created for each material. regions = [fem.RegionQuad(m) for m in container] fields = [fem.Field(r, dim=1).as_container() for r in regions] - -# top level temperature field mesh = container.stack() region = fem.RegionQuad(mesh) -temperature = fem.Field(region, dim=1, values=20.0) # initial temperature 10 °C +temperature = fem.Field(region, dim=1, values=20.0) # initial temperature 20 °C field = fem.FieldContainer([temperature]) materials = [] @@ -141,6 +143,7 @@ temperature=20.0, # °C ) +# %% # For the top and bottom surfaces, the detailed calculation approaches defined # in :class:`~felupe.thermal.SolidBodySurfaceConvection` and # :class:`~felupe.thermal.SolidBodySurfaceRadiation` are used for convection @@ -199,12 +202,8 @@ # thermal solid body, which returns the integrated surface heat flux for a given # boundary region and time step. # -# The convection coefficient values are calculated by ... -# -# The radiation coefficient values are calculated by ... -# -# All values are stored in the -# ``tstep_data`` dictionary, which is passed to the callback function as an argument. +# All values are stored in the ``tstep_data`` dictionary, which is passed to +# the callback function as an argument. def callback(stepnumber, substepnumber, substep, tstep_data): """Save mean surface heat flux at internal and external boundaries.""" @@ -242,13 +241,14 @@ def callback(stepnumber, substepnumber, substep, tstep_data): # %% -# The time step item is created with the thermal solid bodies. It must be located as the -# first item in the step to properly update the time step in the materials. The internal -# and external heat transfer item values are defined in the ramp, which specifies how -# their values change over time. Finally, a job is created with the step and the -# callback function, and evaluated with the top-level temperature field. A result file -# is created for visualization in Paraview, and the temperature field is saved as point- -# data in the result file. +# The time step item is created with the thermal solid bodies. It must be located +# as the first item in the step to properly update the time step in the materials. +# The side, top and bottom heat transfer item values as well as the pipe flux +# values are defined in the ramp, which specifies how their values change over +# time. Finally, a job is created with the step and the callback function, and +# evaluated with the top-level temperature field. A result file is created for +# visualization in Paraview, and the temperature field is saved as point-data +# in the result file. model_list = [*materials, side1_heat_transfer, side2_heat_transfer, top_convection, bottom_convection, top_radiation, bottom_radiation] @@ -307,11 +307,9 @@ def callback(stepnumber, substepnumber, substep, tstep_data): ax.legend() ax.set(xlim=(tmin, tmax), xlabel="time in h", ylabel=r"surface heat flux in W/m$^2$") -plt.savefig("_ex23a.png") - # %% -# Top and bottom surface heat transfer coefficients and pipe heat flux are -# plotted over time. +# Top and bottom convection and radiation surface heat transfer coefficients +# and pipe heat flux are plotted over time. fig, ax = plt.subplots() fig.subplots_adjust(right=0.75) @@ -319,9 +317,9 @@ def callback(stepnumber, substepnumber, substep, tstep_data): twin2 = ax.twinx() ax.set_xlabel("Time (s)") -ax.set_ylabel("Convection coefficient (W/(m2 K))") -twin1.set_ylabel("Temperature (°C)") -twin2.set_ylabel("Pipe heat flux (W/m2)") +ax.set_ylabel("Convection coefficient in W/(m$^2$ K)") +twin1.set_ylabel("Temperature in °C") +twin2.set_ylabel("Pipe heat flux in W/m$^2$") time_steps_h = time_steps / 3600 @@ -340,10 +338,7 @@ def callback(stepnumber, substepnumber, substep, tstep_data): twin2.spines['right'].set_position(('outward', 45)) -plt.savefig("_ex23b.png") - - # %% # A view on the temperature field at the end of the simulation period visualizes # the temperature distribution. -# field.plot("Field", scalar_bar_vertical=True).show() +field.plot("Field", scalar_bar_vertical=True).show() From e4f55d3ac01e73e5b4270c96e5b30e8cbc189d96 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 14 Jun 2026 10:41:42 +0200 Subject: [PATCH 70/71] Add SolidBodySurfaceConvection to test_thermal.py. Not sure how this is used, not tested. --- tests/test_thermal.py | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/tests/test_thermal.py b/tests/test_thermal.py index 814de6767..6f039efac 100644 --- a/tests/test_thermal.py +++ b/tests/test_thermal.py @@ -27,6 +27,10 @@ def test_thermal(): temperature = fem.Field(region, dim=1, values=20.0) field = fem.FieldContainer([temperature]) + region_right = fem.RegionQuadBoundary(mesh, mask=mesh.x == 1.0) + temperature_right = fem.Field(region_right, dim=1) + field_right = fem.FieldContainer([temperature_right]) + region_bottom = fem.RegionQuadBoundary(mesh, mask=mesh.y == 0.0) temperature_bottom = fem.Field(region_bottom, dim=1) field_bottom = fem.FieldContainer([temperature_bottom]) @@ -37,7 +41,7 @@ def test_thermal(): boundaries = fem.BoundaryDict( left=fem.Boundary(temperature, fx=0, value=20.0), - right=fem.Boundary(temperature, fx=1, value=20.0), + # right=fem.Boundary(temperature, fx=1, value=20.0), ) solid = fem.thermal.SolidBodyThermal( @@ -63,6 +67,12 @@ def test_thermal(): temperature=10.0, # °C ) + heat_convection = fem.thermal.SolidBodySurfaceConvection( + field=field_right, + convection_coefficient=7.69, # W/(m2 K) + temperature=10.0, # °C + ) + heat_radiation = fem.thermal.SolidBodySurfaceRadiation( field=field_top, emissivity=0.8, # dimensionless, between 0 and 1 @@ -80,6 +90,8 @@ def test_thermal(): heat_transfer.assemble.matrix(field) heat_flux.assemble.vector(field) heat_flux.assemble.matrix(field) + heat_convection.assemble.vector(field) + heat_convection.assemble.matrix(field) heat_radiation.assemble.vector(field) heat_radiation.assemble.matrix(field) @@ -87,11 +99,15 @@ def test_thermal(): heat_flux.assemble.vector(field) heat_flux.assemble.matrix(field) + heat_convection.time_step = 0.0 + heat_convection.assemble.vector(field) + heat_convection.assemble.matrix(field) + heat_radiation.time_step = 0.0 heat_radiation.assemble.vector(field) heat_radiation.assemble.matrix(field) - time = fem.thermal.TimeStep([solid, heat_transfer, heat_flux, heat_radiation]) + time = fem.thermal.TimeStep([solid, heat_transfer, heat_flux, heat_convection, heat_radiation]) table = fem.math.linsteps([0, 0, 1], num=2) table_emissivity = fem.math.linsteps([1, 1, 1], num=2) * 0.8 ramp = { @@ -104,7 +120,7 @@ def test_thermal(): heat_radiation["emissivity"]: table_emissivity, } step = fem.Step( - items=[time, solid, heat_transfer, heat_flux, heat_radiation], + items=[time, solid, heat_transfer, heat_flux, heat_convection, heat_radiation], ramp=ramp, boundaries=boundaries, ) From bad6edda4d6d3157d1ca4bec6092de1514c146e3 Mon Sep 17 00:00:00 2001 From: AGeissler Date: Sun, 14 Jun 2026 10:43:37 +0200 Subject: [PATCH 71/71] Remove convection tester. --- src/felupe/thermal/_sbsc_tester.py | 115 ----------------------------- 1 file changed, 115 deletions(-) delete mode 100644 src/felupe/thermal/_sbsc_tester.py diff --git a/src/felupe/thermal/_sbsc_tester.py b/src/felupe/thermal/_sbsc_tester.py deleted file mode 100644 index 588d48f45..000000000 --- a/src/felupe/thermal/_sbsc_tester.py +++ /dev/null @@ -1,115 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Tester version for _solid_body_surface_convection.py example case. -""" - -import numpy as np -import felupe as fem - -hc_fun = np.vectorize(fem.FreeConvection(5, 5).hc_fun) - - -mesh = fem.Rectangle(b=(1.0, 0.25), n=(11, 11)) # rectangle w/ 10x10 cells -region = fem.RegionQuad(mesh) -temperature = fem.Field(region, dim=1, values=30.0) -field = fem.FieldContainer([temperature]) - -region_convection = fem.RegionQuadBoundary(mesh, mask=mesh.y == 0.25) -temperature_convection = fem.Field(region_convection, dim=1) -field_convection = fem.FieldContainer([temperature_convection]) - -boundaries = fem.BoundaryDict( - bottom=fem.Boundary(temperature, fy=0, value=30.0), -) - -solid = fem.thermal.SolidBodyThermal( - field=field, - mass_density=1400.0, # kg / m^3 - specific_heat_capacity=1000.0, # J / (kg K) - time_step=720.0, # s - thermal_conductivity=1.0, # W / (m K) -) - -convection_function = fem.thermal.SolidBodySurfaceConvection( - field=field_convection, - convection_coefficient=hc_fun, - temperature=20.0, # °C -) - -def callback(stepnumber, substepnumber, substep, tstep_data): - """Save time step data at top (convective) boundary. - """ - tamb = list(ramp.values())[1][substepnumber] - ts = list(convection_function.field.extract(grad=False)[0])[0][0:1][0].mean() - heat_flux = solid.heat_flux_boundary - qc = heat_flux(region=region_convection) - tstep_data["tstep.s"].append(list(ramp.values())[0][substepnumber]) - tstep_data["qc_top.W.m-2"].append(qc) - tstep_data["tamb.degC"].append(tamb) - tstep_data["ts_top.degC"].append(ts) - tstep_data["hc_top.W.m-2.K-1"].append( - convection_function.results.convection_coefficient.mean()) - tstep_data["hc_fun_top.W.m-2.K-1"].append(hc_fun(ts, tamb).mean()) - tstep_data["hc_calc_top.W.m-2.K-1"].append(abs(qc/(ts-tamb))) - - -n_steps = 20 -time = fem.thermal.TimeStep([solid]) -table = fem.math.linsteps([0, 1], num=n_steps) -air_temperature = fem.math.linsteps([15, 25], num=n_steps) - -ramp = { - time: 18000 * table, # five hours - convection_function: air_temperature, -} - -step = fem.Step( - items=[time, solid, convection_function], ramp=ramp, boundaries=boundaries -) - -tstep_data = {"tstep.s": [], "tamb.degC": [], "ts_top.degC": [], - "hc_top.W.m-2.K-1": [], "qc_top.W.m-2": [], - "hc_fun_top.W.m-2.K-1": [], - "hc_calc_top.W.m-2.K-1": []} - -job = fem.Job(steps=[step], callback=callback, tstep_data=tstep_data).evaluate( - verbose=False -) - -# mesh.view( -# point_data={"Temperature 2 in °C": temperature.values} -# ).plot("Temperature 2 in °C").show() - -# Plot h_c / surface temp. / air temp. vs. time. -# https://matplotlib.org/stable/gallery/axes_grid1/parasite_simple.html -# or -# https://stackoverflow.com/questions/9103166/multiple-axis-in-matplotlib-with-different-scales -import matplotlib.pyplot as plt - -fig, ax = plt.subplots() -fig.subplots_adjust(right=0.75) - -twin1 = ax.twinx() -twin2 = ax.twinx() - -ax.set_xlabel("Time (s)") -ax.set_ylabel("Convection coefficient (W/(m2 K))") -twin1.set_ylabel("Temperature (°C)") -twin2.set_ylabel("Heat flux (W/m2)") - -p1 = ax.plot(tstep_data["tstep.s"], tstep_data["hc_top.W.m-2.K-1"], - label="hc", color='lightblue') -p2 = ax.plot(tstep_data["tstep.s"], tstep_data["hc_fun_top.W.m-2.K-1"], - label="hc_fun(ts_top - t_amb)", color='darkblue') -p3 = twin1.plot(tstep_data["tstep.s"], tstep_data["tamb.degC"], - label="t_amb", color='blue') -p4 = twin1.plot(tstep_data["tstep.s"], tstep_data["ts_top.degC"], - label="ts_top", color='red') -p5 = twin2.plot(tstep_data["tstep.s"], tstep_data["qc_top.W.m-2"], - label="qc_top", color='magenta') - -ax.legend(handles=p1+p2+p3+p4+p5, labelcolor="linecolor") - -twin2.spines['right'].set_position(('outward', 45)) - -plt.savefig("_test.png")