Analyse the heat losses of a Sprinter van

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 dimensionSecondary dimensionThicknessblindflow_insurface_inL_inflow_outsurface_outL_outAreaR
roof1.84.80.0007Truefreehorizontal_cold_bottom0.654545externalparallel4.88.641.253861
sides11.11.90.0007Truefreevertical1.9externalparallel4.821.091.253861
floor1.85.50.005Truefreehorizontal_cold_top0.654545externalparallel4.89.91.253955
windchill11.790.002Falsefreevertical1externalperpendicular1.791.790.257839
backdoors1.841.560.0007Truefreevertical1.9externalperpendicular1.562.87041.253861
vents0.80.80.002Falsefreehorizontal_cold_bottom0.654545externalparallel4.80.640.001429
windows0.653.390.002Falsefreevertical1.9externalparallel4.82.20350.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 dimensionSecondary dimensionThicknessblindflow_insurface_inL_inflow_outsurface_outL_outAreaRT_inT_siT_soT_outh_inh_outQ
roof1.84.80.0007Truefreehorizontal_cold_bottom0.654545externalparallel4.88.641.253861293.16289.732271.488268.164.244044.371730-125.711
sides11.11.90.0007Truefreevertical1.9externalparallel4.821.091.253861293.16289.082272.131268.163.315083.404733-285.12
floor1.85.50.005Truefreehorizontal_cold_top0.654545externalparallel4.89.91.253955293.16285.62275.687268.161.050551.052377-78.4204
windchill11.790.002Falsefreevertical1externalperpendicular1.791.790.257839293.16283.692277.512268.162.531522.562974-42.9039
backdoors1.841.560.0007Truefreevertical1.9externalperpendicular1.562.87041.253861293.16289.082272.131268.163.315083.404733-38.8055
vents0.80.80.002Falsefreehorizontal_cold_bottom0.654545externalparallel4.80.640.001429293.16280.514280.437268.164.244044.371730-34.3498
windows0.653.390.002Falsefreevertical1.9externalparallel4.82.20350.257839293.16284.321276.766268.163.315083.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]:
Tvepslatlongenergy
amsterdam275.666.9166670.652.3691544.8840811490.554194
berlin272.165.5555560.652.51992613.413937-3464.307777
bologna276.163.0555561.544.55387511.336663-2568.901429
bordeaux280.163.3333331.344.868523-0.599140536.687334
bourges277.664.7222221.047.0747982.41111489.040630
brest280.167.7777780.948.390624-4.4819556653.469701
bruxelles275.665.8333330.650.8509904.353178-183.896297
clermont276.664.0833331.245.7809203.127468-1361.511554
dijon275.164.5000001.047.3220045.065966-2179.859117
franckfurt274.664.5833330.750.1180818.739873-2700.363591
genes281.163.3333331.544.4598548.9526301407.975653
gottingen274.665.1944440.651.5884479.920693-1942.834074
groningue275.667.2777780.553.2356056.5671441960.530596
hambourg274.665.8333330.553.5343189.995789-1074.834425
hanover274.665.7777780.652.3698699.808084-1082.486243
koln275.665.5555560.750.9384986.959838-528.492285
leipzig274.165.3888890.751.36263412.392617-1991.258555
lille276.165.8333330.750.6424603.060244297.825261
london279.165.9444440.651.505463-0.1210102793.678106
luxembourg274.665.0833330.749.5963256.125290-2023.971953
lyon277.163.8888891.245.7628294.856874-1200.239019
milano275.162.0277781.445.4748599.190029-4104.094454
montpellier280.164.4444441.743.6907853.8800492168.927898
munchen272.163.6111111.048.13642811.569793-5557.890226
munster275.665.7777780.651.9785757.586453-268.749552
nancy274.364.8611110.848.6911366.187181-2496.303007
nantes279.165.2777781.047.231137-1.5469872066.153971
paris275.665.2222220.848.8562342.363314-943.996588
pau278.663.0555561.543.324249-0.382647-684.814011
reims276.565.0555560.749.2576324.028681-540.849771
rennes278.665.3333330.948.115630-1.6761851682.592872
rotterdam275.666.6944440.651.9284434.4693471148.549708
saragosse280.164.7222222.041.704775-0.8550592783.504877
sarrebruck274.064.7222220.749.2371686.997716-2999.683337
strasbourg273.163.8888890.948.5819457.751026-4561.466030
torino276.161.7500001.645.1303197.689202-3348.018649
toulouse278.663.8888891.443.6251411.455303102.833903
troyes276.664.8611110.848.2886134.088183-659.632022
vannes279.166.1111111.047.665761-2.7609763363.642864
zurich274.163.0555561.147.3615228.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 dimensionSecondary dimensionThicknessblindflow_insurface_inL_inflow_outsurface_outL_outAreaRT_inT_siT_soT_outh_inh_outQ
roof1.84.80.0007Truefreehorizontal_cold_bottom0.654545externalparallel4.88.641.253861293.16290.339277.265276.39753.6969612.021236-90.0929
sides11.11.90.0007Truefreevertical1.9externalparallel4.821.091.253861293.16289.733277.227276.39752.9100912.021236-210.344
floor1.85.50.005Truefreehorizontal_cold_top0.654545externalparallel4.89.91.253955293.16285.861276.986276.39750.96955112.021236-70.0636
windchill11.790.002Falsefreevertical1externalperpendicular1.791.790.257839293.16283.797278.274276.39752.2878711.417640-38.3444
backdoors1.841.560.0007Truefreevertical1.9externalperpendicular1.562.87041.253861293.16289.737277.247276.39752.9100911.729014-28.5931
vents0.80.80.002Falsefreehorizontal_cold_bottom0.654545externalparallel4.80.640.001429293.16280.392280.324276.39753.6969612.021236-30.2106
windows0.653.390.002Falsefreevertical1.9externalparallel4.82.20350.257839293.16284.747278.434276.39752.9100912.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 [ ]:
 
2018-05-29T21:25:19+00:005 May 2018|Categories: Non classé|Tags: , , |Comments Off on Analyse the heat losses of a Sprinter van

About the Author:

Développeur d'outils logiciels de traitement d'image pour darktable. Spécialiste calcul et modélisation thermodynamique chez Cellier Domesticus. Photographe. Pianiste. Développeur spécialisé en Python pour le calcul et la modélisation. Auteur de bouquins et de blog sur les sciences et la technologie. Expériences précédentes dans la fonction publique territoriale, les moteurs électriques industriels, les voitures solaires en fibre de carbone et le non-sens académique (maths sup, DUT).
ut Phasellus libero. venenatis, fringilla adipiscing facilisis elit. ipsum nec risus.