Introduction

General goal

The design of a camping van involves an estimation of the heat losses that will lead to the subsequent design of both the insulation and the heating system.

This study will aim at evaluating the average heat flux through the van surfaces in steady mode.

Tools of the trade

For that matter, we will use mostly Numpy, the scientific computing module for vector and matrice algebra in Python.

As a wrapper for Numpy, we will use the package Pandas which offers a nice way to index rows and columns of Numpy arrays,like an abstraction of a spreadsheet, with SQL-like tools to get, mask and insert data.

Surfaces description

All dimensions in meter in the code, and in mm in the picture, the main dimension is the dimension in axis of the air flow (when the wind is blowing on the van). We use a 2018 Mercedes Sprinter van as a reference.

sprinter © Mercedes-Benz

Notice that the frontdoors are integrated to the sides en the model, and the area below the windchill is ignored. We consider the front windows, 2 roof vents and a side vent on the bottom part of one side as an air inlet (see the section "taking account of the vents").

In [1]:
import numpy as np
import pandas as pd
import plotly.graph_objs as go
from plotly.plotly import plot, iplot
import numba
import colorsys
#init_notebook_mode(connected=False)
In [2]:
# Create the dataframe
surfaces = pd.DataFrame(
    columns=[
        #Geometry
        "Main dimension", "Secondary dimension", "Thickness", "blind",
        # Nature and properties of the flow - internal side of the body
        "flow_in", "surface_in", "L_in",
        # Nature and properties of the flow - external side of the body
        "flow_out", "surface_out", "L_out"
    ],
    index=["roof", "sides", "floor", "windchill", "backdoors", "vents", "windows"])

# thickness of the sheet metal in m : wild guess
sheet = 0.0007

# width of the average surface between external and internal sides
w = 1.8

# length of the sides
L = 4.8

# caracteristic length for horizontal plates (roof and floor)
L_c = (w * L) / (2 * w + 2 * L)

# height of the sides
h = 1.900

surfaces.loc["roof"] = [
    w, L, sheet, True,
    "free", "horizontal_cold_bottom", L_c,
    "external", "parallel", L]

surfaces.loc["sides"] = [
    L * 2 + 0.75 * 2, h, sheet, True,  # 2 sides at once
    "free", "vertical", h,
    "external", "parallel", L]

surfaces.loc["floor"] = [
    w, 5.5, 0.005, True,
    "free", "horizontal_cold_top", L_c,
    "external", "parallel", L]

surfaces.loc["windchill"] = [
    1, 1.79, 0.002, False,
    "free", "vertical", 1,
    "external", "perpendicular", 1.79]

surfaces.loc["backdoors"] = [
    1.84, 1.56, sheet, True,  # 2 doors at once
    "free", "vertical", h,
    "external", "perpendicular", 1.56]

surfaces.loc["vents"] = [
    0.8, 0.8, 0.002, False,  # 2 roof vents at once
    "free", "horizontal_cold_bottom", L_c,
    "external", "parallel", L]

surfaces.loc["windows"] = [
    0.65, 1.79 + 1.6, 0.002, False,  # 2 windows at once
    "free", "vertical", h,
    "external", "parallel", L]

surfaces["Area"] = surfaces["Main dimension"] * surfaces["Secondary dimension"]
surfaces["Area"].sum()
Out[2]:
47.133900000000004

Insulation description

The insulation will be 3 cm of polyurethan foam and 3 mm of fiberglass/epoxy composite on blind panels. On windows and 1 cm of multilayered sheets on windows. Vents are supposed uncovered.

In [4]:
# Thermal Conductivity coeff in W/m/K
k_steel = 46
k_wool = 0.026 # polyurethan actually
k_fiber = 0.03
k_composite = 0.039
k_glass = 1.4

# Thicknesses in m
th_wool = 0.03
th_composite = 0.01
th_fiber = 0.003
In [5]:
# Compute the resistance of the compound walls in m²K/W
surfaces["R"] = 0
surfaces.loc[surfaces["blind"] == 1, "R"] = surfaces.loc[surfaces["blind"] == 1, "Thickness"] / k_steel + th_wool / k_wool + th_fiber / k_fiber
surfaces.loc[surfaces["blind"] == 0, "R"] = surfaces.loc[surfaces["blind"] == 0, "Thickness"] / k_glass + th_composite / k_composite

surfaces.loc["vents", "R"] = surfaces.loc["vents", "Thickness"] / k_glass
In [6]:
surfaces
Out[6]:
Main dimension Secondary dimension Thickness blind flow_in surface_in L_in flow_out surface_out L_out Area R
roof 1.8 4.8 0.0007 True free horizontal_cold_bottom 0.654545 external parallel 4.8 8.64 1.253861
sides 11.1 1.9 0.0007 True free vertical 1.9 external parallel 4.8 21.09 1.253861
floor 1.8 5.5 0.005 True free horizontal_cold_top 0.654545 external parallel 4.8 9.9 1.253955
windchill 1 1.79 0.002 False free vertical 1 external perpendicular 1.79 1.79 0.257839
backdoors 1.84 1.56 0.0007 True free vertical 1.9 external perpendicular 1.56 2.8704 1.253861
vents 0.8 0.8 0.002 False free horizontal_cold_bottom 0.654545 external parallel 4.8 0.64 0.001429
windows 0.65 3.39 0.002 False free vertical 1.9 external parallel 4.8 2.2035 0.257839

Convection calculus

Hypotheses

The conductive only models is a worse case scenario that doesn't describe the fact that the sourrounding air is not a perfect conductor. Adding the convection will take account of the steady air (which gives less pessimistic results) as well as the wind.

The worse case scenario is when the wind is flowing along the length of the van because the ceiling and both side are "wet" and the air flow becomes turbulent at some point. If the wind was flowing across the van, the chances are the flow would stay laminar everywhere, leading to a decreased heat transfer.

We will asume a steady regime where the inside temperature is forced by a heating device (Dirichlet boundary condition).

Heat transfer equations

We assume the van is parallel to a forced flow at a certain speed on the rooftop and the sides, and windchill and rear doors are considered stopping points where the average flow speed is zero, but turbulences happen. We use the 1D model of the equivalent resistances for the composite wall because the thickhess of the wall is neglictible regarding its length:

Flow

where :

  • $R = \sum_i \frac{L_i}{k_i} = \frac{L}{k}$, with $R$ the thermal resistance of the solid wall, already computed at the previous section (insulation), in m²K/W,
  • $R_{in} = \frac{1}{h_{in}}$, with $h_{in}$ the convection coefficient of the internal air volume in W/m²/K
  • $R_{out} = \frac{1}{h_{out}}$, with $h_{out}$ the convection coefficient of the extrenal air flow in W/m²/K

N.B : Recalling our previous analysis, we actually assumed that $T_{in} = T_{s, i}$ and $T_\infty = T_{s, o}$

Finally:

$$ Q = A \frac{T_i - T_{\infty}}{R + R_i + R_o} $$

where $Q$ is the heat transfer in W over a surface which area is $A$ in m².

Now the problem is to compute the convection coefficents $h$.

Convection factors equations

The reference book of this section is Fundamentals on Heat and Mass Transfer, by Bergman, Lavine, Incropera & Dewitt, 7th edition, at Wiley.

$$ \bar h = \frac{\overline{Nu_L} k}{L}$$

where :

  • $\bar h$ is the average convection coefficient on a surface where the flow is laminar on a portion and turbulent on the rest,
  • $L$ is the characteristic length of the 1D problem,
  • $k$ is the thermal conduction coefficient of the flow

We use an average coefficient since we don't know the critical length where the flow becomes turbulent, and we don't really care for that matter. $k$ is the

$\overline{Nu_L}$ is the Nusselt adimensional parameter given by the following equations :

For the inside control volume, we use the external free convection model assuming quiescent air :

  • for vertical panels, we use the equation 9.26 of the reference book : $$\overline{Nu_L} = \left(0.825 + \frac{0.387 Ra_L^{1/6}}{[1 + (0.492/Pr)^{9/16}]^{8/27}}\right)^2$$
  • for the roof, we use the equation 9.32 : $$\overline{Nu_L} = 0.52 Ra_L^{1/5}$$
  • for the floor, we use the equations :
    • 9.30 : $$\overline{Nu_L} = 0.54 Ra_L^{1/4}$$ if $Ra_L < 10^7$ and $Pr > 0.7$
    • 9.31 : $$\overline{Nu_L} = 0.15 Ra_L^{1/3}$$ if $Ra_L > 10^7$

For the outside flow:

  • for the back doors and the windchill, we use equation 7.53 of the reference book : $$\overline{Nu_L} = 0.3 + [0.62 Re_D^{1/2} Pr^{1/3} [1 + (0.4 / Pr)^{2/3}]^{-1/4}][1 + (Re_d / 282000)^{5/8}]^{4/5}$$
  • for the side panels and roof, we use the equation 7.38 : $$\overline{Nu_L} = (0.037 Re_L ^{4/5} - (0.037 Re_c^{4/5} - 0.664 Re_c^{1/2}) Pr^{1/3}$$

with :

  • $Ra_L$ the Rayleigh adimensional number such as $Ra_L = Gr_L Pr = Pr \frac{g \beta (T_s - T_\infty) L^3} {\nu \alpha}$
  • $Pr$ the Prandtl adimensional number,
  • $Gr = \frac{g \beta (T_s - T_\infty) L^3} {\nu}$ the Grashof adimensional number,
  • $Re_L$ the Reynolds adimensional number such as $Re_L = VL/\nu$ ,
  • $Re_c = 50000$ the critical Reynolds number representing the laminar/turbulent transition,
  • $V$ the velocity of the flow,
  • $g$ the gravitational constant,
  • $T_s$ the temperature of the solid surface, either $T_{s, i}$ or $T_{s, o}$ to be consistent with the previous section,
  • $T_\infty$ the temperature of the flow,
  • $L$ the characteristic length,
  • $\nu$ the kinematic viscosity of the fluid
  • $\beta$ the expansion coefficient in 1/K
  • $\alpha$ the diffusion coefficient in m²/s

For the backdoors, we take the model of the cylinder in cross flow, even if we know it is valid for shallow objects. We use the same model for the windchill, which is in the drag zone. Again, we know this is an inaccurate approximation, but it is the most accurate way to stay in 1D while avoiding FEM.

The parameter $L$ should change depending the nature of the flow :

  • in an external flow, it is the length of the surface in the axis of the flow,
  • in an internal quiescent air volume, it is the vertical dimension for vertical panels and the ratio between the area and the perimiter $A_s / P$ for horizontal panels.

The parameters have to be evaluated at the film temperature $T_f = (T_s + T_\infty) / 2$.

Its is useful to remember that these equations are correlations of the solutions of the heat transfer differential equations solved for classic problems in standard conditions. Not only are they valid inside a given range of $Ra$, $Re$ and $Pr$ values, but they lead to an error between 2 % (laminar models) to 20 % (turbulent models) compared to the numeric solution of the differential equations. So we will add 10 % on top of $Q$ to take this error into consideration.

A more accurate study would have to use FEM, but what we want here is more an idea, and as the surfaces are already majorized, we think the study is pessimistic enough for its purpose.

Solving

We notice that as $Ra_L$ depends on $T_s$ on both sides, which are unknown (isothermal formulation), so we need an iterative solver to refine a preliminary assumption :

  1. assume $T_{s, i} = T_{s, o} = (T_{in} + T_{out})/2$
  2. compute $Gr_L$ then $Ra_L$ then $\overline{Nu_L}$ on both sides of the wall
  3. compute $\bar h$ on both sides of the wall
  4. compute $R_i$, $R_o$ and finally $Q$
  5. using $Q$ and $\bar h$, update $T_{s, i}$ and $T_{s, o}$ with the Fourier law $Q = \bar h A (T_s - T_\infty)$
  6. measure convergence with the error between $Q$ for the composite wall and $Q$ for the conduction part in L2 norm
  7. goto 2. and iterate until convergence
  8. output $Q$

Running the computations, we first notice that the equations for external forced flow are inaccurate at low wind speed. As a refinement, we then compute the outside $h$ with the free convection equations, then with the external forced convection equations, and chose the higher value.

Properties

Let's begin by turning the properties tables of air into interpolated functions. From the table A.4 of the reference book. Notice that, as this functions are basic and may re ru-used in further work, we take some time to code them in a library style with debugging symbols.

What we do essentially is interpolated the tables with polynoms, then store the coefficients hard set in the code, and define debug functions to ensure the consistency between the hard-set values and the interpolated values, as well as the plots to ensure the data consistency. We raise Python exceptions in case of a problem to simplify the debugging. These functions can either take vectors or scalar as inputs

In [7]:
# Show the debugging symbols instead of outputs
DEBUG = False
GRAPH = False

# Global variable !!!!
temperatures = np.array([100, 150, 200, 250, 300, 350, 400, 450, 500])  # K

@numba.jit(cache=True)
def debug_print(values, result):
    """Print a graph of interpolated values versus tabulated to ensure the quality of the output"""

    if GRAPH:
        data = [
            go.Scatter(x=temperatures, y=values, name="Tabulated"),
            go.Scatter(x=temperatures, y=result, name="Interpolated")
        ]

        iplot(data, show_link=False, filename='test')


@numba.jit(cache=True)
def polyfit(values, T, w=np.ones_like(temperatures), order=4):
    """Least-squares weighted polynomial fit with a 5th order"""

    coef = np.polyfit(T, values, order, w=w)
    result = np.polyval(coef, T)

    print("Polynomial fit coefficients:",
          np.array2string(coef, precision=12, separator=", "))

    result = np.polyval(coef, T)

    return result, coef


@numba.jit(cache=True)
def input_check(T):
    """Check if the temperature is in the tabulated range"""

    if np.min(T) < temperatures.min() or np.max(T) > temperatures.max():
        raise ValueError(
            "The temperature is of of range. T should be between %i and %i K" %
            (temperatures.min(), temperatures.max()))


@numba.jit(cache=True)
def Pr(T):
    """Prandtl number at asmospheric pressure - 28.97 kg/kmol
        * T : temperature in K
    """

    input_check(T)
    coef = [
        1.878350289902e-12, -2.983250069892e-09, 2.294193742479e-06,
        -9.993034435813e-04, 8.656918371440e-01
    ]
    result = np.polyval(coef, T)

    if DEBUG:
        # Do the actual fitting
        values = np.array(
            [0.786, 0.758, 0.737, 0.720, 0.707, 0.700, 0.690, 0.686, 0.684])

        # Remove the point at 350 which looks like an outliner
        weights = np.ones_like(temperatures)
        weights[5] = 0

        debug_result, coef_updated = polyfit(values, temperatures, w=weights)
        debug_print(values, debug_result)

        if not np.allclose(coef, coef_updated):
            raise RuntimeError(
                "The hard-coded fitting parameters do not match the updated fitting"
            )

        print("Pr = ", result)

    return result


@numba.jit(cache=True)
def nu(T):
    """Kinematic viscosity of air at asmospheric pressure - 28.97 kg/kmol, in m²/s
        * T : temperature in K
    """

    input_check(T)
    coef = [
        5.519813519814e-17, -1.380761460762e-13, 2.047104895105e-10,
        3.336099456098e-09, -2.466666666666e-07
    ]
    result = np.polyval(coef, T)

    if DEBUG:
        # Do the actual fitting
        values = np.array(
            [2, 4.426, 7.590, 11.44, 15.89, 20.92, 26.41, 32.39, 38.79]) * 1e-6
        debug_result, coef_updated = polyfit(values, temperatures)
        debug_print(values, debug_result)

        if not np.allclose(coef, coef_updated):
            raise RuntimeError(
                "The hard-coded fitting parameters do not match the updated fitting"
            )

        print("nu = ", result)

    return result


@numba.jit(cache=True)
def k(T):
    """Thermal conduction coefficient of air at asmospheric pressure - 28.97 kg/kmol, in W/m/K
        * T : temperature in K
    """

    input_check(T)
    coef = [
        3.076923076924e-14, -2.924630924632e-11, -2.250116550116e-08,
        9.619432789433e-05, -3.333333333325e-05
    ]
    result = np.polyval(coef, T)

    if DEBUG:
        # Do the actual fitting
        values = np.array(
            [9.34, 13.8, 18.1, 22.3, 26.3, 30.0, 33.8, 37.3, 40.7]) * 1e-3
        debug_result, coef_updated = polyfit(values, temperatures)
        debug_print(values, debug_result)

        if not np.allclose(coef, coef_updated):
            raise RuntimeError(
                "The hard-coded fitting parameters do not match the updated fitting"
            )

        print("k = ", result)

    return result


@numba.jit(cache=True)
def alpha(T):
    """Thermal diffusivity coefficient of air at asmospheric pressure - 28.97 kg/kmol, in m²/s
        * T : temperature in K
    """

    input_check(T)
    coef = [
        -8.764568764569e-17, -3.812483812483e-14, 2.638554778555e-10,
        2.580160580161e-09, -3.111111111111e-07
    ]
    result = np.polyval(coef, T)

    if DEBUG:
        # Do the actual fitting
        values = np.array(
            [2.54, 5.84, 10.3, 15.9, 22.5, 29.9, 38.3, 47.2, 56.7]) * 1e-6
        debug_result, coef_updated = polyfit(values, temperatures)
        debug_print(values, debug_result)

        if not np.allclose(coef, coef_updated):
            raise RuntimeError(
                "The hard-coded fitting parameters do not match the updated fitting"
            )

        print("alpha = ", result)

    return result


## Reminder
# temperatures = np.array([100, 150, 200, 250, 300, 350, 400, 450, 500])  # K


@numba.jit(cache=True)
def enthalpy(T):
    """Enthalpy of air at atmospheric pressure - 28.97 kg/kmol, in kJ/kg
        * T : temperature in K
        
        These values come from Çengel & Boles: Thermodynamics, an engineer approach, Table A.17 
    """

    input_check(T)
    coef = [
        7.878787878779e-11, 3.030303030430e-09, -1.783333333341e-05,
        1.005474242424e+00, -5.614285714302e-01
    ]

    result = np.polyval(coef, T)

    if DEBUG:
        # Do the actual fitting

        ## WARNING : no properties below 200 K
        values = np.array(
            [0, 0, 199.97, 250.05, 300.19, 350.49, 400.98, 451.8, 503.02])

        # Remove properties below 200 K from fitting
        weights = np.array([0, 0, 1, 1, 1, 1, 1, 1, 1])

        debug_result, coef_updated = polyfit(values, temperatures, w=weights)
        debug_print(values, debug_result)

        if not np.allclose(coef, coef_updated):
            raise RuntimeError(
                "The hard-coded fitting parameters do not match the updated fitting"
            )

        print("enthalpy = ", result)

    return result


@numba.jit(cache=True)
def rho(T):
    """Density of air at asmospheric pressure - 28.97 kg/kmol, in kg/m³
        * T : temperature in K
    """

    input_check(T)
    coef = [
        -1.383487179487e-12, 2.440946386946e-09, -1.708066317016e-06,
        6.019370279720e-04, -1.115476571096e-01, 1.016748333333e+01
    ]
    result = np.polyval(coef, T)

    if DEBUG:
        # Do the actual fitting
        values = np.array([
            3.5562, 2.3364, 1.7458, 1.3947, 1.1614, 0.9950, 0.8711, 0.7740,
            0.6964
        ])
        debug_result, coef_updated = polyfit(values, temperatures, order=5)
        debug_print(values, debug_result)

        if not np.allclose(coef, coef_updated):
            raise RuntimeError(
                "The hard-coded fitting parameters do not match the updated fitting"
            )

        print("alpha = ", result)

    return result


# Display graphs to debug fittings
if DEBUG:
    alpha(temperatures)
    Pr(temperatures)
    nu(temperatures)
    k(temperatures)
    enthalpy(temperatures)
    rho(temperatures)
    
    
R = 287.058 # J/kg/K - ideal gases constant for air

def zero_centered_heatmap(dataset):
    """Draw a red/blue color scale where 0 is grey"""
    
    MAX = np.amax(dataset)
    MIN = np.amin(dataset)
    RANGE = abs(MAX) + abs(MIN)
    ZERO = abs(MIN) / RANGE

    max_saturation = abs(MAX) / RANGE
    min_saturation = abs(MIN) / RANGE

    maxR, maxG, maxB = colorsys.hsv_to_rgb(0.03, max_saturation + 0.2,
                                           1 - max_saturation / 2)
    minR, minG, minB = colorsys.hsv_to_rgb(0.60, min_saturation + 0.2,
                                           1 - min_saturation / 2)

    zeroR, zeroG, zeroB = colorsys.hsv_to_rgb(0.57 * min_saturation, 0. * 1,
                                              0.9)

    return [[0.0,
             'rgb(%i, %i, %i)' % (minR * 255, minG * 255, minB * 255)], 
            [ZERO,
             'rgb(%i, %i, %i)' % (zeroR * 255, zeroG * 255, zeroB * 255)],
            [1.0,
             'rgb(%i, %i, %i)' % (maxR * 255, maxG * 255, maxB * 255)]]

Functions

All the variables that are not tabulated.

In [8]:
@numba.jit(cache=True)
def beta(T):
    """Thermal expansion coefficient of air at asmospheric pressure - 28.97 kg/kmol, in 1/T
        * T : temperature in K
    """

    # Ideal gas model
    result = 1 / T
    if DEBUG:
        print("beta = ", result)

    return result


@numba.jit(cache=True)
def Gr(Ts, Tinf, L):
    """Grashof adimensional number
        * Ts, Tinf in K
        * L in m
    """

    T_mean = (Ts + Tinf) / 2

    result = 9.81 * beta(T_mean) * np.abs(Ts - Tinf) * L**3 / nu(T_mean)

    if DEBUG:
        print("Gr = ", result)

    return result


@numba.jit(cache=True)
def Re(v, L, T):
    """Reynolds adimensional number
        * v in m/s
        * L in m
        * T in K
    """

    result = v * L / nu(T)

    if DEBUG:
        print("Re = ", result)

    return result


@numba.jit(cache=True)
def Ra(Ts, Tinf, L):
    """Rayleigh adimensional number
        * Ts, Tinf in K
        * L in m
    """

    T_mean = (Ts + Tinf) / 2

    result = Gr(Ts, Tinf, L) / alpha(T_mean)

    if DEBUG:
        print("Ra = ", result)

    return result


Re_c = 5e5  # Critical Reynolds number - limit turbulent/laminar


def Nu(Ts, Tinf, L, v, case, surface):
    """Average Nusselt adimensional number
    """

    def free(T_mean, Ra_L, Pr_L, Re_L, surface):
        """Free convection"""

        def vertical(T_mean, Ra_L, Pr_L, Re_L):
            # Vertical plate - L is the height

            if Ra_L <= 2e9:
                return 0.68 + (0.67 * Ra_L**(1 / 4)) / (1 + (0.492 / Pr_L) **
                                                        (9 / 16))**(4 / 9)
            elif Ra_L == 0:
                return 1e-15  # Trick to avoid division by 0
            else:
                return (0.825 + (0.387 * Ra_L**(1 / 6)) /
                        (1 + (0.492 / Pr_L)**(9 / 16))**(8 / 27))**2

        def horizontal_hot_top(T_mean, Ra_L, Pr_L, Re_L):
            # Horizontal plate with the hot part on top - L is the ratio Area / Perimeter

            if Ra_L >= 5e3 and Ra_L <= 1e7 and Pr_L >= 0.7:
                return 0.54 * Ra_L**(1 / 4)
            elif Ra_L >= 1e7 and Ra_L <= 2e11:
                return 0.15 * Ra_L**(1 / 3)
            elif Ra_L == 0:
                return 1e-15  # Trick to avoid division by 0
            else:
                return np.nan

        def horizontal_hot_bottom(T_mean, Ra_L, Pr_L, Re_L):
            # Horizontal plate with the hot part on bottom - L is the ratio Area / Perimeter

            if Ra_L >= 1e4 and Ra_L <= 2e9 and Pr_L >= 0.7:
                return 0.52 * Ra_L**(1 / 5)
            elif Ra_L == 0:
                return 1e-15
            else:
                return np.nan

        surfaces = {
            "vertical": vertical,
            "horizontal_hot_top": horizontal_hot_top,
            "horizontal_hot_bottom": horizontal_hot_bottom,
            "horizontal_cold_bottom": horizontal_hot_top,
            "horizontal_cold_top": horizontal_hot_bottom,
        }

        return surfaces[surface](T_mean, Ra_L, Pr_L, Re_L)

    def external(T_mean, Ra_L, Pr_L, Re_L, surface):
        """Convection on a surface in an external flow"""

        def parallel(T_mean, Ra_L, Pr_L, Re_L):
            # Plate in parallel flow - L is the lenght

            if Pr_L >= 0.6 and Re_L < Re_c:
                return 0.664 * Re_L**(1 / 2) * Pr_L**(1 / 3)
            elif Pr_L >= 0.6 and Pr_L <= 60 and Re_L >= Re_c and Re_L <= 2e8:
                A = 0.037 * Re_c**(4 / 5) - 0.664 * Re_c**(1 / 2)
                return (0.037 * Re_L**(4 / 5) - A) * Pr_L**(1 / 3)
            else:
                return np.nan

        def perpendicular(T_mean, Ra_L, Pr_L, Re_L):
            # Cylinder in cross-flow - L is the diameter

            if Pr_L >= 0.2:
                return 0.3 + (0.62 * Re_L**(1/2) * Pr_L**(1/3) * (1 + (0.4/Pr_L)**(2/3))**(-1/4)) * (1 + (Re_L/282000)**(5/8))**(4/5)
            else:
                return np.nan

        surfaces = {"parallel": parallel, "perpendicular": perpendicular}

        return surfaces[surface](T_mean, Ra_L, Pr_L, Re_L)

    cases = {
        "free": free,  # free convection
        "external": external,  # external flow
    }

    T_mean = (Ts + Tinf) / 2
    Ra_L = Ra(Ts, Tinf, L)
    Pr_L = Pr(T_mean)
    Re_L = Re(v, L, T_mean)
    result = pd.Series(index=Ts.index)

    for i in range(Ts.size):
        result[i] = cases[case[i]](T_mean[i], Ra_L[i], Pr_L[i], Re_L[i],
                                   surface[i])

    if DEBUG:
        print("Nu = ", result)

    return result


@numba.jit(cache=True)
def h(Ts, Tinf, L, v, case, surface):
    """Average convection coefficient number in W/m²/K"""
    T_mean = (Ts + Tinf) / 2

    result = Nu(Ts, Tinf, L, v, case, surface) * k(T_mean) / L

    if DEBUG:
        print("h = ", result)

    return result


@numba.jit(cache=True)
def deg2K(T):
    """Convert °C to K"""
    return T + 273.16


@numba.jit(cache=True)
def HeatFlux(Tin, Tout, R, A):
    result = (Tout - Tin) * A / R

    if DEBUG:
        print("Q = ", result)

    return result


@numba.jit(cache=True)
def SurfaceTemp(Tinf, Q, R, A):
    result = -Q * R / A + Tinf

    if DEBUG:
        print("T_s = ", result)

    return result


def inverse_surface(surface):
    """Map the nature of the outside surface condition knowing the inside surface condition"""
    
    surfaces = {
        "vertical": "vertical",
        "horizontal_hot_top": "horizontal_cold_bottom",
        "horizontal_hot_bottom": "horizontal_cold_top",
        "horizontal_cold_bottom": "horizontal_hot_top",
        "horizontal_cold_top": "horizontal_hot_bottom",
    }
    
    return surfaces[surface]

v_inverse_surface = np.vectorize(inverse_surface)


def solver(T_in, T_out, v, surfaces):

    T_mean = (T_in + T_out) / 2

    # Temperature vector
    surfaces["T_in"] = T_in
    surfaces["T_si"] = T_mean
    surfaces["T_so"] = T_mean
    surfaces["T_out"] = T_out

    # Iterations
    error = 1
    error_new = 1
    iterations = 0

    while (1 - np.abs(error_new - error) /
           (error + error_new)) >= 0.001 and iterations < 50:

        error = error_new

        # h
        surfaces["h_in"] = h(surfaces["T_in"], surfaces["T_si"],
                             surfaces["L_in"], 0, surfaces["flow_in"],
                             surfaces["surface_in"])
        
        ## Compute the outside h as natural convection
        h_out_quiescent = h(surfaces["T_so"], surfaces["T_out"],
                              surfaces["L_in"], v, surfaces["flow_in"],
                              v_inverse_surface(surfaces["surface_in"]))
        
        ## Compute the outside h as a forced flow
        h_out_flow = h(surfaces["T_so"], surfaces["T_out"],
                              surfaces["L_out"], v, surfaces["flow_out"],
                              surfaces["surface_out"])
            
        ## Chose the max because at lower speeds, the forced flow model is inaccurate
        surfaces["h_out"] = pd.concat([h_out_quiescent, h_out_flow], axis=1).max(axis=1)
        
        # Q of the composite wall
        surfaces["Q"] = HeatFlux(
            surfaces["T_in"], surfaces["T_out"],
            surfaces["R"] + 1 / surfaces["h_in"] + 1 / surfaces["h_out"],
            surfaces["Area"])

        # Update Surface Temp
        surfaces["T_si"] = SurfaceTemp(surfaces["T_in"], -surfaces["Q"],
                                       1 / surfaces["h_in"], surfaces["Area"])
        surfaces["T_so"] = SurfaceTemp(surfaces["T_out"], surfaces["Q"],
                                       1 / surfaces["h_out"], surfaces["Area"])

        # Convergence of (Q_conduction - Q_composite) in L2 norm
        Q = HeatFlux(surfaces["T_si"], surfaces["T_so"], surfaces["R"],
                     surfaces["Area"])
        
        error_new = np.linalg.norm(surfaces["Q"] - Q)

        ++iterations

    return surfaces["Q"].sum()

v_inverse_surface(surfaces["surface_in"])
Out[8]:
array(['horizontal_hot_top', 'vertical', 'horizontal_hot_bottom',
       'vertical', 'vertical', 'horizontal_hot_top', 'vertical'],
      dtype='<U21')

Solution and Parametric study

Having done all this work, we are now able to compute the heat transfer for various wind velocities and temperatures. So let's run it for temperatures between -15°C and 20°C and wind velocities between 0 and 60 km/h, assuming a comfortable inside temperature of 20°C.

Be carefull : the free convection equations for horizontal planes depend on the orientation of the hot and cold surfaces and already contain the gravity factor. This model is set for hot inside and cold outside only. To study what happens when the outside is hot and the inside is cold, one would have to change the model for the free convection of horizontal planes.

In [15]:
# Build the parametric tables
N = 40
T_in = deg2K(20)

T_out = deg2K(np.linspace(-15, 20, N))
v_out = np.linspace(0.1, 60e3/3600, 2*N)


Q = []

# Solve
for T in T_out:
    out_row = []
    for v in v_out:
        out_row.append(solver(T_in, T, v, surfaces))
        
    Q.append(out_row)
In [16]:
# Plot
trace = go.Heatmap(z=Q, 
                   x=v_out, 
                   y=T_out-273, 
                   colorscale=zero_centered_heatmap(Q),
                   colorbar = dict(
                    title = 'Surface Heat Transfer (W)',
                    titleside = 'top'
                    )
                  )

layout = go.Layout(
    title='Heat transfer through the boundaries by conduction and convection',
    xaxis = dict(title="Wind speed (m/s)"),
    yaxis = dict(title="Temperature (°C)" )
)

iplot(go.Figure(data=[trace], layout=layout), show_link=False, filename="convection")
Out[16]:

Results validation

The bad thing here is we don't have measurements data to validate the model accuracy. The first thing is to validate the consistency of the model behavior as parameters varies and the intermediate values.

We already know that most van heaters have a 3 kW design power, and that some vanlifers find them more than suitable. Looking above at the parametric study, we find a maximum loss of 1200 W @(-15°C; 16.67 m/s), so it seems reasonable.

We then want to ensure the strict monotony of the results.

In [13]:
# Build the parameters tables
N = 100
T_in = deg2K(20)

T_out = deg2K(np.linspace(-20, 20, N))
v_out = 1

data=[]

for v_out in [0.1, 0.5, 1, 2, 5, 10, 20]:
    
    Q = []
    
    # Solve
    for T in T_out:
        Q.append(solver(T_in, T, v_out, surfaces))

    # Plot
    trace = go.Scatter(
                       x=T_out-273, 
                       y=Q, 
                    name="%f m/s" % v_out
                      )
    data.append(trace)

layout = go.Layout(
    title='Heat transfer through the boundaries by conduction and convection at constant wind speed',
    xaxis = dict(title="Temperature (°C)"),
    yaxis = dict(title="Heat Losses (W)" )
)

    
iplot(go.Figure(data=data, layout=layout), show_link=False, filename="validation-speed")
Out[13]:
In [14]:
# Build the parameters tables
N = 100
T_in = deg2K(20)
v_out = np.linspace(0.1, 40, N)

data = []

for T_out in [-15, -10, -5, 0, 5, 10, 15]:
    
    Q = []

    # Solve
    for v in v_out:
        Q.append(solver(T_in, deg2K(T_out), v, surfaces))

    # Plot
    trace = go.Scatter(
                       x=v_out, 
                       y=Q, 
                        name="%i °C" % T_out
                      )
    
    data.append(trace)

layout = go.Layout(
    title='Heat transfer through the boundaries by conduction and convection at constant temperature',
    yaxis = dict(title="Heat Losses (W)" ),
    xaxis = dict(title="Wind speed (m/s)" )
)

iplot(go.Figure(data=data, layout=layout), show_link=False, filename="validation-temp")
Out[14]:

Overall, the behaviour is as expected but the curves show some bumps whereas they are supposed to be smooth. So we have some computation errors.

In [13]:
solver(deg2K(20), deg2K(-5), 0.1, surfaces)
surfaces
Out[13]:
Main dimension Secondary dimension Thickness blind flow_in surface_in L_in flow_out surface_out L_out Area R T_in T_si T_so T_out h_in h_out Q
roof 1.8 4.8 0.0007 True free horizontal_cold_bottom 0.654545 external parallel 4.8 8.64 1.253861 293.16 289.732 271.488 268.16 4.24404 4.371730 -125.711
sides 11.1 1.9 0.0007 True free vertical 1.9 external parallel 4.8 21.09 1.253861 293.16 289.082 272.131 268.16 3.31508 3.404733 -285.12
floor 1.8 5.5 0.005 True free horizontal_cold_top 0.654545 external parallel 4.8 9.9 1.253955 293.16 285.62 275.687 268.16 1.05055 1.052377 -78.4204
windchill 1 1.79 0.002 False free vertical 1 external perpendicular 1.79 1.79 0.257839 293.16 283.692 277.512 268.16 2.53152 2.562974 -42.9039
backdoors 1.84 1.56 0.0007 True free vertical 1.9 external perpendicular 1.56 2.8704 1.253861 293.16 289.082 272.131 268.16 3.31508 3.404733 -38.8055
vents 0.8 0.8 0.002 False free horizontal_cold_bottom 0.654545 external parallel 4.8 0.64 0.001429 293.16 280.514 280.437 268.16 4.24404 4.371730 -34.3498
windows 0.65 3.39 0.002 False free vertical 1.9 external parallel 4.8 2.2035 0.257839 293.16 284.321 276.766 268.16 3.31508 3.404733 -64.5658

There are some things we want to ensure on the results to validate the calculus :

  • $h_{in}$ is expected between 1 to 5 W/m²/K (OK)
  • $h_{out}$ is expected between 5 to 40 W/m²/K (OK)
  • there is a temperature gradient between the inside surfaces of the roof and the floor (OK)
  • the inside surfaces of transparent bodies are colder than the blind bodies (OK)

Taking account of the vents

The AFNOR standard makes a permanent ventilation mandatory in camping van, with a 100 cm² opening on the roof and a 15 cm² opening on the bottom of the side. This creates a cold air flow inside the van runned by the pressure difference between hot and cold air.

We can consider that the air at the outlet is at the same temperature as the roof and will go through the vent vertically, leading to a mass debit $\dot m$. Using the conservation of the mass, an equal incoming debit will appear at the vent on the bottom.

vent

A natural airflow appears as soon as a pressure gradient or a density gradient occures. We suppose here, the wind is parallel to the vents thus the flow is entirely buyoancy-driven.

We use the equations adapted from the MIT:

$$ \Delta P = P_{in} - P_{out} = g \frac{\rho_{out} - \rho_{in}}{\rho_{in}} \Delta H$$

where $\Delta H$ is the vertical distance between both vents (1.8 m).

Then, the volumetric rate is, considering no discharge :

$$ \dot V = A_{avg} \sqrt{2 \Delta P}$$

where $A_{avg}$ is the average area of the inlet and the outlet, as a rough approximation.

Then, the massic rate is :

$$ \dot m = \dot V \rho_{avg}$$

where $\rho_{avg}$ is the average density of air at the inlet and at the oulet, as a rough approximation.

Once $\dot m$ is found, we can solve the general energy equation between the ouside and the outside again trough the inside :

$$ \hat h_{in} + 1/2 U_{in}^2 = \hat h_{out} + 1/2 U_{out}^2 - q $$

where :

  • $ \hat h$ are the average massic enthalpies in J/kg/s, considering the air flow is composed of dry air and saturated vapour.
  • $q$ is the heat transferred to the air in J/kg/s

We will consider a relative humidity of $\phi = 80 \%$ in the air. The absolute humidity is :

$$ \omega = \frac{0.622 \phi P_g}{P - \phi P_g}$$

where :

  • $P$ is the atmospheric pressure
  • $P_g$ is the saturation pressure of water, 2.3392 kPa at 20°C

Then :

$$ \hat h = h_{air} + \omega h_{vapour}$$

where we use the correlation $h_{vapour}(T) = 2500.9 + 1.82 T$ kJ/kg.

Finally :

$$ -q = \hat h_{in} - \hat h_{out} - 1/2 U_{avg}^2$$

and the heat losses that we search are $Q = \dot m q$.

In [17]:
phi = 0.9
omega = 0.622 * phi * 2339.2 / (101500 - phi * 2339.2)

def enthalpy_vapour(T):
    return 2500.9 + 1.82 * T

omega
Out[17]:
0.013174584726432147
In [18]:
z_out = 1.8 # m
A_out = 0.01 # m²
A_in = 0.0015 # m²

def vent_losses(T_in, T_out):
    dp = 9.81 * (rho(T_out) - rho(T_in)) * z_out / rho(T_in)
    u = np.sqrt(2 * dp)
    m = (A_out * rho(T_in) + A_in * rho(T_out)) /2 * u
    
    h_in = enthalpy(T_in) + omega * enthalpy_vapour(T_in)
    h_out = enthalpy(T_out) + omega * enthalpy_vapour(T_out)
    
    return m * (h_out -  h_in + 0.5 * u**2)

# Estimate the losses for +20°C inside and -20°C outside
vent_losses(293, 253)
Out[18]:
-0.62688148631011342
In [19]:
# Build the parameters tables
N = 40
T_out = deg2K(np.linspace(-15, 20, N))
v_out = np.linspace(0.1, 60e3/3600, 2*N)
T_in = deg2K(20)

Q = []

for T in T_out:
    out_row = []
    for v in v_out:
        out_row.append(solver(T_in, T, v, surfaces) + vent_losses(T_in, T))
        
    Q.append(out_row) 

# Plot
trace = go.Heatmap(z=Q, 
                   x=v_out, 
                   y=T_out-273, 
                   colorscale=zero_centered_heatmap(Q),
                   colorbar = dict(
                    title = 'Surface Heat Transfer (W)',
                    titleside = 'top'
                    )
                  )

layout = go.Layout(
    title='Heat transfer through the boundaries by conduction and convection',
    xaxis = dict(title="Wind speed (m/s)"),
    yaxis = dict(title="Temperature (°C)" )
)
data=[trace]
iplot(go.Figure(data=[trace], layout=layout), show_link=False)
Out[19]:

We notice that the vents change almost nothing in the global heat losses. As an additional study, we wonder how efficient these vents are to renew the inside air depending on the temperature.

In [20]:
Q = []

T_out = deg2K(np.linspace(-20, 18, N))

for T in T_out:
    Q.append(vent_losses(T_in, T)*24)

        
# Plot
trace = go.Scatter(
                   x=T_out-273, 
                   y=Q, 
                  )

layout = go.Layout(
    title='Total energy lost per day due to the ventilation',
    yaxis = dict(title="Energy per day (Wh)" ),
    xaxis = dict(title="Outside temperature (°C)" )
)

iplot(go.Figure(data=[trace], layout=layout), show_link=False)
Out[20]:
In [21]:
V = 15 # m³

def vents_efficiency(T_in, T_out):
    dp = 9.81 * (rho(T_out) - rho(T_in)) * z_out / rho(T_in)
    u = np.sqrt(2 * dp)
    v = (A_out * rho(T_in) + A_in * rho(T_out)) /2 * u / (rho(T_out) + rho(T_in))
    return V/v/3600

time = []

T_out = deg2K(np.linspace(-20, 18, N))

for T in T_out:
    time.append(vents_efficiency(T_in, T))

        
# Plot
trace = go.Scatter(
                   x=T_out-273, 
                   y=time, 
                  )

layout = go.Layout(
    title='Efficiency of the natural ventilation depending on the outside temperature',
    yaxis = dict(title="Time needed to renew the all inside air volume (h)" ),
    xaxis = dict(title="Outside temperature (°C)" )
)

iplot(go.Figure(data=[trace], layout=layout), show_link=False)
Out[21]:

Studying energy production

On YouTube, just one vanlifer said he used an windturbine to produce his electricity (plus the now classic solar panels) and said it provided "more than enough" energy for his needs. We would like to study in this section the combined production and consumption of energy to determine what combination of renewable energies are required depending on the weather.

Windmills models

We will model two different windturbines from the manufacturer datasheet (one 400 W, one 600 W), using their wind speed vs. generated power characteristic and consider 5 m² of solar panels having 20 % efficiency. We consider the general efficiency of the batteries and electrical system is 85 %.

In [22]:
speeds = np.array([2, 3.5, 4.5, 6, 7, 8.5, 9.5, 10, 11])

DEBUG = False
GRAPH = False

def windmill_600(v):
    """Power generation in W depending of the wind speed in m/s - 600 W windmill"""

    coef = [  3.489621259632e-02,  -2.211788459198e+00,   3.213943087260e+01,
  -7.872620986471e+01,   1.003488153562e+02]
    
    if v >= 11:
        result = 690
    elif v <= 2:
        result = 0
    else:
        result = np.polyval(coef, v)
    
    if DEBUG:
        # Do the actual fitting
        values = np.array([250, 500, 1000, 1500, 2000, 2500, 2750, 3000, 3000]) * 1e3 / (365*12)
        
        print(speeds.shape[0], values.shape[0])
        debug_result, coef_updated = polyfit(values, speeds, order=4)
        
        data = [
            go.Scatter(x=speeds, y=values, name="Tabulated"),
            go.Scatter(x=speeds, y=debug_result, name="Interpolated")
        ]

        iplot(data, show_link=False, filename='test')

        if not np.allclose(coef, coef_updated):
            raise RuntimeError(
                "The hard-coded fitting parameters do not match the updated fitting"
            )

        print("power = ", result)

    return result

windmill_600(1.5)
Out[22]:
0

img

In [23]:
#http://www.naviclub.com/electricite/eolienne/sun-watts/EO-420-Hybride/

speeds = np.array([2, 5.8, 8, 9, 11, 14, 16, 17, 18])

DEBUG = False

def windmill_400(v):
    """Power generation in W depending of the wind speed in m/s - 400 W windmill"""

    coef = [  5.317592058356e-03,  -2.417665862664e-01,   3.549419511456e+00,
  -1.761809621471e+01,   4.971569263348e+01,  -5.341557066638e+01]
    
    if v >= 18:
        result = 500
    elif v <= 2:
        result = 0
    else:
        result = np.polyval(coef, v)
    
    if DEBUG:
        # Do the actual fitting
        values = np.array([0, 100, 200, 300, 400, 500, 500, 500, 500])
        
        print(speeds.shape[0], values.shape[0])
        debug_result, coef_updated = polyfit(values, speeds, order=5)
        
        data = [
            go.Scatter(x=speeds, y=values, name="Tabulated"),
            go.Scatter(x=speeds, y=debug_result, name="Interpolated")
        ]

        iplot(data, show_link=False, filename='test')

        if not np.allclose(coef, coef_updated):
            raise RuntimeError(
                "The hard-coded fitting parameters do not match the updated fitting"
            )

        print("power = ", result)

    return result

windmill_400(5)
Out[23]:
63.901284831132507

img

Climate probability

Using several cities statistics during January (average wind speed, average solar radiation power, average temperature), we will determine the most probable using conditions during winter and estimate the difference between the expected energy consumption and expected energy production to see if a such van could be self-sustainable durinng winter and at what cost.

This probabilities assume temperature and wind speed are independant and follow a normal law. We already know that's not accurate, the problem is we don't have access to free datasets on climate in France and Europe to perform statistics on them before computing probabilities.

In [28]:
# https://fr.weatherspark.com/

winter = {
    "nancy":{
        "T": deg2K(1.2),
        "v": 17.5e3/3600,
        "eps": 0.8, #kWh/day/m²
        "lat": 48.691136187962854,
        "long": 6.187180950316247
    },
    "sarrebruck" : {
        "T": deg2K(0.9),
        "v": 17e3/3600,
        "eps": 0.7, #kWh/day/m²
        "lat": 49.23716799059577,
        "long": 6.997716098316573
    },
    "strasbourg": {
        "T": deg2K(0),
        "v": 14e3/3600,
        "eps": 0.9, #kWh/day/m²
        "lat": 48.58194484068934,
        "long": 7.751026153564453
    },
    "berlin": {
        "T": deg2K(-1),
        "v": 20e3/3600,
        "eps": 0.6, #kWh/day/m²
        "lat": 52.51992557040081,
        "long": 13.413937169115968
    },
    "lyon": {
        "T": deg2K(4),
        "v": 14e3/3600,
        "eps": 1.2, #kWh/day/m²
        "lat": 45.762829178117336,
        "long": 4.856873533886812
    },
    "brest": {
        "T": deg2K(7),
        "v": 28e3/3600,
        "eps": 0.9, #kWh/day/m²
        "lat": 48.390624357030404,
        "long": -4.4819551910398445
    },
    "munchen": {
        "T": deg2K(-1),
        "v": 13e3/3600,
        "eps": 1,
        "lat": 48.13642753535,
        "long": 11.569792690803752
    },
    "paris": {
        "T": deg2K(2.5),
        "v": 18.8e3/3600,
        "eps": 0.8,
        "lat": 48.856233862689706,
        "long": 2.363313934141388
    },
    "lille": {
        "lat": 50.642460107903126,
        "long": 3.0602441970460745,
        "T": deg2K(3),
        "v": 21e3/3600,
        "eps": 0.7
    },
    "zurich": {
        "T": deg2K(1),
        "v": 11e3/3600,
        "eps": 1.1,
        "lat": 47.36152183929631,
        "long": 8.565799653723275,
    },
    "bruxelles": {
        "lat": 50.85098962441577,
        "long": 4.353178480013639,
        "T": deg2K(2.5),
        "v": 21e3/3600,
        "eps": 0.6,
    },
    "reims": {
        "T": deg2K(3.4),
        "v": 18.2e3/3600,
        "eps": 0.7,
        "lat": 49.257632205302635,
        "long": 4.0286811837193
    },
    "rotterdam": {
        "lat": 51.92844318259155,
        "long": 4.469347284592459,
        "T": deg2K(2.5),
        "v": 24.1e3/3600,
        "eps": 0.6
    },
    "amsterdam": {
        "lat": 52.36915443865181,
        "long": 4.884081171311209,
        "T": deg2K(2.5),
        "v": 24.9e3/3600,
        "eps": 0.6
    },
    "koln": {
        "lat": 50.93849840403057,
        "long": 6.95983817280603,
        "T": deg2K(2.5),
        "v": 20e3/3600,
        "eps": 0.7
    },
    "hambourg": {
        "lat": 53.53431786340632,
        "long": 9.995789491488722,
        "T": deg2K(1.5),
        "v": 21e3/3600,
        "eps": 0.5
    },
    "luxembourg": {
        "lat": 49.596325051090766,
        "long": 6.1252895287357205,
        "T": deg2K(1.5),
        "v": 18.3e3/3600,
        "eps": 0.7
    },
    "london": {
        "lat": 51.50546328240563,
        "long": -0.12100990291219205,
        "T": deg2K(6),
        "v": 21.4e3/3600,
        "eps": 0.6
    },
    "franckfurt": {
        "lat": 50.11808109168232,
        "long": 8.739872702608466,
        "T": deg2K(1.5),
        "v": 16.5e3/3600,
        "eps": 0.7
    },
    "troyes": {
        "lat": 48.288612602856055,
        "long": 4.088183101551181,
        "T": deg2K(3.5),
        "v": 17.5e3/3600,
        "eps": 0.8
    },
    "dijon": {
        "lat": 47.322004442358306,
        "long": 5.065966304676181,
        "T": deg2K(2),
        "v": 16.2e3/3600,
        "eps": 1
    },
    "gottingen": {
        "lat": 51.58844654881783,
        "long": 9.920693406967644,
        "T": deg2K(1.5),
        "v": 18.7e3/3600,
        "eps": 0.6
    },
    "leipzig": {
        "lat": 51.36263409824682,
        "long": 12.392617235092644,
        "T": deg2K(1),
        "v": 19.4e3/3600,
        "eps": 0.7
    },
    "hanover": {
        "lat": 52.36986923691711,
        "long": 9.808083543686394,
        "T": deg2K(1.5),
        "v": 20.8e3/3600,
        "eps": 0.6
    },
    "munster": {
        "lat": 51.97857460538177,
        "long": 7.586453404233225,
        "T": deg2K(2.5),
        "v": 20.8e3/3600,
        "eps": 0.6
    },
    "milano": {
        "lat": 45.47485948464145,
        "long": 9.190029332820814,
        "T": deg2K(2),
        "v": 7.3e3/3600,
        "eps": 1.4
    },
    "clermont": {
        "lat": 45.780920018053436,
        "long": 3.127467577843845,
        "T": deg2K(3.5),
        "v": 14.7e3/3600,
        "eps": 1.2
    },
    "groningue": {
        "lat": 53.235604792206985,
        "long": 6.5671440363573765,
        "T": deg2K(2.5),
        "v": 26.2e3/3600,
        "eps": 0.5
    },
    "rennes": {
        "lat": 48.11562988832791,
        "long": -1.6761851834323807,
        "T": deg2K(5.5),
        "v": 19.2e3/3600,
        "eps": 0.9
    },
    "nantes": {
        "lat": 47.23113657192925,
        "long": -1.5469866037800557,
        "T": deg2K(6),
        "v": 19e3/3600,
        "eps": 1
    },
    "vannes": {
        "lat": 47.66576117807731,
        "long": -2.7609758615925557,
        "T": deg2K(6),
        "v": 22e3/3600,
        "eps": 1
    },
    "bourges": {
        "lat": 47.07479801076057,
        "long": 2.411113662197522,
        "T": deg2K(4.5),
        "v": 17e3/3600,
        "eps": 1        
    },
    "bordeaux": {
        "lat": 44.8685225722807,
        "long": -0.599140244052478,
        "T": deg2K(7),
        "v": 12e3/3600,
        "eps": 1.3
    },
    "toulouse": {
        "lat": 43.625140844192494,
        "long": 1.455303115322522,
        "T": deg2K(5.5),
        "v": 14e3/3600,
        "eps": 1.4
    },
    "montpellier": {
        "lat": 43.6907848718809,
        "long": 3.8800485808694702,
        "T": deg2K(7),
        "v": 16e3/3600,
        "eps": 1.7
    },
    "pau": {
        "lat": 43.32424902668016,
        "long": -0.38264673163052976,
        "T": deg2K(5.5),
        "v": 11e3/3600,
        "eps": 1.5
    },
    "saragosse": {
        "lat": 41.704775479596734,
        "long": -0.8550588410055298,
        "T": deg2K(7),
        "v": 17e3/3600,
        "eps": 2
    },
    "torino": {
        "lat": 45.130319345557346,
        "long": 7.689202294523625,
        "T": deg2K(3),
        "v": 6.3e3/3600,
        "eps": 1.6
    },
    "genes": {
        "lat": 44.45985400803417,
        "long": 8.952630028898625,
        "T": deg2K(8),
        "v": 12e3/3600,
        "eps": 1.5
    },
    "bologna": {
        "lat": 44.55387494722501,
        "long": 11.336663232023625,
        "T": deg2K(3),
        "v": 11e3/3600,
        "eps": 1.5
    }
    
}
In [29]:
def gaussian(x, mu, sigma):
    """Probability of random variable x considering a normal process of average mu and std sigma"""
    
    return np.exp(-(x - mu)**2 / (2 * sigma**2))

def probability(T, Tavg, v, vavg, sigma_T, sigma_v):
    """Construct a probability matrix for 2 independent variables"""
    
    proba_T = gaussian(T, Tavg, sigma_T)
    proba_v = gaussian(v, vavg, sigma_v)
    result = np.outer(proba_T, proba_v)
    
    return result / result.sum()


city = "london"
proba = probability(T_out, winter[city]["T"], v_out, winter[city]["v"], 3, 2)
trace = go.Heatmap(z=proba, 
                   x=v_out, 
                   y=T_out-273, 
                   colorscale=zero_centered_heatmap(proba),
                   colorbar = dict(
                    title = 'Probability of happening',
                    titleside = 'top'
                    )
                  )

layout = go.Layout(
    title='Probability of the using conditions in %s, January' % city,
    xaxis = dict(title="Wind speed (m/s)"),
    yaxis = dict(title="Temperature (°C)" )
)
data=[trace]
iplot(go.Figure(data=[trace], layout=layout), show_link=False)
Out[29]:
In [77]:
# Build the parameters tables
N = 40
T_out = deg2K(np.linspace(-15, 15, N))
v_out = np.linspace(0.1, 60e3/3600, 2*N)
T_in = deg2K(np.array([15, 17, 19]))

Q = []       
P_600 = []
P_400 = []
             
for T_i in T_in:
    Q_i = []
    for T_o in T_out:
        out_row = []
        out_row_400 = []
        out_row_600 = []

        for v in v_out:
            E_600 = windmill_600(v) * 0.85
            E_400 = windmill_400(v) * 0.85

            out_row_400.append(E_400)
            out_row_600.append(E_600)
            out_row.append(solver(T_i, T_o, v, surfaces) + vent_losses(T_i, T_o))


        Q_i.append(out_row) 
        P_600.append(out_row_600) 
        P_400.append(out_row_400) 
             
    Q.append(Q_i)
    
Q = np.array(Q)
P_400 = np.array(P_400)
P_600 = np.array(P_600)

P_400 = P_400[0:N, ...]
P_600 = P_600[0:N, ...]
In [78]:
# Plot 400
trace = go.Heatmap(z=P_400 + Q[2, ...], 
                   x=v_out, 
                   y=T_out-273, 
                   colorscale=zero_centered_heatmap(P_400 + Q[1, ...]),
                   colorbar = dict(
                    title = 'Difference of power(W)',
                    titleside = 'top'
                    )
                  )

layout = go.Layout(
    title='Difference between power production and heat losses - 400 W',
    xaxis = dict(title="Wind speed (m/s)"),
    yaxis = dict(title="Temperature (°C)" )
)
data=[trace]
iplot(go.Figure(data=[trace], layout=layout), show_link=False)

# Plot 600
trace = go.Heatmap(z=P_600 + Q[2, ...], 
                   x=v_out, 
                   y=T_out-273, 
                   colorscale=zero_centered_heatmap(P_600 + Q[1, ...]),
                   colorbar = dict(
                    title = 'Difference of power(W)',
                    titleside = 'top'
                    )
                  )

layout = go.Layout(
    title='Difference between power production and heat losses - 600 W',
    xaxis = dict(title="Wind speed (m/s)"),
    yaxis = dict(title="Temperature (°C)" )
)
data=[trace]
iplot(go.Figure(data=[trace], layout=layout), show_link=False, filename="consumption-production")
Out[78]:
In [79]:
print('{:<15s}{:>12s}{:>12s}'.format("CITY", "PROD 400 W", "PROD 600 W"))
for key in winter.keys():
    proba = probability(T_out, winter[key]["T"], v_out, winter[key]["v"], 3, 2)
    
    E_400 = np.average(P_400, weights=proba) * 24 # Wh
    E_600 = np.average(P_600, weights=proba) * 24 # Wh
    
    Q_avg = Q[0, ...] * 8 + Q[1, ...] * 8 + Q[2, ...] * 8
    E_loss = np.average(Q_avg, weights=proba) # Wh
    E_people = 200 * 24
    
 
    print('{:<15s}{:>12d}{:>12d}{:>8s}'.format(key,
                                               int(E_400 + E_loss + E_people),
                                               int(E_600 + E_loss + E_people),
                                               "Wh/day"))
CITY             PROD 400 W  PROD 600 W
nancy                 -5182       -1508  Wh/day
sarrebruck            -5483       -1935  Wh/day
strasbourg            -6449       -3649  Wh/day
berlin                -6615       -2323  Wh/day
lyon                  -3317        -518  Wh/day
brest                  1813        7564  Wh/day
munchen               -7283       -4722  Wh/day
paris                 -3955          44  Wh/day
lille                 -3164        1362  Wh/day
zurich                -5761       -3650  Wh/day
bruxelles             -3569         957  Wh/day
reims                 -3326         523  Wh/day
rotterdam             -2886        2289  Wh/day
amsterdam             -2687        2631  Wh/day
koln                  -3756         536  Wh/day
hambourg              -4384         142  Wh/day
luxembourg            -4834        -959  Wh/day
london                 -683        3934  Wh/day
franckfurt            -5057       -1635  Wh/day
troyes                -3345         328  Wh/day
dijon                 -4690       -1344  Wh/day
gottingen             -4776        -801  Wh/day
leipzig               -5073        -926  Wh/day
hanover               -4422          58  Wh/day
munster               -3608         872  Wh/day
milano                -4985       -3575  Wh/day
clermont              -3650        -679  Wh/day
groningue             -2346        3178  Wh/day
rennes                -1504        2594  Wh/day
nantes                -1147        2901  Wh/day
vannes                 -552        4198  Wh/day
bourges               -2624         924  Wh/day
bordeaux              -1188        1142  Wh/day
toulouse              -2167         631  Wh/day
montpellier            -827        2468  Wh/day
pau                   -2343        -232  Wh/day
saragosse              -694        2853  Wh/day
torino                -4225       -2972  Wh/day
genes                  -470        1860  Wh/day
bologna               -4227       -2116  Wh/day

This table shows the prospective difference between energy consumption and energy generation in various cities in Western Europe depending the average temperature, wind and solar power during January and 2 people assumed te be 100 W heaters each. The 600 W windturbine shows that it may be self-sustainable in average in some cities, but the 400 W one is not enough.

In [80]:
for key in winter.keys():
    proba = probability(T_out, winter[key]["T"], v_out, winter[key]["v"], 3, 2)
    
    E_solar = winter[key]["eps"] * .15 * 6 * 1000 * 0.85
    E_people = 200 * 16
    E_600 = np.average(P_600, weights=proba) * 24 # Wh
    Q_avg = Q[0, ...] * 8 + Q[1, ...] * 8 + Q[2, ...] * 8
    E_loss = np.average(Q_avg, weights=proba) # Wh
    
    winter[key]["energy"] = E_600 + E_solar + E_loss + E_people

cities = pd.DataFrame.from_dict(winter, orient="index")
cities
Out[80]:
T v eps lat long energy
amsterdam 275.66 6.916667 0.6 52.369154 4.884081 1490.554194
berlin 272.16 5.555556 0.6 52.519926 13.413937 -3464.307777
bologna 276.16 3.055556 1.5 44.553875 11.336663 -2568.901429
bordeaux 280.16 3.333333 1.3 44.868523 -0.599140 536.687334
bourges 277.66 4.722222 1.0 47.074798 2.411114 89.040630
brest 280.16 7.777778 0.9 48.390624 -4.481955 6653.469701
bruxelles 275.66 5.833333 0.6 50.850990 4.353178 -183.896297
clermont 276.66 4.083333 1.2 45.780920 3.127468 -1361.511554
dijon 275.16 4.500000 1.0 47.322004 5.065966 -2179.859117
franckfurt 274.66 4.583333 0.7 50.118081 8.739873 -2700.363591
genes 281.16 3.333333 1.5 44.459854 8.952630 1407.975653
gottingen 274.66 5.194444 0.6 51.588447 9.920693 -1942.834074
groningue 275.66 7.277778 0.5 53.235605 6.567144 1960.530596
hambourg 274.66 5.833333 0.5 53.534318 9.995789 -1074.834425
hanover 274.66 5.777778 0.6 52.369869 9.808084 -1082.486243
koln 275.66 5.555556 0.7 50.938498 6.959838 -528.492285
leipzig 274.16 5.388889 0.7 51.362634 12.392617 -1991.258555
lille 276.16 5.833333 0.7 50.642460 3.060244 297.825261
london 279.16 5.944444 0.6 51.505463 -0.121010 2793.678106
luxembourg 274.66 5.083333 0.7 49.596325 6.125290 -2023.971953
lyon 277.16 3.888889 1.2 45.762829 4.856874 -1200.239019
milano 275.16 2.027778 1.4 45.474859 9.190029 -4104.094454
montpellier 280.16 4.444444 1.7 43.690785 3.880049 2168.927898
munchen 272.16 3.611111 1.0 48.136428 11.569793 -5557.890226
munster 275.66 5.777778 0.6 51.978575 7.586453 -268.749552
nancy 274.36 4.861111 0.8 48.691136 6.187181 -2496.303007
nantes 279.16 5.277778 1.0 47.231137 -1.546987 2066.153971
paris 275.66 5.222222 0.8 48.856234 2.363314 -943.996588
pau 278.66 3.055556 1.5 43.324249 -0.382647 -684.814011
reims 276.56 5.055556 0.7 49.257632 4.028681 -540.849771
rennes 278.66 5.333333 0.9 48.115630 -1.676185 1682.592872
rotterdam 275.66 6.694444 0.6 51.928443 4.469347 1148.549708
saragosse 280.16 4.722222 2.0 41.704775 -0.855059 2783.504877
sarrebruck 274.06 4.722222 0.7 49.237168 6.997716 -2999.683337
strasbourg 273.16 3.888889 0.9 48.581945 7.751026 -4561.466030
torino 276.16 1.750000 1.6 45.130319 7.689202 -3348.018649
toulouse 278.66 3.888889 1.4 43.625141 1.455303 102.833903
troyes 276.66 4.861111 0.8 48.288613 4.088183 -659.632022
vannes 279.16 6.111111 1.0 47.665761 -2.760976 3363.642864
zurich 274.16 3.055556 1.1 47.361522 8.565800 -4408.523404
In [81]:
data = [ dict(
    lat = cities["lat"],
    lon = cities["long"],
    text = cities.energy.astype(str) + " Wh/day",
    hoverinfo = "text",
    marker = dict(
        color = cities.energy,
        opacity = 1,
        size = 20,
        colorscale = zero_centered_heatmap(cities.energy),
        colorbar = dict(
            ticksuffix = " Wh/day",
        ),
    ),
    type = 'scattergeo',
) ]

layout = dict(
    geo = dict(
        scope = 'europe',
        projection = dict(
            type = 'mercator',
            rotation = dict(
                lon = 0
            )
        ),
        landcolor = "rgb(250, 250, 250)",
        #subunitcolor = "rgb(255, 255, 255)",
        #countrycolor = "rgb(255, 255, 255)",
        #lakecolor = "rgb(255, 255, 255)",
        showsubunits = True,
        showcountries = True,
        showland=True,
        showocean=True,
        showcoastlines=True,
        showlakes=True,
        lonaxis = dict(
            range= [ -5, 15],
        ),
        lataxis = dict (
            range= [ 40.0, 55 ],
        )
    ),
    title = 'Difference between expected heating energy consumption and expected renewable energy generation',
)
fig = { 'data':data, 'layout':layout }
iplot(fig, filename='precipitation')
Out[81]:

Conclusion

This work was primarily focused on migrating 1D heat transfer evaluations from MS Excel spreadsheet to Python in order to perform a parametric study.

We recall that the study was conducted assuming blinded windows (for the night), open vents and no direct sunlight (no heating by radiation and greenhouse effect), assuming 2 persons inside as 100 W heaters. This is then a worse case scenario. The final energy estimation is conducted assuming 8 hours of heating to 21°C (morning and evening), 8 hours of heating at 18°C (night), and 8 hours of heating at 15°C (business hours assuming people are not there).

We notice that the main responsible for the heat losses are the sides (40 %), the roof (17 %), the floor (13 %) and the windows (10 %). They are easy to over-insulate, so all the effort should be put there. The insulation of the floor is rather a matter of comfort, so in a conflict to keep a maximal headroom, the roof should have precedence over the floor. Adding just 1 cm of insulation improves the losses through the roof and the sides by 18 to 20 %, and adding 2 cm of insulation on the blinded windows and windchill (for the night) improves the losses through the glass by 55 %. Overall, it's a 21 % improvement on total heat losses.

It is noticeable that the wind adds no more than 25 % to the heat losses but changes severely the potential of energy generation. Being mobile in a van should encourage to do energetic tourism, looking for more wind to be self-sustainable during winter. In a such scenario, wind should be prefered to warmer temperatures, but there is a trade-off. Looking at the map, we see that the littoral is preferable and regions close to the Alps should be avoided as these mountains are a wind barrier.

The solar panels won't be of much use during winter, except close to the Mediterranean Sea. It is then surprising to see many campers automatically add them on the roof while a windmill would make such a bigger difference : colder Northern Europe shows better energy generation possibilities (Holland, Belgium, England) than the Northern of Italy (except the mountains area).

As bad as the idea of using permanents vents sounded intuitively, we see that the heat loss it produces is less than 1 W even at very low temperatures. One should not miss this opportunity to get healtier air for energy savings motives, since it is worth 0.06 % of the total heat losses. This is of course different if the wind blows directly through the inlet vent. Assuming no friction on the vents (which is wrong because they have grids), the inside air volume of the van should be totally renewed between 30 min (at -15 °C) to 2 h (at 15°C). In colder weathers, it would be a good idea to reduce the area of the vents to improve the comfort while keeping a good ventilation.

This work shows that it is possible to make a self-sustainable van even during winter, depending the region of use. The price is not cheap though, the windturbine itself costs 2500 CAD(1600 €) and the solar panels would be 500 CAD/m², so 3000 CAD in total (2000 €). This is without the battery, MPPT, BMS, etc. Considering the kWh is at 0.15 € in Europe, it would be profitable after 24 000 kWh, and considering an average of 6.25 kWh produced/day (in January), that would be 10.5 years. If the solar panels last that long…

In [82]:
# Average losses during January in Europe
T_mean = cities["T"].mean()
v_mean = cities["v"].mean()
solver(deg2K(20), T_mean, v_mean, surfaces)
surfaces
Out[82]:
Main dimension Secondary dimension Thickness blind flow_in surface_in L_in flow_out surface_out L_out Area R T_in T_si T_so T_out h_in h_out Q
roof 1.8 4.8 0.0007 True free horizontal_cold_bottom 0.654545 external parallel 4.8 8.64 1.253861 293.16 290.339 277.265 276.3975 3.69696 12.021236 -90.0929
sides 11.1 1.9 0.0007 True free vertical 1.9 external parallel 4.8 21.09 1.253861 293.16 289.733 277.227 276.3975 2.91009 12.021236 -210.344
floor 1.8 5.5 0.005 True free horizontal_cold_top 0.654545 external parallel 4.8 9.9 1.253955 293.16 285.861 276.986 276.3975 0.969551 12.021236 -70.0636
windchill 1 1.79 0.002 False free vertical 1 external perpendicular 1.79 1.79 0.257839 293.16 283.797 278.274 276.3975 2.28787 11.417640 -38.3444
backdoors 1.84 1.56 0.0007 True free vertical 1.9 external perpendicular 1.56 2.8704 1.253861 293.16 289.737 277.247 276.3975 2.91009 11.729014 -28.5931
vents 0.8 0.8 0.002 False free horizontal_cold_bottom 0.654545 external parallel 4.8 0.64 0.001429 293.16 280.392 280.324 276.3975 3.69696 12.021236 -30.2106
windows 0.65 3.39 0.002 False free vertical 1.9 external parallel 4.8 2.2035 0.257839 293.16 284.747 278.434 276.3975 2.91009 12.021236 -53.9484
In [ ]:
# Average daily energy production

T_mean = cities["T"].mean()
v_mean = cities["v"].mean()
eps_mean = cities["eps"].mean()
proba = probability(T_out, T_mean, v_out, v_mean, 3, 2)

E_solar = eps_mean * .20 * 6 * 1000 * 0.85
E_600 = np.average(P_600, weights=proba) * 24 # Wh

energy = E_600 + E_solar
energy

Disclaimer - Licence

This work is shared hoping it can help, in the spirit of opensource and opendata movements, so people could review it, criticize and correct it. Should you use the code and computations here, you have to report any mistake or inconsistency you may find. You use the results and calculations displayed here at your own risks and under your responsibility.

In [ ]: