search

loonatick-src / SolarHeatingSim

1 stars 0 forks source link

readme

Rudimentary Solar Water Heating System Simulation

Running this repo

The only requirement is julia, preferably 1.9.

Assumptions and Simplifications

We make the following assumptions and simplifications

Model Parameters

@kwdef struct SWHSModel{T}
    Ns::Int= 10       # stratifications count for tank
    Nc::Int= 100      # number of points in the discretized mesh of collector    
    ρp::T  = 8.0e3    # density of plate
    ρf::T  = 1.0e3    # working fluid density
    ρe::T  = 1.5e3    # effective density (tank + wokring fluid)
    δ::T   = 0.1      # plate thickness
    W::T   = 1.0      # plate width
    L::T   = 2.0      # plate length
    H::T   = 10.0     # height of tank
    ΔH::T  = H/Ns     # height of each stratified layer
    Δy::T  = L/Nc
    dt::T  = 5.0      # tank diameter         
    dr::T  = 0.1      # riser pipe diameter
    dd::T  = dr       # down comer pipe diameter    
    kp::T  = 50.0     # thermal conductivity of plate
    Cpp::T = 450.0    # specific heat capacity of plate
    Cpf::T = 4.2e3    # specific heat capacity of working fluid
    Cpe::T = 5000.0   # effective specific heat of tank + working fluid
    S::T   = 800.0    # radiation flux
    hpf::T = 1000.0   # heat transfer coefficient (plate and working fluid)
    hpa::T = 100.0    # heat transfer coefficient (plate and atmosphere)
    hra::T = 100.0    # heat transfer coefficient (riser pipe and ambient)
    hda::T = hra      # heat transfer coefficient (down pipe and ambient)
    hta::T = 500.0    # loss coefficient (storage tank)    
    Ta::T  = 300.0    # ambient temperature
    T∞::T  = 295.0    # sky temperature
    α::T   = 5.5e-8   # radiation coefficient
    Af::T  = W*L*0.2  # transverse cross section area pipes in collector
    mdot::T = 0.5 * Af * ρf  # mass flow rate (assuming 0.5m per s in collector pipes
    A::T   = π*dt*ΔH  # ambient contact area per stratified layer
    V::T   = ΔH*π*(dt^2)/4  # tank volume per stratified layer
    ρCl::T = 1.0e3    
    mCr::T = 1.0e3   
    mCd::T = mCr
    ...
end

Subsystems

Collector

We consider a flat plate collector with parallel riser tubes modelled as an absorber plate transfering heat to the fluid in the risers. This model is based on the work of [4]. Heat transfer between plate and fluid modelled as convection using Newton's law of cooling [1].

Plate

We have the following processes

$$\rhop \delta C{pp}\frac{\partial T_p}{\partial t} = S + \delta k_p \frac{\partial^2 Tp}{\partial y^2}- h{pf}(T_p - Tf) - h{pa}(T_p - T_a) - \alpha(Tp^4 - T{\text{sky}}^4)$$

where $\alpha$ is a constant that includes the Stefan-Boltzmann constant, emissivities etc. This equation equation models heat transfer rate per unit area of the collector. Since we assume no dynamics along the width of the collector, we multiply the equation by the width $W$ of the collector to get the effective dynamics along the length of the conductor ($T_p = T_p(t, y)$).

$$W\rhop \delta C{pp}\frac{\partial Tp}{\partial t} = WS - Wh{pf}(T_p - Tf) - Wh{pa}(T_p - T_a) - W\alpha(Tp^4 - T{\text{sky}}^4)$$

Note that we only need the width factor to get the source term for the 1D advection equation used for simulating the working fluid in the collector.

$$\frac{\partial T_p}{\partial t} = \frac{1}{\rhop \delta C{pp}}\left( S + \delta k_p \frac{\partial^2 Tp}{\partial y^2}- h{pf}(T_p - Tf) - h{pa}(T_p - T_a) - \alpha(Tp^4 - T{\text{sky}}^4) \right)$$

We use Neumann boundary conditions:

$$\frac{\partial T_p}{\partial T_y} = 0\quad \text{at }y = 0, y = L$$

Working Fluid

We model heat transfer through and by the fluid using a 1D advection equation with the source term being the plate-fluid convective transfer term.

$$\rho_fAf C{pf}\frac{\partial T_f}{\partial t} + \rho_f Af C{pf} v_f\frac{\partial Tf}{\partial y} = Wh{pf}(T_p - T_f),$$

where $A_f$ is the transverse flow cross section. Note that mass conservation requires $\rho_f A_f v_f = \dot m$ to be constant.

$$\rho_f Ac C{pf}\frac{\partial Tf}{\partial t} + \dot{m} C{pf}\frac{\partial T_f}{\partial y} = Wh_pf(T_p - T_f),$$ $$\implies \frac{\partial T_f}{\partial t} = \frac{1}{\rho_f Ac C{pf}}\left(W h_{pf}(T_p - Tf) - \dot m C{pf} \frac{\partial T_p}{\partial y}\right)$$

Storage Tank

We use a stratified, well-mixed tank model as presented in [3,7].

(\rho C)_l V \frac{dT_{l,i}}{dt} = \dot{m} C_p(T_{l,i-1} - T_{l,i}) - h_{ta}A(T_{l,i} - T_a)
\implies \frac{dT_{l,i}}{dt} = \frac{1}{(\rho C)_l V}\left(\dot{m} C_p(T_{l,i-1} - T_{l,i}) - h_{ta}A(T_{l,i} - T_a)\right)

for $i \in 1 \ldots N_s$, where $Ns$ is the number of stratification layers, (model parameter). We let $T{l,0}(t) = T_r(t)$, where $T_t$ is the temperature of the up-riser pipe node connecting the collector outlet to the tank inlet (more on the connecting pipes later). Note that $T_l$ is the temperature of the same working fluid as in the collector, but we use a different subscript to distinguish between the two systems. This model is based on lecture 29 of [3], and a more detailed model is developed in [7]. We can write the system of equations for the storage tank in vector form

\dot{\overline{T}}_l = \frac{1}{(\rho C)_l}(D_l \overline{T}_l + d_l) - h_{ta}A(\overline T_l - T_a)

where

D_l\overline{T}_l + d_l = \begin{bmatrix}
    -1 \\
    1 & -1\\
      &  1 & -1\\
      &    &  1 & -1\\
      &    &    &   &\ddots\\
      &    &    &   &    1  & -1\\
      &    &    &   &       & 1  & -1
\end{bmatrix}
\begin{bmatrix}
T_{l,1}\\
T_{l,2}\\
T_{l,3}\\
T_{l,4}\\
\vdots\\
T_{l,N_s-1}\\
T_{l,N_s}
\end{bmatrix} +
\begin{bmatrix}
    T_r\\
    0\\
    0\\
    0\\
    \vdots\\
    0\\
    0
\end{bmatrix}

Connecting Pipes

We assume that the temperature drop across the connecting pipes is small and model them as single nodes between the collector and the tank, similar to as done in (TODO: cite). Denote $\text{d}$ as the down comer pipe that leads flow into the collector, and $\text{r}$ as the up-riser pipe leading from the collector to the tank. $Tr = T{l,0}$, and $T_d$ will be used in a similar manner with the collector when discretizing the equations.

(mC)_r\frac{dT_r}{dt} = \dot m C_pf(T_f(L) - T_r) - \pi d_r h_{ra}(T_r - T_a)
(mC)_d\frac{dT_d}{dt} = \dot m C_pf(T_{l,N_s} - T_d) - \pi d_dh_{da}(T_d - T_a)

Numerical Considerations

\frac{\partial T}{\partial y} \approx \frac{T_{i+1} - T_{i-1}}{2\Delta y} \quad \forall i \in 1,\ldots N_c

where $T_{i} = T(t, i\Delta y)$, and $N_c$ is the number of degrees of freedom in the discretized mesh. We also need to take into account the boundary conditions.

Plate

Consider first $T_p$.

\begin{bmatrix}
T'_{p,1}\\
T_{p,2}'\\
T_{p,3}'\\
\vdots\\
T_{p,N_c-1}'\\
T_{p,N_c}'\end{bmatrix} \approx \frac{1}{2\Delta y}\begin{bmatrix}
0 & 0 & 0  &       & &\\
-1 & 0 & 1 &       & &\\
  & -1 & 0 & 1     & &\\
  &   &   &        \ddots & &\\
  &   &   &        -1      & 0 & 1\\
  &   &   &        0      & 0 & 0  
\end{bmatrix}
\begin{bmatrix}
T_{p,1}\\
T_{p,2}\\
T_{p,3}\\
\vdots\\
T_{p,N_c-1}\\
T_{p,N_c}
\end{bmatrix}

Where the prime $'$ represents the spatial derivative, and the top and bottom rows are zero because of the Neumann conditions on the plate boundaries. Call this matrix $D_{p1}$ and write the discrete derivative as

\overline{T'}_p = D_{p1} \overline{T}_p

We need the second derivative for the conductive term as well.

\begin{bmatrix}
T_{p,1}''\\
T_{p,2}''\\
T_{p,3}''\\
\vdots\\
T_{p,N_c-1}''\\
T_{p,N_c}''
\end{bmatrix} \approx \frac{1}{\Delta y ^2}
\begin{bmatrix}
-2 & 2  & 0\\
1  & -2 & 1\\
   & 1  & -2 &  1\\
   &    &    &  \ddots &\\
   &    &    &         & 1 & -2 & 1\\
   &    &    &         &   &  2 & -2
\end{bmatrix}
\begin{bmatrix}
T_{p,1}\\
T_{p,2}\\
T_{p,3}\\
\vdots\\
T_{p,N_c-1}\\
T_{p,N_c}
\end{bmatrix}

Where we estimated the second derivative at the boundaries using the Neumann conditions.

\frac{T_{p,N_c+1} - T_{p,N_c-1}}{2\Delta y} = 0 \implies T_{p,N_c+1} = T_{p,N_c-1}
\implies T''_{p,N_c} \approx \frac{T_{p,N_c+1} + T_{p,N_c-1} - 2T_{p,N_c}}{\Delta y ^2} = \frac{2T_{p,N_c-1} - 2T_{p,N_c}}{\Delta y^2}

and similarly for $T''_{p,1}$.

Denote this matrix as $D_{p2}$ and use this compact form for discretized second derivative.

\overline{T''}_p = D_{p2}T_p

Putting everything together, we get the following system of ODEs for the discretized plate subsystem.

\dot{\overline{T}_p} = \frac{1}{\rho_p \delta C_{pp}}\left( S + \delta k_p D_{p2}\overline{T}_p - h_{pf}(\overline T_p - \overline T_f) - h_{pa}(\overline T_p - T_a) - \alpha(\lVert \overline T_p \rVert^4  - T_{\text{sky}}^4) \right)

Fluid in Collector

We can take $T_{f,0} = Td$ and $T{f,N_c+1} = T_r$ and use the central difference approximation as usual to get

\begin{bmatrix}
T'_{f,1}\\
T'_{f,2}\\
T'_{f,3}\\
\vdots\\
T'_{f,N_c-1}\\
T'_{f,N_c}
\end{bmatrix} \approx \frac{1}{2\Delta y}
\begin{bmatrix}
0 & 1\\
-1 & 0 & 1\\
  & -1 & 0 & 1\\
  &   &   &  \ddots\\
  &   &   &        & -1 & 0 & 1\\
  &   &   &        &   & -1 & 0\\
\end{bmatrix}
\begin{bmatrix}
T_{f,1}\\
T_{f,2}\\
T_{f,3}\\
\vdots\\
T_{f,N_c-1}\\
T_{f,N_c}
\end{bmatrix}
+
\frac{1}{2\Delta y}
\begin{bmatrix}
-T_d\\
0\\
0\\
\vdots\\
0\\
T_r
\end{bmatrix}

Again, compact notation

\overline{T}_f' = D_f\overline{T}_f + d_f

So that

\dot{\overline{T}}_f = \frac{1}{\rho_f A_c C_{pf}}\left(W h_{pf}(\overline{T}_p - \overline  T_f) - \dot m C_{pf} (D_f \overline T_f + d_f)\right)

Final System of Equations

\dot{\overline{T}_p} = \frac{1}{\rho_p \delta C_{pp}}\left( S + \delta k_p D_{p2}\overline{T}_p - h_{pf}(\overline T_p - \overline T_f) - h_{pa}(\overline T_p - T_a) - \alpha(\lVert \overline T_p \rVert^4  - T_{\text{sky}}^4) \right)
\dot{\overline{T}}_f = \frac{1}{\rho_f A_c C_{pf}}\left(W h_{pf}(\overline{T}_p - \overline  T_f) - \dot m C_{pf} (D_f \overline T_f + d_f)\right)
\dot{\overline{T}}_l = \frac{1}{(\rho C)_l}\left((D_l \overline{T}_l + d_l) - h_{ta}A(\overline T_l - T_a)\right)
(mC)_r\frac{dT_r}{dt} = \dot m(T_f(L) - T_r) - \pi d_r h_{ra}(T_r - T_a)
(mC)_d\frac{dT_d}{dt} = \dot m (T_{l,N_s} - T_d) - \pi d_dh_{da}(T_d - T_a)

The implementation can found here in the source.

Layout in Memory

We set ourselves up for contiguous memory accesses.

\begin{bmatrix}
\overline{T}_p\\
T_d\\
\overline{T}_f\\
T_r\\
\overline{T}_l
\end{bmatrix}

If we use an array/tensor programming language/library with GPU support (JuliaGPU, Halide, ArrayFire, TensorFlow etc), GPU acceleration should be simple to incorporate).

Solving

Once we have the final system of equations, we should be able to plug it into ODE and even steady-state solvers of our choice. We implement a simple fourth order explicit Runge-Kutta method.

function runge_kutta_4(f!, tspan, dt, prob::SWHSProblem)
    @unpack u, du = prob
    @unpack cache = prob
    T = eltype(u)
    k1 = zeros(T, length(u))
    k2 = zeros(T, length(u))
    k3 = zeros(T, length(u))
    k4 = zeros(T, length(u))
    t = first(tspan)
    while t < last(tspan)
        f!(k1, u, cache)
        f!(k2, u .+ dt .* k1 ./ 2, cache)
        f!(k3, u .+ dt .* k2 ./ 2, cache)
        f!(k4, u .+ dt .* k3, cache)
        @. u += dt/6 * (k1 + 2k2 + 2k3 + k4)
        t += dt
    end
    u
end

Results

Temperature of collector plate Temperature of fluid in collector Temperature distribution across tank

References

Limitations and Possible Refinements