diff --git a/.gitignore b/.gitignore
index 80af874f..7bc3083a 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
diff --git a/docs/felupe/thermal.rst b/docs/felupe/thermal.rst
index 4bbec539..1a7a9842 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/examples/ex23_solid_body_thermal-tabs.py b/examples/ex23_solid_body_thermal-tabs.py
new file mode 100644
index 00000000..187dbc1f
--- /dev/null
+++ b/examples/ex23_solid_body_thermal-tabs.py
@@ -0,0 +1,344 @@
+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`
+
+   * 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 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
+model and geometry. The model is two-dimensional. The system is set up with two
+:class:`solids <felupe.thermal.SolidBodyThermal>`. The temperature boundary
+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 231 W/m2 and directly
+injected at the internal concrete surfaces (no pipe material is modelled).
+
+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
+
+import felupe as fem
+
+# %%
+# 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)
+
+# %%
+# One mesh per material is set up. If a material consists of multiple areas, these
+# are collected in a :class:`mesh container <felupe.MeshContainer>` 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.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.translate(0.02, axis=0),  # left
+        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,
+    decimals=6,
+).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.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
+insulation = fem.MeshContainer(
+    [
+        insulation_1,
+        insulation_1.translate(1.0, axis=0),
+    ],
+    merge=True,
+    decimals=6,
+).stack()
+
+container = fem.MeshContainer([concrete, insulation], merge=True, decimals=6)
+
+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 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]
+mesh = container.stack()
+region = fem.RegionQuad(mesh)
+temperature = fem.Field(region, dim=1, values=20.0)  # initial temperature 20 °C
+field = fem.FieldContainer([temperature])
+
+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,
+        )
+    )
+
+# %%
+# 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])
+
+side_region2 = fem.RegionQuadBoundary(mesh, mask=mesh.x == mesh.x.max())
+side_temperature2 = fem.Field(side_region2, dim=1)
+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)
+bottom_field = fem.FieldContainer([bottom_temperature])
+
+top_region = fem.RegionQuadBoundary(mesh, mask=mesh.y == mesh.y.max())
+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)
+    temperature=20.0,  # °C
+)
+side2_heat_transfer = fem.thermal.SolidBodySurfaceHeatTransfer(
+    field=side_field2,
+    coefficient=7.69,  # W/(m^2 K)
+    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_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_top,  # W/(m^2 K)
+    temperature=20.0,  # °C
+)
+bottom_convection = fem.thermal.SolidBodySurfaceConvection(
+    field=bottom_field,
+    convection_coefficient=hc_bottom,  # W/(m^2 K)
+    temperature=20.0,  # °C
+)
+
+top_radiation = fem.thermal.SolidBodySurfaceRadiation(
+    field=top_field,
+    emissivity=0.9,
+    temperature=20.0,  # °C
+)
+
+bottom_radiation = fem.thermal.SolidBodySurfaceRadiation(
+    field=bottom_field,
+    emissivity=0.9,
+    temperature=20.0,  # °C
+)
+
+# %%
+# 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]])
+
+pipe_region = []
+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)]))
+    pipe_flux.append(fem.thermal.SolidBodyHeatFlux(
+        field=pipe_field[idx],
+        heat_flux=-231.25,  # W / m^2, 74/(4*4*0.02)
+    ))
+
+# %%
+# 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.
+#
+# 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."""
+
+    heat_flux = materials[0].heat_flux_boundary
+    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_)
+    tstep_data["pipes"].append(pflux)
+
+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_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 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]
+
+time = fem.thermal.TimeStep(model_list)
+ramp = {
+    time: time_steps,
+    side1_heat_transfer: t_air,
+    side2_heat_transfer: t_air,
+    top_convection: t_air,
+    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,
+    ramp=ramp,
+)
+
+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, tstep_data=tstep_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,
+)
+
+# %%
+# Top and bottom surface heat flux values are plotted over time.
+#
+# .. note::
+#
+#    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")
+
+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$")
+
+# %%
+# 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)
+
+twin1 = ax.twinx()
+twin2 = ax.twinx()
+
+ax.set_xlabel("Time (s)")
+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
+
+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))
+
+# %%
+# 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/__init__.py b/src/felupe/__init__.py
index 1a6597ed..477cb49a 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 5f6e1af5..a685c559 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 00000000..17b0509f
--- /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 00000000..35907b20
--- /dev/null
+++ b/src/felupe/constitution/heat_transfer/_free_convection.py
@@ -0,0 +1,194 @@
+# -*- 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 <http://www.gnu.org/licenses/>.
+"""
+
+import math
+
+from pyfluids import HumidAir, InputHumidAir
+from scipy.constants import g
+
+
+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 %)
+        Relative humidity of air in (%).
+
+    Notes
+    -----
+    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`),
+    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
+    given in :eq:`l-horiz-plate` (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
+    ----------
+    ..  [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, 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
+
+        # Characteristic length for horizontal plate.
+        self.length = self.plate_width*self.plate_length/\
+                           (2*(self.plate_width+self.plate_length))
+
+    def _pyfluids_units(self):
+        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(self, 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.p_abs),
+                  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
+
+    def _nusselt_horizontal(self, 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(self, ts, tamb):
+        """
+        Calculate convection coefficient for horizontal plate.
+        """
+        dtk, rhf = self._pyfluids_units()
+        tm_c = (ts + tamb)/2
+
+        air = HumidAir().with_state(
+                  InputHumidAir.pressure(self.p_abs),
+                  InputHumidAir.temperature(tm_c + dtk),
+                  InputHumidAir.relative_humidity(self.rh/rhf),
+              )
+
+        ra = self._rayleigh(ts_c=ts, ti_c=tamb, length_=self.length)
+        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)
diff --git a/src/felupe/thermal/__init__.py b/src/felupe/thermal/__init__.py
index 30eb6e1f..9ceede79 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
new file mode 100644
index 00000000..ad2a630c
--- /dev/null
+++ b/src/felupe/thermal/_solidbody_surface_convection.py
@@ -0,0 +1,305 @@
+# -*- 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 <http://www.gnu.org/licenses/>.
+"""
+import numpy as np
+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 or callable
+        Convection heat transfer coefficient :math:`h_c` in W/(m^2 K). A
+        callable requires the parameters 'surface temperature' and 'ambient
+        temperature'.
+    temperature : float
+        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 convective 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` [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
+    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
+
+       q_c = h_c\,\left(\theta_s - \theta_\infty\right)
+
+    .. math::
+       :label: example-horizontal-plate
+
+       L &= \frac{A}{P}
+
+       T_m &= 0.5 \left(T_s + T_\infty\right)
+
+       \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}
+
+       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 &= \frac{Nu\,\lambda_\text{air}}{L}
+
+
+    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:
+        :force_static:
+
+        >>> import math
+        >>> import numpy as np
+        >>>
+        >>> 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
+        ...     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 (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.33)
+        ...     else:
+        ...         nu = 0.15 * math.pow(1E11, 0.33)
+        ...     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::
+        :context:
+        :force_static:
+
+        >>> 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)
+        >>> 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)
+        ... )
+
+    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: 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", 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,
+        ...     convection_coefficient=hc_fun,
+        ...     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_function: air_temperature,
+        ... }
+
+    ... and run.
+
+    ..  pyvista-plot::
+        :context:
+        :force_static:
+
+        >>> step = fem.Step(
+        ...     items=[time, solid, convection_function], ramp=ramp, boundaries=boundaries
+        ... )
+        >>> job = fem.Job(steps=[step]).evaluate(verbose=False)
+        >>>
+        >>> mesh.view(
+        ...     point_data={"Temperature hc_fun() in °C": temperature.values}
+        ... ).plot("Temperature hc_fun() in °C", off_screen=True).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
+    --------
+    felupe.thermal.TimeStep : A time step item.
+    felupe.thermal.SolidBodyThermal : A thermal solid body for heat conduction.
+
+    """
+
+    def __init__(self, field, convection_coefficient, temperature):
+        self.field = field
+        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(
+                    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
+        )
+
+    def __getitem__(self, key):
+        return UpdateItem(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_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
+                )
+        else:
+            self.results.convection_coefficient = self.convection_coefficient
+
+    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.convection_coefficient
+            * (temperature - self.results.temperature)
+        ]
+
+        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.convection_coefficient
+            * 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
diff --git a/src/felupe/thermal/_solidbody_surface_heat_transfer.py b/src/felupe/thermal/_solidbody_surface_heat_transfer.py
index e46a85e9..70d9b61d 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 ca946fe9..97338b69 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:
 
@@ -112,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
 
@@ -160,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)
         ]
 
diff --git a/src/felupe/thermal/_solidbody_thermal.py b/src/felupe/thermal/_solidbody_thermal.py
index 0455d908..edef6d85 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}
 
@@ -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.
 
diff --git a/tests/test_thermal.py b/tests/test_thermal.py
index 814de676..6f039efa 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,
     )