Finite Element Analysis of Heat and Mass Exchange Between Human Skin and Textile Structures

. The study is aimed at the development of finite element models for the analysis of numerical physical systems, applying them to the simulation of heat and mass transfer processes occurring between human skin and textile structures. The proposed methodology and numerical analysis of textile structures can be applied to predict the physical properties and performance of future products at the design stage. For example, the finite element models can be used in the development of protective clothing, outdoor clothing, passive and active cooling systems. In this work, finite element analysis (FEA) was done by using Comsol Multiphysics and Matlab software.


Introduction
It is widely accepted that the main function of protective clothing is to protect the human body from a dangerous environment. According to (Li, 2001) this protection covers many functions such as maintaining the right thermal environment to the body, preventing the body from being injured by abrasion, radiation, wind, chemical substances. This study focuses on thermal comfort between human skin and textiles. According to (Tessier, 2017), (Umair et al., 2016), the most essential components of thermal comfort are air permeability and thermal resistance.
In recent decades, three-dimensional textile structures are widely used as moisture and thermal regulating layer in different applications. For example, automotive textiles, outdoor textiles, protective textiles, as well as medical bandages, mattresses (Orlik et al., 2018). Many researchers developed and investigated various aspects of heat and mass transfer through textile structures. (Zhang et al., 2022) developed a personal microclimate management system for maintaining thermal comfort. The clothing system consisted of several micro-fans placed in the garment side seam to provide cooling air. The computational fluid dynamics (CFD) was used to simulate the influence of the fan's number and air gap distance on a macro-scale. It was found that the convection effect is enhanced by the increase in the fan's numbers, but the fans' cooling effect varies a lot because of various air gap distances (Zhang et al., 2022). (Shen et al., 2021) created a torso thermal manikin model to solve the heat and flow transfer through cold protective clothing under various ambient conditions. The simulations were done using the CFD approach. The simulations allowed researchers to investigate the effect of environmental temperature and clothing thickness on heat and mass transfer. (Siddiqui and Sun, 2017) proposed conjugated heat transfer finite element model to predict the effective thermal conductivity and thermal resistance of plain weft knitted fabric on a micro-scale. Finite element analysis was performed by using Abaqus/CAE software. (Kangro, 2018) analytically and numerically investigated 1-D non-stationary diffusion-convection problem in a layered domain. During the study, the conservative averaging method that was reduced to the initial value problem of ODEs using the integral exponential type splines was applied. The study was performed using Matlab software. (Puszkarz and Machnowski, 2020) numerically studied the safety and thermal comfort of protective clothing used by firefighters. Two multilayer structures that are used in thermal protective clothing were analyzed by using Solidworks Flow Simulation software. The simulation allowed to predict heat transfer index and heat transmission factor. The outcomes demonstrated a good agreement with experimental data. In research work (Puszkarz and Krucińska, 2018) computational techniques were used to determine the air permeability coefficient of the double-layer knitted fabrics. It was found that the air permeability of the knitted fabric strongly depends on the thickness and geometrical parameters of the yarn. (Angelova et al., 2011) numerically analyzed the transverse permeability of a textile woven structure. The flow in the through-thickness direction of the woven structures was depicted as jet systems. The simulations were performed using Fluent software, which allows solving the set of Reynolds averaged Navier-Stokes equations.
To sum up, modern finite element computing technologies provide a highly realistic representation of the physical processes. In most cases, the outcomes of computer simulations are sufficiently adequate and reliable enough to allow simulations to take the place of many field experiments. In this study, finite element models are developed to predict air permeability and thermal resistance coefficients of different textile structures including 3D textile on a micro-scale. The computational simulations were performed using Comsol Multiphysics and Matlab software.

Model development
Computational models such as finite element models facilitate an understanding of physical processes that are often difficult to achieve in other ways. For example, it is difficult to measure temperature due to the complex internal structure of the textile.
According to (Datta and Rakesh, 2009) the essential step of modeling is a physical problem formulation into a mathematical analog (an equivalent mathematical formulation). The general flowchart of the essential computational modeling steps is shown in Fig. 1. Finite element analysis consists of three main steps: pre-processing, processing, and post-processing. Typically, the required data on what to solve and how to solve in finite element software is covered in the pre-processing step (Datta and Rakesh, 2009). In the processing step, the governing partial differential equations are modified into a system of algebraic equations and the unknown values (such as temperature, velocity, etc.) are calculated. Finally, the post-processing step includes visualizing the solution obtained from the previous step (Datta and Rakesh, 2009). A flowchart for model development of pre-processing, processing, and postprocessing steps are depicted in Fig. 2. For all of these steps, the finite element software COMSOL multiphysics will be used in this investigation. However, due to some limitations in the post-processing step, the final results will be displayed using Matlab software.

Post-processig
In this step, the set of algebraic equations is solved. That include error estimation and convergence check.
During pre-processing, the required data is entered into the FE software.

Select solver and compute
Evaluation of solution using visualization (plotting), data export tools, etc. Acceptance of the results and final display of results.

A representative volume element of the textile structure
In this work, 3D steady-state finite element models were developed to evaluate air permeability (AP) and thermal resistance (R ct ) coefficients of different textile structures. A representative volume element (RVE) was created according to (Zupin et al., 2011) air permeability measurements of one-layer textile structures see Table 1. The air permeability experiment was done with the Air Permeability tester FX 3300 Labotester III (Textest Instruments).
In this study, two plain weave samples were examined in order to recreate an RVE of a one-layer textile. More details are presented in our previous work (Gadeikytė et al., 2022). Fig. 3. Illustrates different RVE that consists of the fluid and textile domain. Model x_a and Model x_c indicate a single-layer textile structure, Model x_b and Model x_d represent two-layer textile structures. Furthermore, Model x_e depicts simplified 3D textile geometry. It should be noted, that in fluid flow modeling, the properties of textile material do not influence the computational solution.

Governing equations of air permeability modeling
The numerical prediction of the air permeability (AP) coefficient through textile was done under the assumption that the fluid flow is incompressible Newtonian single-phase flow. In order to evaluate the air permeability (AP) coefficient of textile structure, the set of Navier-Stokes equations (Eq. 1) with the continuity equation (Eq. 2) were set in the fluid domain (Ω 1 ), and Brinkman equations (Eq. 3) with the continuity equation (Eq. 2) were applied in the textile domain (Ω 2 ). Following (Elgeti and Sauerland, 2014), Eq. 1 consists of the inertia (convection) term ( ⋅ ) , divergence of the stress − + ( + ( ) ), and external forces F. It should be noted, it was assumed that there are no external forces. Furthermore, the term −∇ denotes the pressure difference forces.
Where porosity, K local permeability tensor in m 2 , the term is defined as effective viscosity, ρfluid (air) density, uvelocity. The Brinkman equation (Eq. 3) requires the local permeability K coefficient. The local permeability of a porous medium can be determined according to the analytical formulas proposed by (Gebart, 1992) or from well known Kozeny-Carman equation. In this work, the local permeability K and porosity values were used from literature (Pezzin, 2015).
The boundary conditions ∂Ω that were applied in the finite element model were based on ISO 9237:1995 (E) standard. The difference between inlet and outlet boundaries was set to 200 Pa according to (Zupin et al., 2011) measurement. A more detailed analysis of boundary conditions was presented in our previous work (Gadeikytė et al., 2022), Table 2, and Fig. 4.

Figure 4. Boundary conditions of a unit cell of three-dimensional textile
In this study, the coarse mesh was automatically applied for the approximated geometry domain (Model x_c, Model x_d, and Model x_e) and the coarser mesh was set to close to the real geometry (Model x_a and Model x_b). The set of Navier-Stokes and Brinkman equations was solved by using the stationary nonlinear solver that uses Newton's method combined with a direct PARDISO solver. The first-order discretization for velocity and pressure approximation was applied.
For validation purposes, the air permeability model was simplified and compared with the analytical solution proposed by (Gebart, 1992). The analytical transverse local permeability was obtained according to Eq. 4.
where Rradius of fiber, fiber volume fraction and c denotes equivalent Kozeny constant (Gebart, 1992).  Outlet. Relative pressure 0 = 0 Pa on Ω outlet surface was set on. Physical meaning is that absolute pressure is equal to the reference pressure (1 [atm]). It can be expressed as: Numerical simulations were done by investigating 2D finite element model considering that fibers are impermeable, the Reynolds number is low (Re << 1), there are no external forces, no inertia effect. Furthermore, the Navier-Stokes momentum and mass conversation equations (see Eq. 1, Eq. 2) are reduced to the set of Stokes equations (see Eq. 5) (Geoffre et al., 2020).
The boundary conditions of Stokes flow are demonstrated in Fig. 5. The fluid properties were selected according to (Geoffre et al., 2020), (Nabovati et al., 2010) investigations. The local permeability coefficient was determined by applying Darcy's law. This part of the study can be used to determine the local permeability coefficient K that requires the Brinkman equation. The findings of the local permeability coefficient are shown in Fig. 5. It illustrates a good agreement between the analytical solution and the simplified case of the AP model. It should be noted, that the relative error does not exceed 4 %.

Governing equations of thermal resistance modeling
The experimental thermal resistance R ct coefficient can be obtained according to ISO 11092:2014 standard. According to (Tessier, 2017) ISO 11092 is often referred to as the "skin model." It imitates the heat transfer through the skin. In this work, the numerical modeling of thermal resistance measurement was based on Sweating Guarded Hot Plate M259B method. In this method, the experiment of the R ct coefficient requires that the temperature of the measuring unit be set to 35 °C, the ambient air temperature -20 °C, and relative humidity -65%.
In Comosol Multiphysics software, R ct coefficient can be investigated by using Heat transfer in solids and fluids (.ht) interface. The governing equations of the model are based on energy conservation (see Eq. 6), where C p is specific heat capacity, Q stands for the overall heat transfer, = − , where k denotes the thermal conductivity of the fluid-solid mixture (COMSOL, 2018).
The numerical model was done under the assumption that the flow is a stationary nonisothermal flow. Furthermore, the heat is transferred due to conduction. The fluid is air and the textile is made from polyester. The boundary conditions of thermal resistance are summarized in Table 3. The model's geometry with boundary conditions is illustrated in Fig. 6. Furthermore, during the simulation normal mesh size was automatically applied in all heat transfer simulations. Moreover, R ct coefficient might be obtained from Eq. 7. The equivalent formula in Comsol environment is (aveop1(T)-aveop2(T))/aveop1(ht.ntflux), where a) b) aveop1(), aveop2() are the average operators on Ω inlet and Ω outlet , ht.ntflux stands for normal total heat flux.
Where Q stands for the heat flow, Athe surface area, and ∆ is the temperature difference.

Numerical Results
The post-processing analysis was done using Matlab software. The main results of finite element simulations are average air permeability (AP) and thermal resistance R ct coefficients on a micro-scale. Furthermore, distributions of air velocity and temperature flow through textile structures. The AP coefficients through 3D textile are depicted in Fig. 7. The summary of obtained AP coefficients is shown in Table 4. Column Error indicates relative error between experimental data (Zupin et al., 2011) and numerical simulations. It was found that the relative error was less than 7.53 % when close to the real geometry (Model 2_a and Model 3_a) was used. 3D textile models (Model 2_e and Model 3_e) demonstrate that the resistance to air flow increases by adding a spacer yarn.    Figure 9. The temperature distribution of Model 2_a (solid matrix-polyester, fluid matrix-air), and on the right, the element size distribution in µm Moreover, it was found that thermal resistance R ct coefficients of 3D textile structures made from polyester are 0.112 and 0.109 respectively to Model 2_e and Model 3_e. According to literature (Barauskas and Abraitiene, 2011), the measured thermal resistance coefficient is 0.118 (K·m 2 )/W of 3D textile made from polyester. However, the main reason for different results is the different thickness of the sample. Also, 3D textile structures can be made from different types of yarn, including bioceramic additives. The thermal distributions are illustrated in Fig. 8 and Fig. 9.

Conclusions
The proposed methodology allows predicting the air permeability (AP) and thermal resistance (R ct ) coefficients on a micro-scale. These finite element models can be used in the development of protective clothing in the early design stage.
The finite element models demonstrated a good agreement with experimental data found in the literature. In addition, the models can be easily expanded to include various configurations of geometries or environmental conditions.