Functionally graded optimisation of adsorption systems with phase change materials

Adsorption phenomena are encountered in several engineering applications. One of its uses is in the storage and transport of gas in the form of adsorption tanks. The exothermic nature of the adsorption process decreases adsorption capacity presenting an impetus to understand the thermal characteristics of the gas storage process. Studies using mixtures of phase change materials and adsorbents in adsorption tanks demonstrate potential improvements in the adsorption capacity of the tanks. They also show that the distribution of phase change materials and adsorbent are important. Thus, this work presents two approaches for optimising the adsorbent domain. The first is to use a semi-analytic model to determine the best homogeneous material concentration for the adsorbent and phase change material for the vessel composition. The other is to use a 2D axisymmetric model to perform FGM optimisation to distribute material in the tank. Results for both models are presented and discussed for different conditions. The study shows that, for the cylindrical geometry, FGM optimisation is always, at least, marginally better than the homogeneous distribution from the semi-analytic model. However, FGM optimisation demands more computing time increases the complexity of implementation and results assembling. The semi-analytic approach is a possible alternative for optimising adsorption systems with phase change material mixed with adsorbents.


Introduction
Compressed natural gas (CNG) and liquefied natural gas (LNG) require high pressures or low temperatures (cryogenic), respectively. They are commonly used methods used for the storage and transport of gas (Hirata et al. 2009). An alternative means of storing and transporting gas is the use of the adsorbed natural gas (ANG), which makes it possible to store gas under lower pressures (35 to 50 atmospheres) and ambient temperature, when compared to CNG and LNG respectively, thereby, demanding less energy (Hirata et al. 2009). However, ANG storage density, which is lower than CNG, and the difficulty to retrieve the stored gas from an adsorption system, are challenges that ANG technology presents.
Gas adsorption is exploited in many industrial applications; examples include gas transport and storage (adsorbed natural gas), in the industrial separation process (Cavalcante 2000), in the CCS process (Leung et al. 2014) and in cooling application (Askalany et al. 2012). Adsorption consists of the adhesion of molecules from a fluid phase to the adsorbent surfaces. When the fluid reaches the adsorbed state, it becomes denser and sticks to the adsorbent surface, which is why adsorption may be used to store gas; as the process is exothermic, it releases heat. In contrast, desorption is an endothermic process and takes heat from the system. Figure 1 shows an schematic of this process. ANG technology uses the natural gas adsorption in porous materials to store and transport gas (Ayawei et al. 2017), thereby bringing benefits such as operating in relatively low pressures, which can result in lower operating costs and low energy consumption (Ayawei et al. 2017).
Two important characteristics of the adsorption systems that are related to the performance of the equipment are the cycle time required to complete its charging and discharging cycles and the amount of residual gas left inside the vessel after a full cycle. These two characteristics are affected by the thermal behaviour of the system. Adsorption phenomena liberate heat; however, when the adsorbent material temperature increases, its adsorption capacity decreases. The opposite is the case for desorption, as temperature drops, so does the desorption capacity. The Dubinin-Astakov relation (Sahoo et al. 2010) describes the influence of the temperature in the saturation state. Consequently, thermal management has a significant influence on an adsorbent system as it impacts the adsorption and desorption capacity. The impact varies according to system ability to dissipate heat via thermal conduction and convection. The reduction of temperature changes during the cycles allows the same system to retain the same amount of gas for a shorter cycle time or store more gas during the same cycle time. In this work, phase change materials are used as a mean of passive thermal management.
Phase change materials (PCMs) are substances which have a high phase change enthalpy and, therefore, are considered suitable for storing energy in their internal energy (Sharma et al. 2009). Their elevated phase change enthalpy permits the material to absorb significant amounts of heat at a constant temperature; as a consequence, PCMs are a potential means of passively decreasing Fig. 1 Adsorption and desorption phenomena sketch. Adsorption liberates heat when the adsorbed state is reached, and desorption takes heat to reach gaseous state the temperature difference between the adsorption and desorption cycles of an ANG vessel.
The literature presents models for the adsorption phenomena, most of which involve carbon as adsorbent (Hirata et al. 2009;Goetz and Biloé 2005;Rahman et al. 2010;Sahoo et al. 2011;Sahoo et al. 2010;Mota et al. 2004;Mota 1999;Mota et al. 1997;Sahoo and Ramgopal 2014;Rahimi et al. 2015;Vo et al. 2016;Gautam 2020). The adsorbed gas stored inside the vessel is usually represented by the portion of pores occupied by the gas in the adsorbed phase. The Navier-Stokes equations may be used to solve the fluid problem (Wenger et al. 2008;Sahoo et al. 2010;Vasiliev et al. 2012); however, in gas adsorption problems, and for the porosity range presented in this study, the result from Darcy's law equation provides a good correlation to the volume-averaged Navier-Stokes equation due to the decrease of the inertial and viscous terms (Sahoo et al. 2011). Additionally, Darcy's law eliminates approximation problems regarding the LBB conditions (Polak 1983;Biloé et al. 2002;Song et al. 2010). There are 0D approaches known to present acceptable results for adsorption problems . Similarly to adsorption, there has been a range of models of the application, modelling and simulations involving phase change materials (Sharma et al. 2009;Zhang et al. 2007;Farid et al. 2004;Tyagi and Buddhi 2007;Mondal 2008;Fořt et al. 2017;Riahi et al. 2016;Kozak and Ziskind 2012).
The literature contains experimental results related to phase change materials improving adsorbent vessels. Studies involve the placement of PCM inside various adsorbent systems including for methane storage (Li and Li 2015), CO 2 storage (Toledo et al. 2013) and gas filtering (Zimmermann and Keller 2006). The means of positioning and geometry of the PCM varies from spheres (Toledo et al. 2013), copper tubes (Li and Li 2015). In all reported cases, PCM addition brings improvements related to higher storage/filtering capacity due to the decrease of the thermal amplitude during the operating cycles. The study shows that the phase change material amount and position inside the domain affect the results (Toledo et al. 2013). This suggests a distribution problem that arises from the fact that adding PCM reduces the area available to absorb. When the PCM amount and position are correct, it brings some improvements, and it is verified that in a filtering system, a mixture presenting 25% of PCM increased 17% the amount of gas filtered by the equipment (Zimmermann and Keller 2006). In this work, the employment of a Functionally Graded Material (FGM) distribution is expected to increase the capacity of the store and deliver gas by the adsorption storage systems evaluated.
FGM is related to the concept of topology optimisation (Paulino and Silva 2005). In this work, the topology optimisation method is employed to generate a graded solution for the adsorption domain with PCM. Topology optimisation algorithms are versatile tools applied in many fields and allow the optimised distribution of materials to be determined. Heat transfer problems are one of them, and topology optimisation algorithms have been developed to act in such systems; for instance, the optimisation of twodimensional heat conduction problems (Gersborg-Hansen et al. 2005;Zhuang et al. 2007;Yoon 2010;Turteltaub 2001), multi-physics thermo-fluid problems (Koga et al. 2013;Matsumori et al. 2013), three-dimensional heat conduction problems (Burger et al. 2013), and transient heat transfer problems (Zhuang et al. 2013). In adsorption systems, it is verified that topology optimisation presented improvements when distributing different materials inside the system . The reliability of the obtained solutions was also verified experimentally (Sá et al. 2018).
The system studied in this paper is a cylindrical adsorption vessel with a mixture of adsorbent and phase change material mixed in its interior. This type of continuous mixture of PCM and adsorbent is possible as both may be particulate (Hawlader et al. 2003;Su et al. 2017). Figure 2 shows a cylindrical vessel filled with adsorbent and PCM spheres. These spheres are mixed in the interior and adsorbent stores gas while PCM spheres keep the temperature low. The external walls have contact with the ambient air and lose heat via convection.

Fig. 2
Sketch of an adsorption system studied in this paper. It presents a mixture of adsorbent and phase change material inside its domains An adsorption vessel performance may be improved, by changing its external surrounding, boundary conditions, and geometry . Such vessels may also be enhanced internally by changing material distribution, composition, and any apparatus that aids in thermal management inside the domain. This study focuses on the internal approach, the employment of a (i) 0D semianalytic and a (ii) 2D functionally graded distribution will present a homogeneous solution and a graded solution, respectively, both of which aim to increase the capacity to store and deliver gas. Here, a 0D semi-analytic model for an adsorption system with an adsorbent-PCM mixture in its interior is developed. For the 2D and semi-analytic model, the phase change is modelled using the apparent specific heat capacity method for PCM phase change phenomena. The optimisation of the semi-analytic model is performed by a full factorial sweep varying the ambient temperature, and the phase change temperature. In contrast, the optimisation of the 2D model is performed via FGM. Results from both models are compared. Figure 3 shows the sketch of the 0D homogeneous solution and a graded solution for this system. The graded solution considers the effects caused by boundaries and heat transfer transport mechanisms, whereas the homogeneous solution neglects them. The grey zones determine the portion of phase change material mixed with the adsorbent.
The 2D approach employs a single variable to describe the portion of PCM at each finite element to simulate the PCM behaviour inside the vessel, resulting in the modelling of a porous domain, composed of both PCM and adsorbent. This variable is responsible for adjusting the terms of the governing equations linearly between PCM and adsorbent. The adsorption tank described in Sahoo et al. (2011) is modelled to validate the implemented physics and be optimised. The implementation for PCM is validated using the Stefan problem analytical solution.
The optimised PCM distribution is expected to improve the adsorption vessel capacity to store and deliver gas. Results will determine the optimised PCM ratio in a mixture that reduces the thermal amplitude during the operating cycles and maximise the mass of stored and delivered gas. Additionally, the FG optimisation will show how different an optimised homogeneous PCM distribution is from an optimised FG PCM distribution. This information may be useful when deciding whether it is convenient to develop technology to build such vessels.
This paper conducts the study as follows. Section 2 describes the formulation and assumptions adopted. Furthermore, Section 3 presents the optimisation problem regarding the phase change material distribution inside an adsorption system. Section 4 illustrates the details about the implementation and the variational formulation employed for finite element method solving. Section 5 discusses the Fig. 3 The sketch of homogeneous distribution and an FGM distribution inside an adsorption system results obtained for the proposed cases. Finally, Section 6 concludes the work and presents the research highlights.

Mathematical model
This section describes the physics of the considered problem. Section 2.1 presents the basis assumptions for the model. Section 2.2 presents the governing equations for the adsorption phenomenon isolated for the 2D model. Section 2.3 describes the phase change phenomenon. Section 2.4 combines the equations for adsorption and PCM showing the governing equations for the 2D system with a mixture of adsorbent and PCM inside the domain. Section 2.5 presents the semi-analytic model equations derived by removing the space integral from the 2D model equations.

Assumptions
In order to build mathematical relations for this problem, the following assumptions are made: 1. Adsorbent constant specific heat capacity: the range of temperature fluctuation is expected to be around 100 K for the full carbon vessel and about 10 K for the vessel with PCM; consequently, it is suitable to consider the specific heat constant for an optimisation . Additionally, there is a good correlation between experiment and simulation results when the specific heat is constant (Sahoo et al. 2011. 2. Local thermal equilibrium: it is assumed that the gas, the adsorbent material, and the PCM are in local thermal equilibrium (Mota et al. 2004). For a given position inside the domain, the gas, PCM, and adsorbent can be described by a single temperature. 3. Ideal gases: gases are assumed to exhibit their ideal thermodynamic behaviour (Sahoo et al. 2011); consequently, the equation that expresses density is: where M g is the molecular gas mass, R g is the universal gas constant, P is the gas pressure, and T is the gas temperature. 4. Linear Driving Force (LDF) adsorption model: the Linear Driving Force model is chosen to represent the film diffusion, which is the dominant effect in the adsorption phenomena studied in this work (Mota et al. 1997). 5. Constant adsorption enthalpy H : is assumed to be a constant due to the temperature and pressure ranges as confirmed by Sahoo et al. (2011). 6. Thermal conductivity: the Maxwell-Eucken relation is adopted to describe the thermal conductivity (k eff (ε t )) as a porous domain (Smith et al. 2013): where k s and k g are the thermal conductivity for the solid and gas, respectively, and ε t is the material total porosity. 7. Gas velocity: the gas velocity, u g , can be approximated by the Darcy law due to the low velocity within in the porous media as in ; therefore, it can be expressed by: where μ g is the gas viscosity and K f is the medium permeability. Inertial and viscous forces are assumed to be negligible. 8. Encapsulated PCM: As the PCM is encapsulated, there is no mass exchange between the PCM and the adsorbent (Hassabou 2011). The shell is assumed to be soft; consequently, the pressure inside the spheres is equal to the adsorption domain internal pressure. This assumption also implies that the enclosure can accommodate the density variation between one phase and the other with no pressure increases due to the expansion resistance imposed by the shell, thereby undergoing a phase change at constant pressure. 9. Constant properties for PCM: the specific heat, thermal conduction, density, and phase change enthalpy are considered constant for the phase change material in this work (Hassabou 2011). As observed in the performed analysis, the temperature variation during the simulations is narrow. Additionally, as the system is expected to transition between solid to liquid, the pressure variation will not cause considerable variation in the phase change temperature and thermophysical properties as explained in chapter 6 in Sandlel et al. (2006). Therefore, PCM properties are assumed to be constant during the process.
These general assumptions are used to assemble the basis for the 2D model and the semi-analytic model. However, as the semi-analytic model is obtained by removing the spatial derivatives from the 2D model, the assumptions that concern the spatial gradients such as thermal conduction are not applied to the semi-analytic model. The semianalytic model considers the domain to have uniform characteristic properties. Table 1 illustrates for which model each assumption is made. Most of them are applied to both models; the exceptions are thermal conductivity and gas velocity. These find no meaning in the semi-analytic model because it does not use the spatial derivatives.

Physics formulation for adsorption in the 2D model
This subsection presents the transport equations for the porous adsorbent material used in the 2D model. The state variables for the governing equations as P (pressure), T (temperature), and Q (adsorbed gas mass ratio). The mass equation transport may be written as : where ε t (ads) = ε b (ads) + (1 − ε b (ads) )ε p (ads) is the porosity accessible to the gas phase, composed of the microporosity ε b (ads) and the macroporosity ε p (ads) , u g is the gas velocity, ρ g is the gas density, ρ b = (1 − ε t (ads) )ρ s is the bed density, ρ s is the adsorbent density, and Q is the ratio of adsorbed gas mass to adsorbent mass ( kg kg ). The first term describes the gas present in the macropores and the adsorbed gas in the micropores while the second term describes the gas flow in the location. Figure 4 shows the physical representation of macro and microporosity. Macroporosity is related to the spaces between the spheres while microporosity is related to the micropores in the particles.
The energy balance can be written as : where C pg is the specific heat of gas, and C eff (ads) is the effective specific heat for the adsorbent. C eff (ads) is written as : The LDF model (Mota et al. 1997), which describes the adsorption rate, ∂Q ∂t can be written as: where Q eq is the equilibrium density of the adsorbed gas, given by the Dubinin-Astakhov (DA) equation (Sahoo et al. 2011): where E 0 is the characteristic energy of adsorption, β is an affinity coefficient, and n f is an equation parameter related to the micropore dispersion. G is given by (Vasiliev et al. 2012): where D is the intra-particle diffusion time constant. E [a,d] is either the adsorption activation energy (subscript 'a') or the desorption activation energy (subscript 'd') and depends on whether there is local adsorption or desorption taking place. If δQt >= 0 then E [a,d] = E a and E [a,d] = E d otherwise.
The adsorbed gas density, ρ adg , is given by Sahoo et al. (2011): where α e is the mean value of the thermal expansion of liquified gases and ρ boil is the gas density at boiling point temperature (T boil ). The effective volume of micro-pores, W f , is given by: where ρ f is the adsorbent density and w f is the specific volume of micropores. The adsorption potential, A, is given by Sahoo et al. (2011): where the saturated pressure P sat is given by Sahoo et al. (2011): where P cr is the critical pressure and T cr is the critical temperature of the gas.

Formulation of phase change phenomenon in the 2D model
This section presents the phase change material model for the 2D axisymmetric model. The method for modelling the phase change phenomena for this application is based on the apparent heat capacity (Civan and Sliepcevich 1987). Although there are two other methods (enthalpybased and temperature-based), however, the apparent heat method is chosen because the primary variable, which is the temperature, can be obtained directly from the solution of the finite element problem, and the semi-analytic problem. Furthermore, it does not require the separate solution of single-phase equations for each of the phases. The limitation of this method is the singularity that arises when the phase change interval is set to zero, and a range is not admissible (Civan and Sliepcevich 1987). Therefore, this method is characterised by a small temperature range where the phase change phenomenon occurs, Hassabou (2011), which is 3K in this work. The energy equation, for a pure substance undergoing a phase change, as shown in image 32, may be written as Civan and Sliepcevich (1987): where the subscripts 1 and 2 refer to phase 1 and 2 respectively, ρ is density, h is enthalpy, and φ 1 is the volume fraction (volume of one phase divided by the volume occupied by both phases).
Equation (14) may be simplified by applying the PCM constant properties assumption. Thus, (14) may be rewritten as (Civan and Sliepcevich 1987): where: where the subscript designates each phases properties, one refers to phase 1 and 2 to phase 2, L pc is the phase change enthalpy. The φ 1 is similar to a Heaviside function, and it assumes value 1 when the temperature is above phase change temperature and value 0 otherwise. In this work, a hyperbolic tangent function is chosen to ensure a continuous behaviour for the phase change physics. Therefore, φ 1 becomes: where γ pcm is a factor which adjust the curve behaviour, here assuming the value 1.4.
The chosen method to model phase change assumes the phenomenon takes place in a small temperature range. The derivative of φ 1 may be written as: Figure 5 presents the plots for φ and ∂φ ∂T . In this work, the phase change material is modelled as a porous body composed of many spheres. The same porous medium approach, which is applied for the adsorbent, is applied for the PCM. The gas may not be adsorbed by phase change material; however, it may flow through the spaces between spheres; therefore, the mass balance for the PCM becomes: ∂(ε t (pcm) where ε t (pcm) = ε bpcm + (1 − ε b(pcm) )ε p (pcm) and describes the total porosity brought by the spaces between the PCM spheres, ρ g is the gas density, and u g is the gas velocity. The first term represents the air between the spheres and the second term represents the gas stream that travels through these spaces. Equation (15) also can be written as the energy balance of this porous medium.
where C eff (pcm) is written as: where k eff (pcm) is (2) applied to PCM properties.

Formulation of the mixed adsorbent-PCM domain in the 2D model
In this subsection, the adsorption model and the PCM model are coupled using the design variable ε pcm . The design variable defines the element composition; it assumes value 0 when no PCM is present and 1 when the element is composed 100% of PCM. The equations terms are adjusted according to the design variable. Figure 6 illustrates the role of the design variable in determining the tank.
Summing (4) and (20), the mass equation for the domain becomes: The LDF model remains the same, and no adjustments are made.

Semi-analytic model
The semi-analytic model is obtained from (23) and (24) with spatial derivatives set to zero. It considers the whole system as a control volume and neglects local phenomena such as thermal conduction, internal convection, and pressure gradients, thereby removing spatial derivatives.
The spatial derivatives are replaced by their respective full domain integral resulting in the additional terms for the tank volume (V ) and the inlet area (A ) being introduced. For the semi-analytic model the mass equation becomes: where V is the system volume and A is the inlet area. Similarly, the energy balance suffers more changes as the thermal conduction terms are removed. The spatial derivatives are also replaced by the vessel volume and the inlet area by setting ∇ = 0 and integrating over V + A . Equation (24) becomes: The LDF model, from (7) is unchanged.

Optimisation
It is often desirable to reduce cycle time. This study employs the FGM approach to systematically distribute phase change material inside the adsorption storage medium in a manner that improves the thermal behaviour of such vessels, thus, improving their capacity for storing and delivering gas. Such a method implies a trade-off between the available region for adsorption and the amount of PCM to improve the vessel temperature. The adsorbent material replaced by phase change material must imply a temperature amplitude reduction that is enough to improve the adsorption capacity and compensate for the loss of adsorption area.

Optimisation of an ANG storage medium
The goal for optimally distributing phase change material inside of ANG storage media is to increase its capacity for storing and delivering gas. In this work, this goal is reached by maximising the amount of adsorbed gas inside the vessel. The semi-analytic model assumes a homogeneous solution, whereas the FGM optimisation method tackles the problem of generating a material distribution. The design variable, ε pcm , denotes the phase change material (PCM) distribution. The proposed problem considers a system in which phase change material and adsorbent composition can vary continuously in the range of ε pcm = 0 to 1. Here, the FGM optimisation methodology is employed to produce a graded distribution inside the domain.
The goal for the optimisation in this study is to maximise the mass of gas inside the vessel by the end of the adsorption cycle, which results in the optimisation problem written in (27): (27) where C 1 is: where r is the radial spatial coordinate and V m is the molar volume of gases in STP conditions.

Numerical implementation
This section describes the numerical implementation of the physics governing equations. First, the 2D model is explained in the first four subsections. Finally, the semianalytic model presented in the last subsection.

2D model
In this work, the finite element method (FEM) approximates the solution for pressure and temperature in the domain by solving the governing equations, and numerical integration is employed to solve the adsorption kinetics. Figure 7 shows the triangular element and the order of interpolation for each variable. Pressure, temperature, and adsorption ratio use first-order continuous Lagrange interpolation while ε pcm uses discrete Lagrange of 0 degree. This strategy is possible because adsorption kinetics does not have secondorder derivatives, and it is attractive because it provides a faster solution. FEM and sensitivities are implemented using Python programming language alongside with the FEniCS (Logg 2007) and Dolfin-Adjoint (Farrell et al. 2013) open-source software. Figure 8 summarises the whole process for obtaining the adjoint problem. Unified form language (UFL) handles the forward equations; then, Libadjoint creates the discrete adjoint equations; finally, FEniCS Form Compiler (FFC) package creates the adjoint code. Dolfin interface orchestrates all of these packages and tools.
The optimisation sequence follows the diagram presented in Fig. 9. Initially, a PCM distribution is defined. Then, the governing equations are solved, and Dolfin-adjoint transforms the physical problem solution into linear subproblems. Then, the adjoint problem is generated to make possible the evaluation of sensitivities. The optimiser chosen for this application is the augmented Lagrangian method from ROL (Rapid Optimisation Library) (Rapid Optimization Library 2020). The new solution is determined by quasi-Newton line search method. The second-order derivatives are obtained via the L-BFGS-B method. The optimiser receives the gradients and returns a new distribution of pseudo-porosity's (PCM distribution).

Weak formulation
In this section, the weak form of equations will be presented. The mass and energy balance equations weak form are described in Section 4.2.1. In Section 4.2.2, the numerical integration regarding adsorption kinetics is presented.

Variational equation-fluid and thermal problem
The weak form of the smooth function L(u) is defined by: Also, any Neumann boundary condition inside ∂ can be determined by: The Dirichlet boundary condition for the problem can be expressed as: where I ⊂ ∂ , also φ ∈ H 1 2 ( I ) and H 1 0 ( ) is the set of all functions in H 1 ( ) such that u = 0 on ∂ . From (4), and (5), the weak forms for the mass and energy balance can be written as: where ads, pcm, and both refer to the adsorbent process, phase change material process, and both processes.
Finally, the solution for the problem can be obtained by solving:

Numerical integration
This study treats the adsorption physics separately of the fluid and thermal problem, which variational equations solve. Sorption problem is solved by approximating the time derivative numerically using the backward Euler method. Therefore, first, differential equations will handle the mass and energy balance equation, thereby, updating the pressure and temperature state variables. After that, the numeric differentiation in time will tackle the adsorption problem and solve it for the new set of pressure and temperature, updating the adsorption portion for the next time step. After the numerical differential is applied, (7) becomes: where Q 1 is the new adsorption density, whereas Q 0 is the adsorption density of the last time step. The t is the time step.
The derivatives for P and T may be written as: where the subscripts 0 and 1 are the anterior state and the new state, respectively.

Helmholtz filter
The Helmholtz filter (Remacle et al. 2012) is applied to aid numerical convergence and reduce mesh dependency and checkerboard. The filter updates each design variable according to the local values, thereby acting as an averaging tool. The filter length determines the considered area during the update process. This tool is implemented by solving the following PDE (Remacle et al. 2012): where z is the filter length andε pcm is the function receiving the filter effect (smoothing). The implementation of the Helmholtz filter consists of the solution of a PDE to smooth the design variable. The weak form of the (38) :

Semi-analytic model
For the semi-analytic model, no PDEs are solved as only time derivatives are solved. Figure 10 shows the schematic of the semi-analytic model solving sequence. After material properties, pressure input and initial conditions are set, the mass balance equation calculates the velocity at which gas flows.

Results and discussion
This section presents the correlation and the results obtained by both models. In Section 5.1, the validation for both models is presented. The optimisation based on the adsorption 2D model is compared with the literature and results are presented in Section 5.1.1. In Section 5.1.2, the phase change material model is validated with the analytic solution of the Stefan problem. The semi-analytic model is compared with the 2D model in Section 5.1.3. In Section 5.2, the optimisation results are presented. In Section 5.2.1, the optimisation based on the semianalytic model is presented. In Section 5.2.2, the optimised material distributions are presented and compared with the homogeneous solution.

Model validation
This subsection presents the validation of the two models employed in this work. It will demonstrate the validated adsorption model, show the phase change material numeric model validation with a closed problem analytical solution, and then present the semi-analytic model validation.

Adsorption
In order to validate the adsorption model, it is compared to the experiment presented in Sahoo et al. (2011). The 30 L min case is used for the comparison. This scenario was simulated by assuming ε pcm = 0.0.
For a convection coefficient of h = 700 W m 2 K is assumed, and considering the thermal conductivity of copper (k = 385.0 W mK ), and the vessel external wall thickness as 2 mm, (w t = 0.002m), it is possible to calculate the Biot number (about 0.0036 for the cylinder shape). This Biot value Fig. 11 Bed geometry, finite element mesh, and cycle details indicates negligible thermal resistance for the tank external walls. Therefore, the role of tank wall thickness is neglected in the simulation.
For all charging cycle analysis, the initial conditions for the adsorbent domain are as follows: The results are compared with those found in the literature with a prescribed volumetric rate of 30 L/min (Sahoo et al. 2011). The vessel geometry, boundary conditions, dimensions, and the mesh employed in the simulation are shown in Fig. 11. The boundary condition for this specific analysis is set to prescribed volumetric flow at I , given by: where n is the unit outward on the surfaces.
The temperature in the inlet ( I ) is equal to the ambient temperature, which determines the convective heat transfer on the walls. Temperature boundary conditions are: Figure 12 shows the comparison between this work and (Sahoo et al. 2011). In Fig. 12a, the pressure evolution is presented. The current model shows a good correlation with the literature, presenting a difference of 8% pressure in the final steps. In Fig. 12b, the temperature at the centre of the bed shows good agreement with the experimental and numerical data found in the literature (Sahoo et al. 2011). Comparison reliability improves when the gas flow is verified. Figure 13 shows the total gas volume inside the vessel. The 30 L/min observed in literature (Sahoo et al. 2011) is also verified in the virtual analysis.  Figure 14 shows the temperature profile of the adsorption vessel at 3 min during the charging cycle. It matches the profiles found in the literature (Sahoo et al. 2011). The temperature profile, pressure, and temperature curves are in agreement with the data found in the literature.

Phase change model
The phase change material model is compared to the analytical solution of the second problem described in Stefan (1891), also known as the Stefan problem. The Stefan problem describes a situation in which a phase change phenomenon takes place in a homogeneous domain, and its analytical solution may be used to verify numerical models describing phase change (Ogoh and Groulx 2010). This scenario was simulated by assuming ε pcm = 1.0 and ε t (pcm) = 0.0. Figure 15 illustrates the boundary condition for the proposed problem and the domain. The wall temperature Tw is set to 350 K, whereas the domain temperature receives the value of 300 K. The phase change material properties are described in Appendix D. Phase change temperature is 313 K. The analysis simulates the transient heating phenomenon for 16 h. For all simulations, γ pcm = 1.8, and it is assumed that the phase change phenomenon takes place in a 3 K temperature interval. This scenario simulates the same conditions of the Stefan problem; which is a one-dimensional problem wherein a melting front moves from one side to another due to the melting process caused by a hot wall in one of the extremities.
An analytic model for the phase change phenomena is available in the literature in the form of a Laplace transform of the thermal energy balance equation. Appendix A describes the analytical equation for the Stefan problem. The Stefan problem explores a condition where a full PCM domain starts at the phase change temperature and with a heat source in one of its walls. The heat source allows the domain to undergo the phase change, and a temperature profile is created inside the domain. An analytic equation describes this temperature profile; it serves as a comparative to verify if the value chosen for γ pcm is capable of returning reliable results.   Figure 16 shows the temperature profile at a range of times. Results show a good correlation between data. Although the temperature inside the line is slightly colder than it should be in 3 h of analysis, the phase change wall position is correct. In later periods, the phase change barrier advances faster than it is expected. Nevertheless, for all the adsorption analysis considered in this study, a period shorter than 3 h will be considered; thereby making the designed parameters suitable for these purposes.
This validation is relevant because it adjusts the phase change material capacity. If unsuitable parameters are employed, the material will present altered phase-change enthalpy, then, behaving unrealistically.

Semi-analytic model verification
The semi-analytic model must be compared to the 2D model analysis to verify its reliability and quantify how much natural convection influences the results. All initial and boundary conditions are the same as those applied for the 2D model and are accordingly shown in Fig. 11, except for The thermal conductivity for the adsorbent is decreased due to the porosity; therefore, the only regions affected by the convection are those close enough to the border. Figure 14 depicts the temperature gradient presented in the operating cycle. Figure 17 shows that the semi-analytic model can predict the behaviour of the vessel in the experimented condition for a full carbon vessel well. The deviation in the temperature is caused by natural convection, which is present in the 2D model. When the amount of adsorbed gas is taken into account, it becomes clear that, for this condition, natural convection plays a minor role in the total amount of adsorbed gas. Therefore, the semi-analytic model is aligned not only with the 2D model but also with the literature. Figure 18 shows the same charts of Fig. 17 but for a 10% PCM scenario. It is possible to see the temperature plateau caused by phase change and the temperature deviation due to the convection. As expected, the 2D model has lower temperatures by the end of the analysis. Both models show good agreement in temperature and adsorbed gas.

ANG optimisation
A cylindrical adsorption vessel with adsorbent and PCM mixture in its interior is optimised in this section. The model is axisymmetric; a radial section is modelled, and it represents an adsorbent bed subject to heat exchange with the exterior through natural convection. As thermal Three conditions are analysed in this work: • Case1-Charging only: the vessel is exposed to a single charging cycle. The pressure rises from the initial state to the final state following a given function; • Case2-Charging followed by discharging: The vessel charges in the same fashion as in case 1, but there is a discharge sequence after that. The main characteristic that is present in this scenario is that the discharging process starts at the same pressure level that the charging process finishes and desorption also benefits from the heat generated by adsorption, which remains in the material internal energy or temperature; • Case3-Charging followed by discharging after thermal equilibrium is reached: This case is particularly curious because it may be associated with a transport situation. The vessel is filled like in case 1; however, after this process ends, it is assumed that it remains at ambient temperature until it reaches thermal equilibrium and, then, it performs the discharging cycle. In this situation, the gas that remained in the gaseous form inside the tank during the charging cycle would be adsorbed due to the temperature drop, thereby decreasing the pressure. For this case, the initial pressure for the discharging cycle is obtained by maintaining constant the total gas mass inside the vessel and setting the internal temperature to the same value as the external temperature. The value for the initial pressure of the discharging cycle is defined by solving (42): where T f and P f are, respectively, the final temperature of the charging cycle and the final pressure of the charging cycle, T amb is the ambient temperature, and P i is the initial pressure for the discharging cycle. Finally, the discharging cycle takes place at a different pressure Fig. 18 Comparison between the 2D model and the semi-analytic model built in Wolfram Mathematica for a 10% PCM mixture inside the tank than the end of the charging cycle, and desorption no longer benefits from the heat generated by adsorption.
The pressure input for each case is defined by a timedependant function, as presented in Fig. 19. The time period considered in the analysis was from 0 to 240. For case 2, as indicated in Fig. 19, the considered charging cycle starts in 0 s and ends in 300 s; after that, discharging occurs from that time until 3600 s. For case 3, the time period adopted for charging cycle was from 0 to 360 s and for the discharging cycle from 0 to 3600 s. For all cases, the number of time steps used in the simulation is 1200.

Semi-analytic model full factorial sweep
Solving the semi-analytic model has a low computational cost. Therefore, it is possible to perform a sweep in a range of variable values of interest to visualise the behaviour of merit function of (27). In this study, the swept variables are the phase change temperature for the PCM and its concentration. The sweep was performed for ambient temperatures of 290.0 K, 302.5 K, 315.0 K, and 327.5 K. This procedure was applied in all cases.
The total number of analysis for each case is 28,684. The PCM concentration was tested from 0.0 to 0.5 in steps of 0.005, whereas the phase change temperature was tested from 290 until 340 K in steps of 0.704 K. The total number of steps for each variable are, respectively, 101, 71, and 4.
If the phase change temperature and the ambient temperatures are fixed in 310.65 K and 300.0 K, just the PCM portion may be adjusted. Figure 20 shows how the adsorbed volume and the energy consumed in the phase change phenomenon behave. The optimised PCM portion is, in this case, unique and well defined. The adsorbed volume decreases with increasing PCM ratio beyond the optimised PCM ratio. The problem lies in finding the point where the energy consumed by PCM stops increasing linearly with the increase of PCM ratio. The energy rises, initially, with the increasing of the design variable. Its rising trend   Fig. 22 Case 1 adsorbed gas volume evolution during the filling cycle for full carbon (FC) and optimised PCM mixture for two external temperatures declines after the optimised point because the heat generated by the remaining adsorbent is not enough to saturate the PCM. Consequently, going beyond the optimised point would remove the adsorbent area without significantly increasing the energy stored in PCM, which causes the vessel performance to decline. In this case, the optimised PCM ratio is 14.0%, and the improvement is about 21% when compared to the full carbon vessel. Figure 21 presents the contour plot for the improvement in terms of adsorbed gas for all inlet temperature verified. The black regions indicate conditions where inserting PCM in the domain brings no benefit to the system. The improvement increases as the phase change temperature approach the ambient temperature. Because adsorption is the main phenomena, the phase change temperature must be higher than the initial temperature.
Case 1 reached a 23.2% improvement at 329.29 K of phase change temperature for an ambient temperature of 327.5 K. For all ambient temperatures, the improvement behaves similarly, and the optimised amount of PCM for this case is around 16%. Nevertheless, if more PCM is inserted in the domain beyond the optimised, the system loses adsorption capacity due to the trade-off between the adsorption area and temperature. The chart shows a small decrease in gain if phase change temperature decreases from optimal because PCM would start the analysis partially melted. Results suggest that the higher the ambient temperature is, the more beneficial the addition of PCM becomes. The full-carbon vessel can adsorb more gas at a lower ambient temperature than the same system at higher ambient temperatures. The change in the mass of adsorbed gas is shown in Fig. 22 for two different external temperatures, 327.5 K and 290.0 K and the full carbon vessel and a vessel with optimised PCM mixture. The internal temperature variation is presented in Fig. 23 while Table 2 shows the temperature variation during the filling cycle for two ambient temperatures, considering full carbon and optimised PCM mixture configurations. For both ambient temperatures, PCM causes a similar temperature amplitude reduction. Consequently, the total volume of adsorbed gas increase between the external ambient temperatures is similar, as presented in Table 3. However, as adsorption systems adsorb less gas at high ambient temperature, the increase is proportionally higher, which explains the larger gas increase with PCM with an increase of ambient temperature.  The summary of the optimised points found for each temperature are presented in Table 4: For cases 2 and 3, the improvement is measured by the total amount of gas the system can deliver by the end of the cycles. The total amount of delivered gas if the difference between the total amount of gas by the end of the charging cycle and the amount of gas by the end of the discharging cycle. Case 2 reached 21.7% of improvement at 331.43 K of phase change temperature for an ambient temperature of 327.5 K. The best PCM distribution is around 15.5%. The optimised amount of PCM increases as ambient temperature increases. If the ambient temperature is 327.5 K, the optimised amount of PCM is 15.5%. Figures 24 and 25 present the contour plot for the improvement in terms of adsorbed gas for all inlet temperature verified for these cases. As in case 1, PCM gains rise as external temperature rises. The temperature reduction during the filling cycle, due to PCM insertion, causes proportionally higher adsorbed gas increase at higher external temperatures; however, for case 2, the heat stored in PCM during the charging cycle aids desorption. As high ambient temperatures favour desorption, which is a slower phenomenon than adsorption, the differences between PCM gains for various ambient temperatures are significant.
The summary of the optimised points fount for each temperature is presented in Table 5. Figures 21 and 24 present a white region that address negative gains. As these gains are undesired, they are removed to identify the desired areas, which present improvements.
Case 3 reached 167.1% of improvement at 290.0K of phase change temperature for an ambient temperature of 290.0K. In this case, the best takes place for low ambient temperatures. The optimised amount of PCM for this case is around 23.5%.
Case 3 presents substantial gains due to the pressure drop in full carbon tanks. After the charging cycle, the vessel is at an elevated temperature; as the temperature drops, the gaseous gas inside the vessel when the temperature was high, is gradually adsorbed. Consequently, the pressure  inside the vessel drops. The result is a vessel with lower internal pressure and at a lower temperature. The lower pressure may be considered less gas that will be delivered without desorption, and the lower temperature implies lower rates of desorption. Therefore, the pressure drop is a double loss in terms of delivered gas, and it tends to lead to more residual gas by the end of the discharging cycle. When PCM is added in case 3, it decreases the pressure drop by reducing the maximum temperature reached during the charging cycle. Consequently, as the vessel reaches thermal equilibrium with the ambient, less gas in the gaseous form will be adsorbed after the charging cycle ended. The result is a vessel with an internal pressure at higher values and able to deliver more gas. Figure 26 shows the pressure drop between the full carbon vessel and the vessel with PCM for two ambient temperatures. At 290.0 K, the full carbon vessel presents a pressure drop of 46% while the 23.5% PCM vessel pressure drops by 8%. At 327.5 K, the full-carbon vessel presents a pressure drop of 42% while the 19.5% PCM vessel pressure drops by 18%. The effect of this pressure variation in the gas delivery capacity is presented in Fig. 27. After an initial increase, the curve stagnates during the filling cycle because the SA model does not consider convection, which means that the system cannot dissipate heat and increase adsorption equilibrium density. The increase of adsorbed gas after 360 s represents the gas that is in the gaseous state and is adsorbed after the charging cycle is finished. The full-carbon vessel has significantly more residual gas than the vessel with PCM by the end of the discharging cycle. The step in the curves after the charging cycle is caused by the pressure drop between the cycles.
Unlike the other cases, the best option for case 3 would be a PCM with phase change temperature equal to or below the ambient temperature. This contrast arises because desorption is the dominant phenomena. Additionally, the gains steadily increase as ambient temperature decreases. For case 3, the trend is the opposite of that observed in cases 1 and 2, where PCM gain increases when external temperature decreases. The explanation is that the dominant effect, in this case, is the pressure loss due to temperature variation between the charging cycle and discharging cycle. The pressure loss effect represents the adsorption of the gas in gaseous form by the end of the charging cycle but in adsorbed form by the start of the discharging process. Naturally, the full carbon vessel adsorbs less gas at higher ambient temperatures than at lower ambient temperatures. Consequently, the temperature variation at high ambient temperature is lower, resulting in a smaller pressure loss. Because PCM reduces the temperature increase during the charging cycle, it reduces the pressure loss, which is higher at lower ambient temperatures, thereby being more effective.
The summary of the optimised points found for each temperature are presented in Table 6: These results regard the homogeneous distribution of PCM inside an adsorption tank. In case 1, adsorption is the only phenomena taking place. Although case 3 starts with adsorption and ends with desorption, being the latter, the limiting phenomena. This dominance determines if the phase change temperature is higher or lower than the ambient temperature. However, for these cases, varying the ambient temperature does not modify the optimised amount of PCM. In case 2, adsorption and desorption are both taking place, and desorption is influenced by adsorption for the discharging cycle starts in the state that the charging cycle ends. In this case, the favoured phenomenon is adsorption. The reason for that might be related to storing as much energy as possible in the PCM internal energy so it can be used to aid the desorption phenomenon.

Graded material optimisation
In this section, the benefits brought by optimised homogeneous PCM distribution inside an adsorption vessel are quantified. As the semi-analytic model does not consider the effects caused by the boundaries convection boundary conditions, it is expected that it is possible to go beyond the homogeneous optimised solution when applying FGM optimisation to distribute material inside the vessel properly. It is expected that the larger the container, the closer the 2D model and the SA model match, while for smaller tanks, the gain brought by the graded solution tends to be higher. FGM implementation considers the parameters shown in Fig. 11. The initial PCM condition is a homogeneous distribution with all design variables set according to the homogenised optimised PCM distribution obtained in the previous subsection. For the optimisation of the 2D model, only case 1 will be run. It will consider an ambient temperature of 300 K, the initial PCM distribution is set to ε = 0.0% for all domain, and the Helmholtz filter length is set to z = 0.01. The ROL (Rapid Optimisation Library) (Rapid Optimization Library 2020) optimiser applies augmented Lagrangian method guided by quasi-Newton line search method (Rapid Optimization Library 2020) with L-BFGS-B.
The optimiser chosen for this application is the augmented Lagrangian method from ROL (Rapid Optimisation Library) (Rapid Optimization Library 2020). The new solution is determined by quasi-Newton line search method. The radial section of the cylindrical vessel considered for optimisation is the same as that used for adsorption validation as depicted in Fig. 11. The 2D FGM optimisation considers the scenario of natural convection and for a convection coefficient of 700 W/m.K, assumed at boundaries, representing a value describing efficient surface cooling.
It is verified that for lower volume rates during the filling cycle, the tanks present the highest efficiency, while for higher volume rates, the efficiency is lower (Sahoo et al. 2010). This happens because, for lower volume rates, the cycle time is higher than for lower volume rates, which allows convection to remove more heat from the system during the filling cycle, resulting in a cooler system. This picture may be addressed as a time-scale problem. For the SA model, which does not consider convection, the timescale does not affect the result. For the 2D model, this time-scale problem is handled by altering the heat transfer rate, instead of filling time so that convection will take more or less heat from the system.  Figure 28 shows the results obtained for the vessel optimisation considering 300.0 K as ambient temperature, a low convection rate, and a cycle time of 240 s. It also depicts the comparison between the temperature and amount of gas for three-vessel configurations, being them, full carbon, homogenised solution (of 14% PCM), and optimised solution. As expected, the distribution is almost homogeneous throughout most of the domain, placing no PCM where the adsorbent has its temperature reduced by the boundary effect. PCM concentration ranges from zero to 14% for the slower filling cycle. When compared to the homogeneous solution, the amount of gas that the graded vessel is capable of storing is about 1% more. This improvement comes from excluding PCM from the boundary regions, which are already cooled by the air intake and convection. The integral of the total amount of PCM inside the vessel reduces from the homogeneous solution (14.0%) to 13.7%, which is expected as the optimiser removes PCM from the regions close to the borders. The variation between the homogeneous and the graded solution is low because the external convection does not remove enough heat to cause a significant change in results.
In Fig. 29 considering 300.0 K as ambient temperature, a high convection rate. As in the low convection case, the border effects such as external convection and inlet gas temperature affects the internal PCM distribution. The optimised PCM distribution in his case stores 2% more gas than the homogenised solution. The low variation between the homogeneous and FGM solution in both cases is related to the low thermal conductivity of the adsorbent and also the adsorption isotherms, which affects how adsorption capacity changes according to the temperature. Figure 30 shows the evolution of the total PCM ratio throughout the optimisation. The vessel with low convection converges to 13.7% while the tank with high convection converges to 8.4%. This variation is caused to the external convection boundary condition. Table 7 summarises the total adsorbed gas volume by the end of the analysis for the three-vessel configurations. The FGM solution adsorbs more gas than the other, especially when boundary condition effects are more pronounced. In Table 8, the improvements of each solution are compared. At low convection conditions, PCM brings more improvement to the system. In this case, the variation between the homogeneous solution and the FGM solution is around 1%. For the high convection scenario, the difference between the solutions doubles, reaching 2%. The variation between both solutions is low in both cases. The FGM solution is more suitable for cases with higher boundary effects magnitude. However, the same table shows that at high external convection values, the benefits brought by PCM are lower than those for low convection. The computational time of both models is presented in Table 9. The 2D model spent 24,736 s during its computation while the SA model spent 1304 s (about 5% of the 2D model). Both models were run on an Intel(R) Xeon(R) CPU E5-2630 v2 @ 2.6 GHz with 24 cores and 70 Gb of RAM.
A virtual analysis for the charging cycle was run with the high convection optimised PCM distribution to compare with the full carbon configuration. Figure 31 presents the evolution of the total amount of gas inside the system during the charging cycle. As expected, the optimised vessel starts the analysis with less gas because there is less adsorbent inside it but finished the process with more gas, because the mixture has a higher adsorbed gas density due to temperature reduction. The temperature and adsorbed gas mass are already presented in Fig. 29. These results are sensitive to the adsorbent material. The materials here used were chosen because their properties are easily found in the literature. The material adsorption characteristics and the adsorption heat play substantial roles in this optimisation; likewise, the phase change temperature and the phase change enthalpy also does.
It is important to remark that the axisymmetrical nature of this analysis hinders optimisation. It makes it impossible for the optimiser to build 3D patterns. A full 3D analysis

Conclusions
Two approaches were presented to distribute PCM material inside an adsorption system optimally. The semi-analytic model describes a homogeneous PCM distribution that is suitable for vessels whose dimensions are large compared to the inlet and external convection effects because it does not consider local effects. The 2D axisymmetric optimisation was compared to the semi-analytic model to verify how much the boundary effects could represent a change in the vessel optimised PCM distribution. Results show that for cases where convection is not negligible, the semi-analytic model has its benefits reduced. The 2D model is indicated for adsorption systems with longer cycle times and/or higher convection rates due to the effect that the boundaries have on the solution. The 2D axisymmetric results use less PCM than the homogeneous result from the semi-analytic model. Consequently, the 2D model may be indicated when the volume of used PCM is a concern.
The semi-analytic model presented a good agreement to the results found in literature when convection is low. The full factorial sweep returned results which provide improvements that are very close to those obtained via FGM optimisation in these cases. Therefore, this approach presents a very attractive cost-benefit when optimising an adsorbent system with mixtures. As the adsorption phenomenon is time transient and complex cycles may demand many time steps, its computational cost may be a  problem for the 2D model, which does not happen in the semi-analytic approach. For all cases, the optimised result is an adsorbent and PCM mixture. The mixture removes part of the adsorbent and increases the adsorbed density of those adsorbent particles that remain. The result is that the vessel stores less gas at standard conditions for temperature and pressure, which is desirable because it means less residual gas by the end of the discharging process, and more gas by the end of the charging process because the temperature amplitude inside the vessel during the cycle is reduced by PCM. For case 3 in SA model, PCM reduced the pressure loss after the charging cycle.
In this study, adsorption simulation was carried out as a transient phenomenon. However, if it could be addressed as a stationary problem, the efficiency and computational cost for the analysis would be reduced, and more complex cycles could be optimised via FG or topology optimisation. The stationary approach is available for further exploration.
Another aspect of evaluating may be whether it is convenient to develop technology to build non-homogeneous adsorbent systems. This technology could employ additive manufacturing or even create several internal chambers inside the vessel, each filled with its respective adsorbent and PCM mixture. Nevertheless, results present small variations between the FG solution and the optimised homogeneous solution. Consequently, the cost and effort to develop the technology for building these vessels are unlikely to be justified.

Replication of results
The implementation in the FEniCS platform is direct from the description provided of the equations and numerical implementation in this article. FEniCS is based on a highlevel description for the variational formulation (UFL); therefore, the generation of the matrix equations is automated.
The pseudocode of the implementation is indicated in Appendix B and C, where the main FEniCS/dolfin-adjoint functions that are being used are represented between Fig. 31 Comparison between full carbon vessel and the optimised design in virtual analysis Ɵ parentheses. When using dolfin-adjoint, the dolfin-adjoint library provides an interface to ROL. In Appendix B, the pseudocode for the semi-analytic model describes the procedure for each case.
Finally, the solution for the temperature along the system is given by:    (Limousin et al. 2007) source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.