# Numerical modeling of a Vuilleumier heat pump and validation using CFD reference data

I Barakat, E Albin, S Xin

<figure><img src="image_1.png" /></figure>

## To cite this version:

I Barakat, E Albin, S Xin. Numerical modeling of a Vuilleumier heat pump and validation using CFD reference data. Turbulence, Heat and Mass Transfer, Jul 2025, Tokyo (Japan), Japan. hal-05188198

HAL Id: hal-05188198

https://hal.science/hal-05188198v1

Submitted on 26 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical modeling of a Vuilleumier heat pump and validation using CFD reference data

I. Barakat¹, E. Albin² and S. Xin¹

¹INSA Lyon, CNRS, CETHIL, UMR5008, 69621 Villeurbanne, France, israa.barakat@insa-lyon.fr

²Universite Claude Bernard Lyon 1, CETHIL, UMR5008, CNRS, INSA Lyon, Villeurbanne, F-69100, France, eric.albin@univ-lyon1.fr

**Abstract** – Vuilleumier heat pumps are promising thermal machines for energy-efficient heating and cooling applications, as they operate without requiring much electrical work input. This study presents a numerical model of a Vuilleumier heat pump based on a one-dimensional approach that solves the conservation equations of mass, momentum, and energy while the working fluid is a perfect gas. The model employs an implicit time-stepping scheme and finite volume spatial discretization to ensure numerical stability and accuracy. The model is validated by comparing its results to reference data obtained from Computational Fluid Dynamics (CFD) simulations under operating conditions documented by Dogkas et al. (2019b). The comparative analysis shows strong agreement between the one-dimensional model and CFD predictions, confirming the reliability of this simplified modeling approach. This work highlights the potential of the numerical model to support the design and optimization of Vuilleumier heat pumps. It provides an efficient tool for exploring innovative configurations and improving performance, contributing to the advancement of sustainable energy solutions.

## 1. Introduction

Regenerative thermal machines are in sharp focus, being important means to optimize energy transfers in minimizing energy losses, and addressing current energy challenges. They do provide a good means of converting thermal energy into mechanical work or useful heat with a thermodynamic cycle. Therefore, they are being applied to power generation or refrigeration. Among them, the Stirling engines and Vuilleumier heat pumps (VHPs) stand out for their efficiency and adaptability for work with a number of various energy sources, including renewables.

VHPs were invented by Vuilleumier (1918) and are characterized by a novel thermodynamic cycle which enables both heating and cooling with a low energy input. Using inert gases such as helium, these systems cover wide ranges of temperatures and are especially well suited for residential heating and low-temperature refrigeration applications.

A wide variety of modeling approaches have been proposed for regenerative thermal machines: from analytical zonal models (Hargreaves, 1991; Formosa and Despesse, 2010; Ahmadi et al., 2016), to more refined numerical models (Rogdakis et al., 2019; Boroujerdi et al., 2011; Wang et al., 2016; Chen et al., 2014). The analytical modeling gave certain theoretical representation about the thermodynamic cycles, although most were less accurate due to vast simplification. Numerical models, especially finite volume method-based numerical models, allow for improved detail regarding the internal flows and heat transfer.

Dogkas et al. (2019b) used CFD-based work to analyze VHP performance, but these methods are still very expensive computationally and very difficult to apply for optimizing new configurations.

In this study, we develop a one-dimensional (1D) numerical model based on the finite volume method to simulate VHP operation. The model has been validated through experiments on Stirling engines, showing good agreement with real-world results, making it suitable for the study of Vuilleumier Heat Pumps (Barakat et al., 2025). Building upon the work of Wang et al.
<figure><img src="image_3.png" /><figcaption>Figure 1: Geometrical representation of the Vuilleumier heat pump.</figcaption></figure>

(2016) and Boroujerdi et al. (2011), our model incorporates an improved treatment of geometry discontinuities and thermal resistances at interfaces. The numerical predictions are validated against theoretical results from Dogkas et al. (2019b,a). This study aims to explore innovative configurations of regenerative thermal machines and optimize their design, particularly for trigeneration applications that integrate heating, cooling, and power generation.

## 2. Model description

A CFD simulation of a β-type Vuilleumier Heat Pump was conducted by Dogkas et al. (2019b). The corresponding geometry, representing 1/8th of the machine, is illustrated in Figure 1a. In this study, the same geometry is modeled using a transient numerical tool to capture its dynamic behavior, employing a 1D discretization approach, as schematically depicted in Figure 1c.

The numerical tool is based on the approach proposed by Wang et al. (2016) and Boroujerdi et al. (2011), with several refinements to enhance both accuracy and computational efficiency. A particular treatment is applied at the interface between components or at the transition through the porous matrix to handle section jumps accurately. Figure 1b illustrates an abrupt change in gas flow cross-section at interface $f_i$, separating cells i and i-1 ($S_{l,f_i}^g \neq S_{r,f_i}^g$ with $l = left$ and $r = right$). This section change also leads to additional heat transfer surfaces ($A_{f_i}^{o,w}$ and $A_{f_i}^{g,w}$), which are nonzero at the interface between the external surface (o = outside) and the solid wall (w = wall).

### 2.1 Wall Discontinuity Management

Unlike conventional approaches in the literature, the wall and the porous matrix (m for matrix) are treated separately, each governed by independent equations.

In the absence of discontinuities, the energy conservation in the wall is expressed by a semi-discrete equation (space discretized, time continuous):

$$ \frac{\partial(m_i^w c_{v,i}^w T_i^w)}{\partial t} = k_{f_i+1}^w S_{f_i+1}^w \frac{T_{i+1}^w - T_i^w}{\delta x_{f_i+1}} - k_{f_i}^w S_{f_i}^w \frac{T_i^w - T_{i-1}^w}{\delta x_{f_i}} + \sum_{j \in \{g,m,o\}} \alpha_i^{w,j} A_i^{w,j} (T_i^j - T_i^w) \quad (1) $$

A second-order implicit time discretization is then applied ($\frac{\partial \Psi_i}{\partial t} \approx \frac{a_k}{\delta t} \Psi_i^{k+1} + \frac{a_{k-1}}{\delta t} \Psi_i^{k-1} + \frac{a_{k-2}}{\delta t} \Psi_i^{k-2}$),
leading to a tridiagonal system to solve:

$$a_{1,i}^w \cdot T_{i-1}^{w,k} + a_{2,i}^w \cdot T_i^{w,k} + a_{3,i}^w \cdot T_{i+1}^{w,k} = b_i^w \quad (2)$$

$$ \text{with} \quad \left\{ \begin{aligned} a_{1,i}^w &= -\frac{k_{fi}^{w,k} S_{r,fi}^w}{\delta x_{fi}^k} \\ a_{2,i}^w &= \frac{a_k}{\delta t} \cdot m^w c_v^w + \frac{k_{fi+1}^{w,k} S_{l,fi}^w}{\delta x_{fi+1}^k} + \frac{k_{fi}^{w,k} S_{r,fi}^w}{\delta x_{fi}^k} + \sum_{j \in \{g,m,o\}} \alpha_i^{w,j,k} A_i^{w,j} \\ a_{3,i}^w &= -\frac{k_{fi+1}^{w,k} S_{l,fi}^w}{\delta x_{fi+1}^k} \\ b_i^w &= -\frac{a_{k-1}}{\delta t} (m^w c_v^w T^w)_i^{k-1} + \frac{a_{k-2}}{\delta t} (m^w c_v^w T^w)_i^{k-2} + \sum_{j \in \{g,m,o\}} \alpha_i^{w,j,k} A_i^{w,j} T_i^{j,k} \end{aligned} \right. $$

The resolution of this equation, coupled with those of the gas and porous matrix temperatures, is performed iteratively until convergence. This equation is not solved in heaters and coolers where wall temperatures are imposed to be constant as in Dogkas et al. (2019b).

The discontinuity at interface $f_i$ introduces an additional thermal resistance in conduction. To account for this, a fictitious temperature cell $T_{fi}^w$ is introduced. The energy balances on cells $i, f_i$, and $i-1$ are written as :

$$ \frac{\partial m_{i-1}^w c_{v,i-1}^w T_{i-1}^w}{\partial t} = \frac{T_{fi}^{w,k} - T_{i-1}^{w,k}}{R_{l,fi}} - k_{fi-1}^{w,k} S_{fi-1}^w \frac{T_{i-1}^{w,k} - T_{i-2}^{w,k}}{\delta x_{fi-1}} + \sum_{j \in \{g,m,o\}} \alpha_{i-1}^{w,j,k} A_{i-1}^{w,j} (T_{i-1}^{j,k} - T_{i-1}^{w,k}) \quad (3a) $$

$$ \frac{\partial m_{fi}^w c_{v,fi}^w T_{fi}^w}{\partial t} = \frac{T_i^{w,k} - T_{fi}^{w,k}}{R_{r,fi}} - \frac{T_{fi}^{w,k} - T_{i-1}^{w,k}}{R_{l,fi}} + \sum_{j \in \{g,m,o\}} \alpha_{fi}^{g,w,k} A_{fi}^{g,w} (\theta_{fi}^{g,k} - T_{fi}^{w,k}) \quad (3b) $$

$$ \frac{\partial m_i^w c_{v,i}^w T_i^w}{\partial t} = k_{fi+1}^{w,k} S_{fi+1}^w \frac{T_{i+1}^{w,k} - T_i^{w,k}}{\delta x_{fi+1}} - \frac{T_i^{w,k} - T_{fi}^{w,k}}{R_{e,fi}} + \sum_{j \in \{g,m,o\}} \alpha_i^{w,j,k} A_i^{w,j} (T_i^{j,k} - T_i^{w,k}) \quad (3c) $$

By incorporating convection resistances $R_{l,fi} = \frac{\delta x_{i-1}/2}{k_{i-1} \cdot S_{i-1}^w} + \frac{|\ln r_{moy}/r_{i-1}|}{2\pi k_{i-1} \cdot e}$ and $R_{r,fi} = \frac{\delta x_i/2}{k_i \cdot S_i^w} + \frac{|\ln r_i/r_{moy}|}{2\pi k_i \cdot e}$ with $r_{moy} = \frac{r_{i-1}+r_i}{2}$ and the thickness $e$ of the discontinuity. The vertical wall at $f_i$ exchanges heat with the gas at temperature $T_i^g$ if the discontinuity is increasing, at $T_{i-1}^g$ if it is decreasing, and does not exchange heat if no discontinuity is present. The heat transfer coefficient is approximated by $\alpha_{fi}^{g,w} = (\alpha_{i-1}^{g,w} + \alpha_i^{g,w})/2$.

The balance (3b) is used to express the averaged temperature of the discontinuity based on the temperatures of adjacent cells $i$ and $i-1$:

$$ T_{fi}^{k,w} = \frac{\frac{T_i^{w,k}}{R_{r,fi}} + \frac{T_{i-1}^{w,k}}{R_{l,fi}} + N_{fi}^w}{D_{fi}^w} \quad (4) $$

$$ \text{with} \quad \left\{ \begin{aligned} D_{fi}^w &= \frac{1}{R_{l,fi}} + \frac{1}{R_{r,fi}} + \sum_{j \in \{g,m,o\}} \alpha_{fi}^{g,w,k} A_{fi}^{g,w} + \frac{a_k}{\delta t} (m^w c_v^w)_fi^k \\ N_{fi}^w &= \sum_{j \in \{g,m,o\}} \alpha_{fi}^{g,w,k} A_{fi}^{g,w} \theta_{fi}^{g,k} - \frac{a_{k-1}}{\delta t} (m^w c_v^w T^w)_fi^{k-1} - \frac{a_{k-2}}{\delta t} (m^w c_v^w T^w)_fi^{k-2} \end{aligned} \right. $$

By injecting the expression (4) of $T_{fi}^w$ into equations (3a) and (3c), $T_{fi}^w$ is eliminated to maintain the same dimension of the linear system to be solved, whether a discontinuity is present or not. The balance at cells $i-1$ and $i$ can be formulated as modified tridiagonal systems:

$$ a_{1,i-1}^{w,*} \cdot T_{i-2}^{w,k} + a_{2,i-1}^{w,*} \cdot T_{i-1}^{w,k} + a_{3,i-1}^{w,*} \cdot T_i^{w,k} = b_{i-1}^{w,*} \quad \text{et} \quad a_{1,i}^{w,*} \cdot T_{i-1}^{w,k} + a_{2,i}^{w,*} \cdot T_i^{w,k} + a_{3,i}^{w,*} \cdot T_{i+1}^{w,k} = b_i^{w,*} \quad (5) $$

with
$$
\left\{
\begin{array}{l}
a_{1,i-1}^{w,*} = a_{1,i-1}^{w} \\
a_{2,i-1}^{w,*} = a_{2,i-1}^{w} - \frac{\delta x_{fi}}{k_{fi} S_{fi}^{w}} + \frac{1}{R_{i,fi}} \left( 1 - \frac{1}{R_{i,fi}} \frac{1}{D_{fi}^{w}} \right) \\
a_{3,i-1}^{w,*} = a_{3,i-1}^{w} + \frac{\delta x_{fi}}{k_{fi} S_{fi}^{w}} - \frac{1}{R_{r,fi}} \frac{1}{R_{i,fi}} \frac{1}{D_{fi}^{w}} \\
b_{i-1}^{w,*} = b_{i-1}^{w} + \frac{1}{R_{i,fi}} \frac{N_{fi}^{w}}{D_{fi}^{w}}
\end{array}
\right.
\quad \text{and} \quad
\left\{
\begin{array}{l}
a_{1,i}^{w,*} = a_{1,i}^{w} + \frac{\delta x_{fi}}{k_{fi} S_{fi}^{w}} - \frac{1}{R_{r,fi}} \frac{1}{R_{i,fi}} \frac{1}{D_{fi}^{w}} \\
a_{2,i}^{w,*} = a_{2,i}^{w} - \frac{\delta x_{fi}}{k_{fi} S_{fi}^{w}} + \frac{1}{R_{r,fi}} \left( 1 - \frac{1}{R_{r,fi}} \frac{1}{D_{fi}^{w}} \right) \\
a_{3,i}^{w,*} = a_{3,i}^{w} \\
b_{i}^{w,*} = b_{i}^{w} - \frac{1}{R_{r,fi}} \frac{N_{fi}^{w}}{D_{fi}^{w}}
\end{array}
\right.
$$

The coefficients of the tridiagonal system are corrected as indicated in (5) for the cells adjacent to the discontinuities at faces $f_i$.

## 2.2 Discretization of the other conservation equations

The porous matrix does not present any discontinuities, and the energy conservation is expressed by a semi-discrete equation (6) which is similar to that of wall (1). The temperatures of the matrix are then solved using the tridiagonal system (2) where the index w is replaced by m.

$$ \frac{\partial m_i^m c_{v,i}^m T_i^m}{\partial t} = \left[ k^m S^m \frac{\partial T^{m,k}}{\partial x} \right]_{f_i}^{f_{i+1}} + \sum_{j \in \{g,w,o\}} \alpha_i^{m,j,k} A_i^{m,j} (T_i^{j,k} - T_i^{m,k}) \quad (6) $$

The helium mass varies in each cell $i$ according to the equation (7a) where $[\dot{m}]_{fi}^{f_{i+1}} = \dot{m}_{fi+1}^k - \dot{m}_{fi}^k$ represents the difference in gas mass flow rates at each face, ensuring mass conservation in the engine. The conservation of momentum and energy in the gas is described by the equation (7b) and (7c). The friction in the cell $f_i$ is calculated by $C_{fi}^F \cdot V_{fi} = C_{fi}^{F,l}V_{fi}^l + C_{fi}^{F,r}V_{fi}^r$ to distinguish the friction in the left and right volumes at the discontinuities: $V_{fi}^l = \pi r_{i-1}^2 \cdot \delta x_{i-1}/2$ and the right volume $V_{fi}^r = \pi r_i^2 \cdot \delta x_i/2$. The friction coefficients $C_{fi}^F$ and the convective heat transfert $\alpha_i^{g,j}$ are defined by correlations described in (Boroujerdi et al., 2011; Wang et al., 2016).

$$ \frac{\partial m_i^g}{\partial t} + [\dot{m}]_{fi}^{fi+1} = 0 \quad (7a) $$

$$ \frac{\partial (m^g u)_{fi}}{\partial t} + [\dot{m}u]_{i-1}^i = (p_i - p_{i-1}) \frac{V_{fi}}{\delta x_{fi}} + C_{fi}^F V_{fi} \dot{m}_{fi} \quad (7b) $$

$$ \frac{\partial m_i^g \left( \frac{u_i^2}{2} + c_{v,i} T_i^g \right)}{\partial t} + \left[ \dot{m} \left( \frac{u_i^2}{2} + c_P^g T^g \right) \right]_{fi}^{fi+1} = \left[ k^g S^g \frac{\partial T^g}{\partial x} \right]_{fi}^{fi+1} + p_i \frac{dV_i}{dt} + \sum_{j \in \{m,w,o\}} \alpha_i^{g,j} A_i^{g,j} (T_i^j - T_i^g) \\ + \alpha_{fi}^{g,w,j} A_{fi}^{g,w} (T_{fi}^w - T_i^g) \delta_{fi}^+ + \alpha_{fi+1}^{g,w,j} A_{fi+1}^{g,w} (T_{fi+1}^w - T_i^g) \delta_{fi+1}^- \quad (7c) $$

The ideal gas law describes the thermodynamic state of helium and allows the system to be closed:

$$ p_i = m_i^g r T_i^g / V_i^g \quad (8) $$

The equation system (7) is solved with a second-order staggered scheme for the velocity flow ($\dot{m}u$) and enthalpy ($\dot{m}c_p T^g$):

$$ (\dot{m}u)_i^k = \max(u_i^k, 0) (G_{1,i} \dot{m}_{fi}^k - G_{2,i} \dot{m}_{fi-1}^k) + \min(u_i^k, 0) (G_{3,i} \dot{m}_{fi+1}^k - G_{4,i} \dot{m}_{fi+2}^k) \quad (9a) $$

$$ (\dot{m}h)_{fi}^{k} = \max(\dot{m}_{fi}^{k} c_{p,fi}^{k}, 0) (G_{1,fi} T_{i-1}^{g,k} - G_{2,fi} T_{i-2}^{g,k}) + \min(\dot{m}_{fi}^{k} c_{p,fi}^{k}, 0) (G_{3,fi} T_{i}^{g,k} - G_{4,fi} T_{i+1}^{g,k}) \quad (9b) $$

The geometric interpolation coefficients are:

$$ G_{1,i} = \frac{\delta x_i^k + 2\delta x_{i-1}^k}{2\delta x_{i-1}^k}, \quad G_{2,i} = \frac{\delta x_i^k}{2\delta x_{i-1}^k}, \quad G_{3,i} = \frac{\delta x_i^k + 2\delta x_{i+1}^k}{2\delta x_{i+1}^k}, \quad G_{4,i} = \frac{\delta x_i^k}{2\delta x_{i+1}^k} \quad (10a) $$

$$ G_{1,fi} = \frac{\delta x_{fi}^k + 2\delta x_{fi-1}^k}{2\delta x_{fi-1}^k}, \quad G_{2,fi} = \frac{\delta x_{fi}^k}{2\delta x_{fi-1}^k}, \quad G_{3,fi} = \frac{\delta x_{fi}^k + 2\delta x_{fi+1}^k}{2\delta x_{fi+1}^k}, \quad G_{4,fi} = \frac{\delta x_{fi}^k}{2\delta x_{fi+1}^k} \quad (10b) $$
The mass conservation (7a) is coupled with momentum (7b) through the ideal gas law (8). After applying a second-order time discretization, it results in a pentadiagonal system to be solved:

$$a_{1,f_i}^g \cdot \dot{m}_{f_i-2}^k + a_{2,f_i}^g \cdot \dot{m}_{f_i-1}^k + a_{3,f_i}^g \cdot \dot{m}_{f_i}^k + a_{4,f_i}^g \cdot \dot{m}_{f_i+1}^k + a_{5,f_i}^g \cdot \dot{m}_{f_i+2}^k = b_{f_i}^g \quad (11)$$

with

$$\left\{
\begin{aligned}
a_{1,f_i}^g &= G_{2,i-1} \cdot \max(u_{i-1}^k, 0) \\
a_{2,f_i}^g &= -G_{1,i-1} \cdot \max(u_{i-1}^k, 0) - G_{2,i} \cdot \max(u_i^k, 0) - \frac{V_{f_i}^k}{\delta x_{f_i}^k} \frac{\delta t}{a_k} \frac{r T_{i-1}^k}{V_i^{k,g}} \\
a_{3,f_i}^g &= -G_{3,i-1} \cdot \min(u_{i-1}^k, 0) + G_{1,i} \cdot \max(u_i^k, 0) + a_k \frac{\delta x_{f_i}^k}{\delta t} + \frac{V_{f_i}^k}{\delta x_{f_i}^k} \frac{r \delta t}{a_k} \left( \frac{T_{i-1}^k}{V_i^{k,g}} + \frac{T_{i-1}^k}{V_i^{k,g}} \right) - C_{f_i}^F V_{f_i}^{g,k} \\
a_{4,f_i}^g &= G_{4,i-1} \cdot \min(u_{i-1}^k, 0) + G_{3,i} \cdot \min(u_i^k, 0) - \frac{V_{f_i}^k}{\delta x_{f_i}^k} \frac{\delta t}{a_k} \frac{r T_{i-1}^k}{V_i^{k,g}} \\
a_{5,f_i}^g &= -G_{4,i} \cdot \min(u_i^k, 0) \\
b_{f_i}^g &= \left\{ \left( \frac{a_{k-1}}{a_k} \frac{p_{i-1}^{k-1} V_{i-1}^{k-1}}{T_{i-1}^{k-1}} + \frac{a_{k-2}}{a_k} \frac{p_{i-1}^{k-2} V_{i-1}^{k-2}}{T_{i-1}^{k-2}} \right) \frac{T_{i-1}^k}{V_i^{k-1}} + \left( \frac{a_{k-1}}{a_k} \frac{p_{i-1}^{k-1} V_{i-1}^{k-1}}{T_{i-1}^{k-1}} + \frac{a_{k-2}}{a_k} \frac{p_{i-1}^{k-2} V_{i-1}^{k-2}}{T_{i-1}^{k-2}} \right) \frac{T_{i-1}^k}{V_i^{k-1}} \right\} \cdot \frac{V_{f_i}^{g,k}}{\delta x_{f_i}^k}
\end{aligned}
\right.$$

Compared to literature, the term in pressure gradient factor ($V_{f_i}/\delta x_{f_i}$) is modified with $V_{f_i} = V_{f_i}^l + V_{f_i}^r$. Without discontinuities, this term is equivalent to the one passage used in the literature.

After time discretization, the energy conservation (7c) leads also to a pentadiagonal system:

$$a_{1,i}^g T_{i-2}^{g,k} + a_{2,i}^g T_{i-1}^{g,k} + a_{3,i}^g T_{i}^{g,k} + a_{4,i}^g T_{i+1}^{g,k} + a_{5,i}^g T_{i+2}^{g,k} = b_{i}^g, \quad (12)$$

$$\left\{
\begin{aligned}
a_{1,i}^g &= G_{2,f_i} \cdot \max(\dot{m}_{f_i}^k c_{p,f_i}^k, 0) \\
a_{2,i}^g &= -G_{1,f_i} \cdot \max(\dot{m}_{f_i}^k c_{p,f_i}^k, 0) - G_{2,f_{i+1}} \cdot \max(\dot{m}_{f_{i+1}}^k c_{p,f_{i+1}}^k, 0) - \frac{k_{f_i}^{g,k} S_{r,f_i}^{g}}{\delta x_{f_i}^k} \\
a_{3,i}^g &= -G_{3,f_i} \cdot \min(\dot{m}_{f_i}^k c_{p,f_i}^k, 0) + G_{1,f_{i+1}} \cdot \max(\dot{m}_{f_{i+1}}^k c_{p,f_{i+1}}^k, 0) + \frac{k_{f_i}^{g,k} S_{r,f_i}^{g}}{\delta x_{f_i}^k} + \frac{k_{f_{i+1}}^{g,k} S_{t,f_i}^{g}}{\delta x_{f_{i+1}}^k} + \frac{a_k}{\delta t} m^g c_v^g + \sum_{j \in \{m,w,o\}} (\alpha_i^{g,j})^k A_i^{g,j} \\
a_{4,i}^g &= G_{4,f_i} \cdot \min(\dot{m}_{f_i}^k c_{p,f_i}^k, 0) + G_{3,f_{i+1}} \cdot \min(\dot{m}_{f_{i+1}}^k c_{p,f_{i+1}}^k, 0) - \frac{k_{f_{i+1}}^{g,k} S_{t,f_{i+1}}^{g}}{\delta x_{f_{i+1}}^k} \\
a_{5,i}^g &= -G_{4,f_{i+1}} \cdot \min(\dot{m}_{f_{i+1}}^k c_{p,f_{i+1}}^k, 0) \\
b_i^g &= -\left(\dot{m}_{\frac{u^2}{2}}\right)_{f_{i+1}}^k - \left(\dot{m}_{\frac{u^2}{2}}\right)_{f_i}^k - \frac{a_k}{\delta t} \left(\frac{m^g u^2}{2}\right)_{i}^k + \frac{a_{k-1}}{\delta t} \left(m^g c_v T^g + \frac{m^g u^2}{2}\right)_{i}^{k-1} + \frac{a_{k-2}}{\delta t} \left(m^g c_v T^g + \frac{m^g u^2}{2}\right)_{i}^{k-2} + \sum_{j \in \{m,w,o\}} (\alpha_i^{g,j})^k A_i^{g,j} T_{i}^{j,k}
\end{aligned}
\right.$$

If the system does not present any discontinuities ($\delta_{f_i}^+ = \delta_{f_{i+1}}^- = 0$), it remains unchanged; However, with increasing discontinuity ($\delta_{f_i}^+ = 1$) at the face $f_i$ or decreasing discontinuity ($\delta_{f_{i+1}}^- = 1$) at the face $f_{i+1}$, it is necessary to modify the coefficients $a_{3,i}^g$ and $b_i^g$:

$$\left\{
\begin{aligned}
a_{3,i}^g &= a_{3,i}^g + \alpha_{f_i}^{g,w} A_{f_i}^{g,w} \delta_{f_i}^+ + \alpha_{f_{i+1}}^{g,w} A_{f_{i+1}}^{g,w} \delta_{f_{i+1}}^- \\
b_i^g &= b_i^g + \alpha_{f_i}^{g,w} A_{f_i}^{g,w} \delta_{f_i}^+ T_{f_i}^w + \alpha_{f_{i+1}}^{g,w} A_{f_{i+1}}^{g,w} \delta_{f_{i+1}}^- T_{f_{i+1}}^w
\end{aligned}
\right.
\quad (13)$$

### 3. Results of Vuilleumier Heat Pump Modeling

The operating conditions of the simulation are set to impose a mean pressure of 44 bar and wall temperatures respectively of 500 and 0°C for hot and cold heaters, and 40°C for coolers, with an operating frequency of 41.68 Hz. The simulation ran for approximately 2000 cycles, which enables to achieve convergence of the mechanical and heat powers.
## 3.1 Temperature

In the hot and cold heaters, the simulation results show that the gas temperature generally remains lower than the wall temperature throughout most of the cycle, leading to heat absorption by the gas, as shown in Figure 2. In contrast, within the hot and cold expansion spaces, the gas temperature occasionally exceeds the heater temperature, reflecting the dynamic interaction between heat transfer and gas movement.

<figure><img src="image_4.png" /><figcaption>Figure 2: Gas temperature evolution in the hot/cold expanders and at the center of the hot/cold heaters: CFD vs. 1D transient model (with and without enhanced convection).</figcaption></figure>

Figure 2a shows the gas temperature profiles in the hot heater and expander, while Figure 2b displays the corresponding behavior in the cold section. Our transient 1D model captures the general trends of the CFD results by Dogkas et al. (2019b), but some amplitude differences are observed, especially in the expansion spaces.

To analyze these differences, we introduce a third curve representing the same transient model, but with the convective heat transfer coefficient multiplied by a factor of 2 throughout the entire machine. This enhanced-convection case exhibits closer agreement with the CFD results in both the hot and cold sections, particularly in matching the peak values.

The improvement suggests that the CFD simulation by Dogkas, which uses the $k-\epsilon$ turbulence model with a coarse mesh, results in more intense heat exchange between the gas and the wall. In our model, the convective heat transfer coefficients are computed using classical empirical correlations, with separate formulations for laminar and turbulent flows. While these correlations are widely used in 1D models, they may underestimate the actual convective exchanges under oscillating flow conditions. As Dogkas pointed out in his PhD thesis, “The experimentally derived equations for the Nusselt of oscillating flow failed to match with CFD results for all heat exchangers because the conditions of the experiments differ from those appearing in the tested Vuilleumier machine” (Dogkas, 2019). This supports the explanation of the observed differences.

Without experimental data, it is difficult to assert which approach is more realistic. However, this analysis provides valuable insight into how convective terms influence the model accuracy.

## 3.2 Pressure

Figure 3a shows the evolution of the four variable volume spaces of the machine during the cycle: hot expansion (①→②), cold expansion (②→③), cold compression (③→④), and hot compression (④→①).

In addition to these spaces, the total compression space is also presented, which is the sum of the hot and cold compression volumes. The total gas volume remains constant throughout the cycle, which is a key aspect of the machine's operation. However, but the gas moves from one space to another, impacting heat transfers and pressures.
Figure 3b presents the numerical pressure-volume diagrams for the three variable-volume spaces. These diagrams illustrate the internal work exchanges within the system: the enclosed areas of the total compression space (i.e., the sum of the hot and cold compression volumes) represent the work required to compress the gas, while the areas associated with the hot and cold expansion spaces correspond to the work produced during gas expansion.

Our 1D transient results were compared with the CFD simulations from Dogkas et al. (2019a), revealing an overall satisfactory trend. However, a noticeable discrepancy appears in the cold expansion space, where our model exhibits lower pressure losses compared to the CFD results. This difference can be attributed to an underestimation of dissipative effects in the 1D model.

After introducing a friction coefficient multiplied by a factor of two, the agreement with the CFD results improved significantly, particularly in the cold and hot parts of the cycle. This highlights the importance of accurately accounting for frictional effects, which appear to be more dissipative in the CFD approach, especially in the cold expansion zone.

<figure><img src="image_5.png" /><figcaption>Figure 3: Pressure and volume evolution during the converged cycle.</figcaption></figure>

## 4. Conclusion

A one-dimensional transient model was developed to simulate the operating processes of a Vuilleumier heat pump (VHP). The model is based on the finite volume method and includes detailed treatment of the main thermodynamic and heat transfer phenomena, particularly at component interfaces where geometric discontinuities occur.

It incorporates classical empirical correlations to evaluate convective heat transfer and frictional effects in tubes and porous media.

A comparative analysis with the CFD results from Dogkas et al. (2019b,a) demonstrated the ability of the model to reproduce the global thermal and dynamic behavior of the machine.

However, some discrepancies were observed. They were reduced by adjusting the convective and friction coefficients, revealing the sensitivity of the 1D model to these empirical parameters.

Overall, the proposed approach offers a reliable and computationally efficient tool for analyzing and optimizing VHPs. It also paves the way for future investigations into advanced configurations systems for heating, cooling, or power production.
Acknowledgements

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under ANR-21-CE05-0007 (project DAVinCCHeaP). They are also grateful to Dr. Dogkas for responding to our requests.

References

Ahmadi, M. H., Ahmadi, M. A., Pourfayaz, F., Bidi, M., Hosseinzade, H., and Feidt, M. (2016). Optimization of powered stirling heat engine with finite speed thermodynamics. *Energy Conversion and Management*, 108:96–105.

Barakat, I., Albin, E., and Xin, S. (2025). Simulation numérique du stirling gpu3 avec gestion améliorée des sauts de section. In *Actes du Congrès de la Société Française de Thermique (SFT)*.

Boroujerdi, A. A., Ashrafizadeh, A., and Naeenian, S. M. M. (2011). Numerical analysis of stirling type pulse tube cryocoolers. *Cryogenics*, 51(9):521–529.

Chen, W.-L., Wong, K.-L., and Chang, Y.-F. (2014). A computational fluid dynamics study on the heat transfer characteristics of the working cycle of a low-temperature-differential gamma-type stirling engine. *International Journal of Heat and Mass Transfer*, 75:145–155.

Dogkas, G. (2019). *Thermodynamic Analysis of Vuilleumier Heat Pumps*. PhD thesis, National Technical University of Athens, Athens. Supervisor: Emmanouil Rogdakis.

Dogkas, G., Bitsikas, P., Tertipis, D., and Rogdakis, E. (2019a). Vuilleumier machine speed-effect investigation with cfd and analytical model. *Applied Thermal Engineering*, 162:114314.

Dogkas, G., Rogdakis, E., and Bitsikas, P. (2019b). 3d cfd simulation of a vuilleumier heat pump. *Applied Thermal Engineering*, 153:604–619.

Formosa, F. and Despesse, G. (2010). Analytical model for stirling cycle machine design. *Energy Conversion and Management*, 51(10):1855–1863.

Hargreaves, C. M. (1991). *The Philips Stirling Engine*. Elsevier.

Rogdakis, E., Bitsikas, P., Dogkas, G., and Antonakos, G. (2019). Three-dimensional cfd study of a β-type stirling engine. *Thermal Science and Engineering Progress*, 11:302–316.

Vuilleumier, R. (1918). Method and apparatus for inducing heat changes.

Wang, K., Dubey, S., Choo, F. H., and Duan, F. (2016). A transient one-dimensional numerical model for kinetic stirling engine. *Applied Energy*, 183:775–790.