Whole-Body Cryostimulation: New Insights in Thermo-Aeraulic Fields Inside Chambers

1. Introduction

Whole body cryotherapy (WBC) is a technique that entails subjecting one or multiple individuals to extremely cold and dry air within temperature-controlled cryochambers, for a duration ranging from one to four minutes [1]. Extensive research has demonstrated the various physiological, psychological, and physical advantages offered by cryotherapy [2,3,4]. Specifically, cryotherapy techniques have been developed to enhance post-exercise recovery, alleviate pain, and alleviate symptoms of depression and anxiety in patients with rheumatism and inflammatory conditions [5,6,7]. In October 2020, the International Institute of Refrigeration released its 39th Information Note on Refrigeration Technologies [8], which emphasized the need to establish a standardized method for measuring the actual temperatures inside whole-body or partial-body cryotherapy equipment. Furthermore, it is evident that there is a scarcity of well-established protocols and appropriate temperature guidelines to ensure optimal and safe cryotherapy treatments [9]. The scientific literature notably lacks comprehensive data regarding the exact temperature conditions within whole-body cryotherapy chambers, as well as the specific effects of multiple individuals being present during a cryotherapy session. This lack of information underscores the need for further research and investigation to enhance our understanding of these specific aspects. The objective of this study is to investigate the temperature distribution within a cryotherapy chamber as per the practices of cryotherapists. Specifically, the study aims to provide a precise understanding of the temperature distribution within empty chambers when the sole criterion for programming sessions is to reach the machine's set temperature, and whether the actual temperature inside the chamber matches the setpoint temperature. Additionally, the study investigates how the temperature field is affected by the presence of one or several people inside the chamber, which is a common scenario when the chamber size permits.

A mixed numerical and experimental approach is adopted, wherein simplified computational fluid dynamics simulations are used to model the global thermo-aerodynamic behavior inside empty and occupied cryotherapy chambers, and an experimental analysis using a network of thermocouples is conducted to determine the resulting 3D thermal field in real situations and provide answers to the research questions raised. To the best of our knowledge, this is the first study in the literature to investigate this topic in the context of whole-body cryotherapy.

2. Materials and Methods

The equipment utilized in this study was a dual WBC chamber (Mecotec^®, 2.10x1.84x2.15 m³) supplied with cold air from a cascade refrigeration system in a closed circuit. The cold air entered the chamber via perforated grids on the ceiling, consisting of 3132 holes of 13 mm diameter, whereas the evacuation occurred from a slotted extraction plate (110 slotted holes 105x20 mm²) located in the lower part of one of the walls with a flow rate of 1500 m³/h. A total of 15 thermocouples (type K 430-2000-2-1) were arranged in the volume, as depicted in Figure 1. and simultaneously connected to a data acquisition system to record local temperature changes. Furthermore, the thermocouples enabled volume mapping restitution through Abaqus^® software in conjunction with a finite element interpolation.

To illustrate how the enclosure heats up in the presence of a person, a simplified CFD steady-state model was developed. The walls are considered adiabatic and the dynamic inlet conditions are relative to a cold air flow of 1500 m³/h. The finite volume software ANSYS^® 2020 R2 was used, combined with a k-ε turbulence model. The mesh was considered structured, varying from 0.6x10⁶ cells (empty chamber) to 2.4x10⁶ cells (one male inside). A male 3D body (Artec^® scanner) was used and initialized at a constant surface temperature of 27°C corresponding to the mean value of male skin temperature during a WBC session [10,11].

Regarding the influence of persons considered as volume heated sources within the actual chamber, three young healthy participants (two females and one male), staff members of a cryocenter participated in this study. Corresponding anthropometric data are given in Table 1.

Body Surface Area and Body Volume were derived from the empirical formulas of Schumm et al. [12] and Sendroy & Collison [13] respectively. Moreover, the dissipated heat power per person was estimated at ambient temperature based on the energy balance equations [14]. Figure 1 illustrates the general positioning of the the individuals within the cryotherapy chamber. They are situated in the center of the chamber walls, approximately 30 cm away from the edges.

3. Results

3.1. Setpoint temperature vs actual temperature

The temperature mentioned (setpoint temperature) often refers to the temperature of the air at the outlet or in the vicinity of the expansion device. Unlike in cryosaunas (Partial Body Cryotherapy), where studies have demonstrated significant differences between the actual temperature inside the cryocabin and the temperatures reported by manufacturers [15,16], no research has been conducted on this topic for Whole Body Cryotherapy. To answer the question regarding the agreement between the setpoint temperature and the actual temperature in the middle of the cryochamber in WBC, Figure 2 presents the time evolution of both these fields. We began recording the two temperatures (setpoint and actual) when the setpoint temperature reached -94°C. It is evidenced that there is a strong temperature difference between the two curves over time. Once the steady state was reached, the difference was 18°C, corresponding to a 16% increase compared to the setpoint temperature. The main reason lies in the very high thermal gradient between the inside and the outside of the chamber, yielding thermal gains from walls, door, floor, and window as well as from thermal bridges at the mechanical junctions between walls and upper/ground floors. Although this study focuses on a specific technology, it can be assumed that other technologies face similar challenges.

3.2. Experimental 3D temperature fields

The experimental approach proved to be challenging to implement within the chamber due to the extreme temperatures and thermal limitations present, which render current electronic measurement devices incompatible. Consequently, we had to restrict our measurements to an actual temperature of -80°C. Figure 3 displays the temperature distribution inside the cryochamber over time for three scenarios, with and without participants, at various timestamps during a typical cryotherapy session of 3 min. The empty prisms in columns b and c indicate the approximate position of each participant. The results indicate that the presence of multiple occupants has a significant impact on the thermal fields inside the cryochamber. It should be noted that the hot thermal convective plumes are shifted towards the outlet, which is located at the bottom of the backside. This approach of accommodating several individuals together inside the cryochamber, driven by economic considerations, must be carefully considered in terms of response/dose balance. In other words, the duration of the protocol should be increased proportionally to the number of occupants present.

To better analyze the temperature distribution inside the chamber, two parameters have been introduced: the thermal gradient $\Delta T$ and the thermal heterogeneity percentage $TH\left(\%\right).$ The thermal gradient $\Delta T$ is defined as the difference between the maximum and minimum temperatures recorded inside the chamber, expressed as:

$$\Delta T=T_{max}-T_{min}$$

(1)

On the other hand, the thermal heterogeneity percentage $TH(\%)$ is defined as the absolute value of the ratio between the thermal gradient inside the whole studied volume and the minimum temperature. Specifically:

$$TH\left(\%\right)=100\left|\frac{\Delta T}{T_{min}}\right|$$

(2)

A completely homogeneous temperature field inside the chamber, where the maximum and minimum temperatures are equal, corresponds to $TH=0\%$, while $TH=100\%$ indicates that the maximum temperature is zero. The time-evolutions of the thermal gradient ΔT and thermal heterogeneity percentage are presented in Figure 4. It can be observed that an empty chamber exhibits a homogeneous thermal situation with limited stratification. However, as the number of occupants in the cryochamber increases, the magnitude of the thermal gradient (up to 10°C) and temperature heterogeneity (up to 13%) also increase. The observation of a stable thermal behavior with three participants, compared to a single participant where the two analyzed parameters increase over time, may be attributed to a more uniform spatial distribution of human heat sources within the volume.

4. Qualitative evidence of numerical thermo-aeraulic fields

4.1. Computational grid

The size of the mesh elements was meticulously determined, taking into account the criteria obtained from a comprehensive grid convergence study, which is extensively described in a recent article [17]. Firstly, a surface mesh is generated on the body, which is then followed by the creation of a Cartesian volumetric mesh for the fluid domain. In the computational domain, the cell size is set to be smaller than or equal to 0.023 m. The cell size assigned to the body ensures a high level of resolution for the boundary layer and ensures sufficient accuracy for thermal transfer and convection calculations within the air. The utilization of a structured Cartesian mesh contributes to improved convergence of results and helps to restrict the number of mesh cells. Overall, the mesh encompasses more than 3.5 million cells. To visually depict this, Figure 5 showcases the mesh structure on the body walls and in its immediate vicinity.

4.2. Numerical methods

In this investigation, the commercially available computational fluid dynamics (CFD) software ANSYS Fluent^® 2020 R2 was employed. This CFD solver, which utilizes the finite volume method, facilitates the solution of the governing equations for fluid flow. The numerical study conducted was three-dimensional, time-dependent, and non-isothermal in nature. For the resolution of pressure-velocity coupling, the widely used SIMPLE algorithm was employed, employing a first-order discretization scheme [18]. To model turbulence effects, the standard k-ε turbulence model was selected, which serves as a closure for the Reynolds-averaged Navier-Stokes equations. Heat transfer and convective mass transfer were represented using the following equations:

The continuity equation:

$$\frac{\partial w}{\partial t}+\frac{{\partial w}_x}{\partial x}+\frac{{\partial w}_y}{\partial y}+\frac{{\partial w}_z}{\partial z}=0$$

(3)

The Navier-Stokes equations:

$$\frac{{\partial w}_x}{\partial \tau }+w_x\frac{{\partial w}_x}{\partial x}+w_y\frac{{\partial w}_x}{\partial y}+w_z\frac{{\partial w}_x}{\partial z}=-\frac{1}{\rho }\frac{\partial P}{\partial x}+\nu \left(\frac{{{\partial }^{2}w}_x}{\partial x^2}+\frac{{{\partial }^{2}w}_x}{\partial y^2}+\frac{{{\partial }^{2}w}_x}{\partial z^2}\right)+g_x$$

(4)

$$\frac{{\partial w}_y}{\partial \tau }+w_x\frac{{\partial w}_y}{\partial x}+w_y\frac{{\partial w}_y}{\partial y}+w_z\frac{{\partial w}_y}{\partial z}=-\frac{1}{\rho }\frac{\partial P}{\partial y}+\nu \left(\frac{{{\partial }^{2}w}_y}{\partial x^2}+\frac{{{\partial }^{2}w}_y}{\partial y^2}+\frac{{{\partial }^{2}w}_y}{\partial z^2}\right)+g_y$$

(5)

$$\frac{{\partial w}_z}{\partial \tau }+w_x\frac{{\partial w}_z}{\partial x}+w_y\frac{{\partial w}_z}{\partial y}+w_z\frac{{\partial w}_z}{\partial z}=-\frac{1}{\rho }\frac{\partial P}{\partial z}+\nu \left(\frac{{{\partial }^{2}w}_z}{\partial x^2}+\frac{{{\partial }^{2}w}_z}{\partial y^2}+\frac{{{\partial }^{2}w}_z}{\partial z^2}\right)+g_z$$

(6)

The energy equation:

$${\rho }_f{C}_{p,f}\left(\frac{{\partial T}_f}{\partial \tau }+w_x\frac{{\partial T}_f}{\partial x}+w_y\frac{{\partial T}_f}{\partial y}+w_z\frac{{\partial T}_f}{\partial z}\right)=\frac{\partial }{\partial x}\left({\lambda }_f\frac{{\partial T}_f}{\partial x}\right)+\frac{\partial }{\partial y}\left({\lambda }_f\frac{{\partial T}_f}{\partial y}\right)+\frac{\partial }{\partial z}\left({\lambda }_f\frac{{\partial T}_f}{\partial z}\right)$$

(7)

With $w$ being the fluid velocity (m/s); $\nu$ the dynamic viscosity (kg/m.s); $x$,$y$,$z$ the spatial coordinates (m); $P$ the Pressure (Pa); $\rho$ the density (kg/m³); $T$ the temperature (K); $C_p$ the specific heat (J/Kg.K); $\lambda$ the thermal conductivity (W/m/K).

Polynomial laws were incorporated into the computational code to accurately capture the variations in thermo-physical properties of air with temperature. These laws enable the calculation of density, dynamic viscosity, thermal conductivity, and specific heat at each time step within the temperature range of -110°C to 33°C, which encompasses the relevant temperature range of interest.

The numerical injected temperature was -80°C, corresponding to the range later chosen to avoid experimental equipment failure due to extreme cold. Modeling the empty chamber using simplified adiabatic wall conditions reveals a 3D primary curvilinear flow stream from inlet to outlet, associated with a large-scale primary vortex cell and secondary spiral-type vortices at the top and side edges (see Figure 6). To evidence both flow patterns and convective effects, the CFD results illustrate how heat released by a person (mean skin temperature 27°C) situated in the center of the room is transported to the exhaust vents. The resulting 3D complex convective flow results in local temperature heterogeneities, maybe leading to both thermal stratification and non-whole-mixing process which are experimentally to be verified.

Figure 7 depicts the temperature fields associated with three subjects in the CCE chamber at t=180 s. Notably, the man's temperature is observed to be higher compared to the two women. When comparing the two women, it can be observed that the woman positioned near the aeraulic outlet experiences a more pronounced cooling effect than the woman in the middle. By examining both Figure 6 and Figure 7, it becomes evident that the flow topology plays a crucial role in influencing the cooling process of the female subject situated on the left side of the Figure. This is attributed to the acceleration of air velocity at the outlet, which enhances heat exchange and subsequently promotes skin cooling, particularly in the lower extremities.

Figure 8 presents a comparison of the time-dependent evolution of average skin temperature for a single individual during a 180-second continuous cooling exposure (CCE) session. The experimental data is contrasted with the numerical simulation, highlighting the numerical method's capability to accurately replicate the heat transfer processes between the human body and its surrounding environment. These findings underscore the reliability and validity of the numerical approach in studying thermal dynamics in human subjects.

5. Conclusions

The present pilot study reveals two major points. Firstly, it appears necessary for a reflection to be initiated among the community of cryogenic chamber users (manufacturers, physicians, cryotherapists, researchers) regarding the concordance between set temperatures and the temperatures actually reached in the chamber volume. Secondly, the common practice of having multiple people in the chamber simultaneously requires an awareness of the induced thermal modifications, leading to an appropriate adaptation of protocol duration, particularly when the analgesic thermal threshold is the targeted value of cutaneous temperatures. The results showed, for example, that an empty chamber exhibits a homogeneous thermal situation with limited stratification, whereas as the number of occupants in the cryochamber increases, the magnitude of the thermal gradient (up to 10°C) and temperature heterogeneity (up to 13%) also increase. Future studies should target this current topic.