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.
© 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").
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)
# 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()
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.
# 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
# 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
surfaces
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:
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 :
- assume $T_{s, i} = T_{s, o} = (T_{in} + T_{out})/2$
- compute $Gr_L$ then $Ra_L$ then $\overline{Nu_L}$ on both sides of the wall
- compute $\bar h$ on both sides of the wall
- compute $R_i$, $R_o$ and finally $Q$
- 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)$
- measure convergence with the error between $Q$ for the composite wall and $Q$ for the conduction part in L2 norm
- goto 2. and iterate until convergence
- 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
# 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.
@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"])
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.
# 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)
# 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")
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.
# 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")
# 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")
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.
solver(deg2K(20), deg2K(-5), 0.1, surfaces)
surfaces
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.
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$.
phi = 0.9
omega = 0.622 * phi * 2339.2 / (101500 - phi * 2339.2)
def enthalpy_vapour(T):
return 2500.9 + 1.82 * T
omega
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)
# 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)
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.
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)
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)
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 %.
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)
#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)
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.
# 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
}
}
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)
# 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, ...]
# 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")
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"))
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.
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
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')
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…
# 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
# 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.
