Topology optimisation of biphasic adsorbent beds for gas storage

Adsorption is a retention mechanism of fluid molecules on solid surfaces and presents a wide range of applications, such as fuel storage, refrigeration and separation processes. This work describes the modelling of gas adsorption on porous media and presents an optimisation approach for the design of adsorption systems based on biphasic adsorbent beds by topology optimisation. A comprehensive formulation for the adsorption phenomenon is presented, detailing the derivation of governing equations and respective weak forms and discretisation for the implementation of the finite element method (FEM). A new topology optimisation material model based on offset hyperbolic tangents is introduced. The derivation of sensitivities is presented in detail, based on a transient adjoint problem. A diverse set of optimised adsorbed natural gas (ANG) tanks, considering real material properties of activated carbon and steel, is presented. Results indicate the suitability of the method in optimising the distribution of phases across adsorbent beds and show that biphasic ANG tanks can perform significantly better than traditional tanks.

of solid particles attracted by disrupted bonds (Mota et al. 2004), while the latter is a bulk phenomenon, in which gas molecules diffuse through another phase, accumulating at the grain boundaries and often forming a new phase (Blanco et al. 2014).
Adsorption is divided into two categories: (i) physical adsorption, based on weak van der Waals forces, which is also know as physisorption and (ii) chemical adsorption, based on covalent bondings, also know as chemisorption. While the former is easily reversible and leaves the adsorbent intact, the latter involves a chemical reaction that affect the surfaces, as is the case of corrosion phenomenon. Figure 1 illustrates the mechanism of adsorption, where molecules concentrate on micro-pores walls.
The capacity of adsorption of porous media is dependent on characteristics such as the distribution and characteristic of micropores on the adsorbent and the affinity of adsorbent and adsorbate (Biloė et al. 2002). The capacity of adsorption is also directly proportional to the pressure and inversely proportional to the temperature. During an adsorption cycle, the adsorbent bed temperature naturally tends to increase, since the adsorption reaction is exothermic, releasing heat when a molecule of gas is adsorbed by the surface. Conversely, the desorption reaction is endothermic, consuming heat when molecules detaches from the surface. Therefore, the major challenge in designing adsorption systems is often the internal temperature variation, which prevents the system to reach its full capacity of gas adsorption and hinders the adsorbent bed purge.
Several engineering applications are based on the adsorption process, such as fluid separation, refrigeration, CO 2 capture and storage (CCS) and gas fuel storage. In the case of natural gas storage, adsorbed natural gas (ANG) tanks consist of a solid matrix of activated carbon (AC) with natural gas molecules accumulated in its pores. Due to the relatively high density of the adsorbed phase, an overall storage capacity around 164 V/V (164 times denser than at atmospheric temperature and pressure) can be achieved at room temperature with pressure ranging from 35 to 50 atm (Hirata et al. 2009). Conversely, compressed natural gas (CNG) is ordinarily stored at 200 atm, offering a capacity of 245 V/V. In turn, liquefied natural gas (LNG) requires cryogenic temperatures to force gas condensation, which happens, for methane, around 113 K and provides a storage capacity of 600 V/V (Hirata et al. 2009). These characteristics place ANG as a promising method to maintain gas confinement, as it works at atmospheric temperature and much lower pressure than CNG.
A common objective when designing ANG tanks is to maximise the mass of gas that can be stored and recovered in reasonable short adsorption and desorption cycles. Such tanks require materials presenting high capacity of adsorption and appropriate thermal management. The amount of gas a bed is able to adsorb is determined by its microstructure. Activated carbons presenting the highest adsorption of natural gas typically have narrow microporosity and micropores width of 2.0 nm (Biloė et al. 2002). Monte Carlo simulations determined that the theoretical maximum capacity of adsorption for methane is 209 V/V on monolithic carbon and 146 V/V on pelletised carbon (Matranga et al. 1992). It has been shown that, due to the temperature variation, real discharge cycles in ANG tanks may present between 8 and 29% loss of capacity when compared to ideal isothermal cycles (Chang and Talu 1996). Thermal management can be addressed by varying the tank shape and boundary conditions. Different flow rates, aspect ratios and convection coefficients can increase significantly the storage capacity of cylindrical tanks when imposing forced convection on relatively slim tanks (Sahoo et al. 2011). Adsorption of hydrogen on activated carbon has been improved about 25% by inserting equally distributed metallic fins on an adsorbent bed connected to a cold source at 233 K (Vasiliev et al. 2014). Active cooling has also been used to improve ANG tanks performance, through the addition of cooling jackets (Mota et al. 2004) and fin-tube type heat exchangers (Rahman et al. 2011).
The use of topology optimisation on the design of adsorption systems has not been reported in literature, although key concepts can be inherited from heat transfer and fluid problems. Topology optimisation has been applied to heat transfer problems for many years. Cases include the cross-section optimisation of a minimum-material fin (Anagnostou et al. 1992), material optimisation in transient heat conduction problems (Turteltaub 2001), electrothermal microactuators subjected to top heat convection (Sigmund 2001) and convection-dominated problems (Bruns 2007). The application to fluid systems is introduced by the design of channels seeking the minimisation of power dissipation (Borrvall and Petersson 2003), which has been extended to incompressible laminar viscous flow at moderate Reynolds numbers (Gersborg-Hansen et al. 2005) and generalised for scaled for Stokes and Darcy equations (Guest and Prėvost 2006). The use of topology optimisation in problems coupling heat transfer and fluids is also reported, as in the optimisation of heat dissipating structures with forced convection (Yoon 2010;Koga et al. 2013;Matsumori et al. 2013;Kontoleontos et al. 2012) and in the design of optimised electro-thermo-mechanical (ETM) actuators by considering their interaction with the surrounding fluid (Yoon 2012). A detailed review of topology optimisation applied to heat transfer and fluid flow systems is presented in Dbouk (2017).
The objective of this work is to design optimised ANG tanks by distributing non-adsorbent materials across the adsorbent bed, as demonstrated in the tank depicted in Fig. 2, to achieve high storage densities faster than full-carbon tanks. Given the complexity of adsorption systems, such design problem is approached by topology optimisation. The objective function is the maximisation of gas intake during a filling cycle of fixed length. The transition between adsorbent and non-adsorbent materials is made through an offset hyperbolic tangent material model. Resulting topologies are based on the optimisation of twodimensional domains, representing axisymmetric or very long tanks. This paper is organised as follows. In Section 2, the theoretical modelling of adsorption systems is presented. In Section 3, the optimisation problem is defined with respective material models and objective functions. In Section 4, the numerical implementation of the adsorption model and sensitivities are presented. Optimised topologies are presented and discussed in Section 5. Concluding remarks are presented in Section 6.

Modelling of adsorption systems
A typical adsorption domain Ω is illustrated in Fig. 3, where Γ W and Γ L denote sub-regions from ∂Ω where boundary conditions of wall and gas inlet/outlet are imposed, respectively.
In this work, the following assumptions are made: Ideal gases: in the systems considered, gases are considered to present ideal thermodynamic behaviour (Sahoo et al. 2011), therefore, their densities are described by: where M g is the molar gas mass, R g is the universal gas constant, P is the gas pressure and T is the gas temperature.
Darcy's Porous media: the Darcy's law governs the fluid flow through porous media, making the velocity to be given by:

Fig. 3 Adsorption domain
Local thermal equilibrium: the adsorbent bed is locally in thermal equilibrium with the gas (Mota et al. 2004). All intra-particle and film resistances in heat and mass transfer are neglected.
Constant specific heat capacity: the specific heat capacity is considered to be constant within the typical range of temperatures observed in adsorption systems. Its variation is neglected since it usually varies around only 8% (Mota et al. 2004) in a range of 100K.
Constant adsorption enthalpy: this work considers an approximation in which ΔH assumes a constant mean value within the range of pressures and temperatures studied (Sahoo et al. 2011).

Lagergren adsorption model:
For the studied cases, it is considered that the film diffusion is dominant and the rate of adsorption is governed by Lagergren model (Vasiliev et al. 2014).

Dependency of material properties on porosities
Effective properties of porous media are dependent on macro-porosity (ε M ) and microporosity (ε m ). Given these two classes of porosity, the total porosity (ε t ) is defined as: Density: it is given by the sum of gaseous and solid phase weighted by the total porosity, ε t , with the addition of the density of adsorbed gas, Q, which is considered to occupy a neglectful volume: Permeability: dependent on the macro-porosity, following the Kozeny-Carman's model for spherical particles (Carman 1937): where ø s is the effective characteristic particle diameter and τ s denotes the pores tortuosity.
Thermal conductivity: it can be given by both Maxwell-Eucken's or Russel's relations (Smith et al. 2013). The former considers the effect of large spherical pores with thin walls and the latter the effect of a cubic pore inclusion in the material matrix by considering a parallel/series resistor.
Here, the conductivity of the thin film of gas adsorbed, k q , is neglected and the Maxwell-Eucken relation is adopted to interpolate the gaseous and solid conductivities (Smith et al. 2013): where k s and k g are the thermal conductivities for the solid and gas, respectively.

Adsorption and desorption kinetics
The concentration q represents the mass of adsorbed phase per mass of adsorbent. Therefore, the density of adsorbed phase in reference to the volume of the adsorbent bed, Q, is given by: The specific volume of micropores, w s , represents the volume of micropores per mass of adsorbent. The ratio of micropores volume per adsorbent bed volume, W s , is then given by: From the equilibrium concentration of adsorbed gas, q eq (Sahoo et al. 2011), the equilibrium density of adsorbed gas is given by: where β is the affinity coefficient of the adsorbentadsorbate pair, E 0 is the characteristic energy of adsorption and n s stands for the effective micro-pores dispersion. The local adsorbed gas density, ρ adg , is given by Sahoo et al. (2011): where ρ boil is the gas density at boiling point temperature (T boil ) and α e is the mean value of the thermal expansion of liquefied gases. The Polanyi's adsorption potential, A, is given by: which saturated pressure P sat is given by: where P cr and T cr are, respectively, the critical pressure and temperature of the gas. The rate of adsorption decreases as more gas is adsorbed, following the LDF model given by the following expression (Mota et al. 1996): where the coefficient G, which dictates the adsorption and desorption kinetics, is given by Vasiliev et al. (2012): where the activation energy of adsorption (E a ) is considered when Q ≤ Q eq and the activation energy of desorption (E d ) otherwise. D is the intra-particle diffusional time constant, given by Sahoo et al. (2011): where ø s is the characteristic particle diameter and D 0 is the gas diffusivity.

Continuity
In adsorption systems, the coupled continuity equations of gaseous, solid and adsorbed phases is expressed by Sahoo et al. (2011): where Q is the mass of adsorbed gas per unit of adsorbent volume and ρ s is the bulk adsorbent density, which is constant. By substituting the ideal gas and Darcy's laws, the balance becomes: where the solid density ρ s is considered constant. A pressure inlet condition is imposed on Γ L , representing the connection of the tank to an infinite reservoir at fixed pressure and temperature or the opening of the tank to the atmosphere. Inlet pressure is varied smoothly, being approximated by the following sinusoidal profile: where P ini and P let denote initial and filling/emptying pressures, t f is the final time and σ is a smoothness factor. Figure 4 illustrates P Γ L curves for the adsorption cycle.

Energy balance
The energy balance of the gaseous phase is given in terms of Cp g and T as Bird et al. (2007): where the specific heat capacity, Cp g , is considered constant, τ denotes the viscous stress tensor, ∂lnV ∂lnT = 1 for ideal gases and D denotes the Lagrangian derivative.
Following the assumptions of negligible viscous dissipation and Darcy's porous media flow, the energy balance becomes: where the continuity equation is observed in C, making these terms to vanish, and the last term is neglected, given the very small pressure gradients observed in the range of permeabilities of interest in this work. Analogously, the energy balance of the solid and adsorbed phases are given by: where ΔH denotes the isosteric heat of adsorption. The energy balance for coupled gaseous, solid and adsorbed phases assuming local thermal equilibrium between phases (Mota et al. 2004) and ideal gas is then given by the following expression: where the effective thermal conductivity, k f , is given by (6). Inlet heat flux is considered on Γ L , which is expressed by the following Robin boundary condition (Sahoo et al. 2011): Assuming that heat convection takes place on the walls, a Robin boundary condition is defined on Γ W : On the region Γ \ (Γ W ∪ Γ L ), since no boundary terms are explicitly defined, the heat flux is null, representing adiabatic condition. Figure 5 illustrates the typical variation of the average gas density, temperature and pressure in an ANG tank used as an on-demand source of fuel, such as in gaspowered vehicles or in residences. Starting from the average density of total gas at the equilibrium, (Q + ρ g ) ini,a , the tank admits gas while pressure and temperature increase during the filling cycle, reaching the final average total density of (Q + ρ g ) f in,a . Such average density is conserved inside the tank until its future depletion. While the tank is kept quiescent, its temperature slowly approaches the equilibrium with the atmospheric temperature, T atm , causing its internal pressure to drop due to the gas phase change from gaseous to adsorbed. The stored gas is then delivered through a nearly isothermal desorption, maintaining equilibrium with T atm during the whole emptying cycle.

Topology optimisation of ANG tanks
The goal when designing ANG tanks applied to ondemand consumption is to maximise the total intake of gas during the filling cycle by mitigating the rise of the adsorbent bed temperature. A possible approach to improve heat transfer is not to occupy the full volume of the tank  However, the temperature reduction must be enough to compensate the lost of adsorbent material volume. Figure 6 illustrates such trade-off for a tank presenting the final average temperature of 323 K at the end of the adsorption cycle. It is shown how much volume is possible to occupy with non-adsorbent materials according to the reduction of the temperature in order to achieve certain gain in the mass of gas adsorbed. For instance, if the final average temperature of tank is 293 K, the same mass of gas could be adsorbed using only 67% of the original adsorbent volume. By using 80% of the original volume, it would be possible to increase the total amount of gas adsorbed in 20%. Topology optimisation is adopted to find optimised distributions of non-adsorbent material across the adsorbent bed. The alternation between adsorbent (ads) and nonadsorbent (non) materials in the domain is expressed by the pseudo-porosity e as follows: where the actual property Π(e) is given by Π ads or Π non depending on e value.
However, such binary law of material defines an illposed problem, with multiple local minima (Bendsøe and Sigmund 1999). Such issue can be avoided by relaxing the definition and assuming intermediate values, ranging from 0 to 1. Since only 0 and 1 represent existent materials, it is desirable to recover the discrete nature of this material law.
In the optimisation of biphasic adsorbent beds, intermediate materials often constitute the best solution from the optimiser perspective. For instance, in the case of seeking an optimised placement of steel on a carbon bed, the conductivity of the non-adsorbent material is usually more than 200 times larger than adsorbent's and steel adsorption capacity is null. Therefore, a material presenting half of steel conductivity and half of carbon adsorption capacity is a good solution, although it is not a real material. 1.23 kg Improvement: 6.03%   (27) is suitable for this sort of necessity. For the property Π A , where the functional of interest J (Π A,ads ) (adsorbent property in the whole domain) is higher than J (Π A,non ) (non-adsorbent property in the whole domain), the phase change is centred on e c = 0.5 + . Conversely, for a property Π B which presents J (Π B,non ) higher than J (Π B,ads ), e c = 0.5 − is adopted. In the absence of an offset, there would be intermediate properties at e = 0.5 regardless of the penalisation value. Figure 7 illustrates such material model for the offset = 0.02. where: with sgn denoting the signum function.
Topology optimisation presents inherent issues such as checker-board instability and mesh dependency, which are mitigated in this work by filtering the design variables before inputting them to the material model based on the solution of the following Helmholtz-type PDE (Lazarov and Sigmund 2011): (29) where the distribution is smoothed to the C 2 function e, according to the filter length parameter z.
The final optimisation problem for the design of ANG tanks for on-demand consumption is expressed as:

Numerical implementation
Modelling and optimisation of adsorption systems are implemented in this work aided by FEniCS (Logg 2007), its module Dolfin-Adjoint (Funke and Farrell 2015) and Scipy's L-BFGS-B optimiser (Byrd et al. 1995;Jones et al. 2011). The numerical formulation of such implementation is presented in this section.

Adsorption problem
The solution of the adsorption governing equations is approximated via finite element method (FEM). In order to do so, the inner product of the strong PDEs and arbitrary test    which is rewritten as: where the state variable w is constituted by pressure, temperature and adsorbed gas: Based on implicit Euler method for time discretisation, the n-th timestep of such system is given by: The following lower-triangular system is then defined: where w ia denotes the initial state conditions and the matrix C is given by: where:

Sensitivities calculation
The sensitivity of the functional J a in respect to the design variable e j is obtained via adjoint method, which is described in detail in Appendix A.3. The final sensitivity expression is given by: where δ denotes the Dirac's delta (centred on time instants t f a and t ia ), λ is the adjoint variable and B is the total gas mass, which is defined in Appendix A.3 by (75). Despite the functional and PDE systems being expressed in terms of filtered pseudo-densities, e j , the optimisation is based on non-filtered pseudo-densities, .
The discrete Helmholtz filtering system is given by: where Z and are defined in Appendix A.2.
The implicit derivative of such system in respect to results in: (42) where the first term is null, given that Z is not dependent on pseudo-porosities, and the right-hand side is: in which δ jk denotes the Kronecker's delta column vector, being null in all lines, except for the j -th, which is equal to 1.
The sensitivity of the functional J to the non-filtered pseudo-porosity is then given by:

Optimisation algorithm
The optimisation algorithm follows the flow represented in Fig. 8. In the optimisation loop, the sensitivities are obtained  via dolfin-adjoint based on the equilibrium equations solved on FEniCS. The optimiser L-BFGS-B (Byrd et al. 1995;Jones et al. 2011) receives the sensitivities and output a new distribution of design variables, aiming to improve the functional of interest. The relative functional variation is defined by considering the last two iterations: When the functional relative variation is non-negative and smaller than a tolerance ω min or when the maximum number of iterations it max is reached, the loop is broken. The penalisation p is then increased by Δp, following

Results and discussion
This section presents examples of the method implemented on a set of different tanks constituted by activated carbon subjected to methane adsorption.

Prismatic tanks
A two-dimensional reduction of a three-dimensional prismatic tank is presented in Fig. 9. It consists on the approximation of a very long tank with fixed cross-section. Based on symmetry, only one quarter of such section is considered. In this domain, the impact of the gas inlet is not present and the pressure is constant on the whole section, making the convective term in the energy balance null. In this and following examples, the adsorbent bed is considered to be wrapped by a thin sheet of metallic material subject to heat exchange with the room through convection. The Table 5 Optimisation parameters for cylindrical tank contribution of this wall to the energy balance is assumed to be neglectful (Sahoo et al. 2011) and it is not included in the model. In fact, being a thin sheet, its Biot number (Bi) is very small and the heat conduction inside it is much faster than the heat convection away from its surface. TOM implementation is based on the parameters shown in Table 1. It consists of a square section 0.5 m wide subjected to forced convection of 100 W/m.K. The adsorption cycle lasts 10 800 s (3 h), in which the pressure is smoothly increased during the first 1 800 s from 0.1 to 1.0 MPa. The material model defined in (28) with an offset = 0.02. The optimisation procedure is based on L-BFGS-B optimiser and starts with the homogeneous distribution of the design variables values at e = 0.5. The penalisation is smoothly increasing from p ini = 1 to p f in = 51, following the continuation method limited to a maximum of 20 iterations per continuation step (Table 1).
The resulting topology depicted in Fig. 10 shows perfectly symmetric steel fins connecting the central region of the section to the walls subjected to forced convection. The threshold defining the transition from activated carbon to steel is set on e = 0.5. As seen in Fig. 11, the material model adopted makes such interface clearly defined. The convergence of the objective function is shown in Fig. 12. Figure 13 shows that the fins made the maximum temperature inside the domain drop from 325 to 310 K, allowing 6.03% more gas to be admitted during the filling cycle (Table 2), despite the adsorbent volume lost (Fig. 14).

Complex axisymmetric tanks
An axisymmetric domain based on the radial section shown in Fig. 15 is considered to demonstrate the implementation of the method to complex shaped tanks. Walls on Γ W are subjected to natural convection (h = 5 W/m.K) while a fixed temperature of 273 K is imposed to the region denoted by Γ F . The full description of the model parameters is presented in Table 3.
The resulting optimised topology is seen in Fig. 16. Metallic limbs occur mainly in the vicinity of the cold wall Γ F , but also around the inlet region Γ L . Such distribution   Fig. 18. The convergence of the objective function is shown in Fig. 19.
The strategy adopted by the optimisation algorithm is to grow metallic limbs across regions with high temperature gradient. Figure 20 shows a reduction of the central zone temperature, which increases the density of gas as observed in Fig. 21. Such strategy leads to an improvement of 9.86 % in the total intake of gas during the filling cycle, as seen in Table 4.

Cylindrical tanks with high heat capacity materials
The impact of the increase of the heat capacity of the non-adsorbent material is studied here by considering the cylindrical domain shown in Fig. 22. Materials presenting high heat capacity can be used to roughly approximate phase-change materials (PCM) (Toledo et al. 2013). This example is based on the parameters listed in Table 5, considering natural heat convection and a non-adsorbent material presenting specific heat capacity four times higher than Steel (Cp non = 1880 J/kg.K). A different strategy of penalisation continuation is adopted here, starting at p ini = 20 and increasing the penaliser by Δ p = 10 every 10 optimiser iterations. Figure 23 shows the resulting topology for this case, where the non-adsorbent regions are not necessarily connected to the gas inlet or to the walls, given that materials with high heat capacity act as heat sinks. The distribution of the optimised filtered pseudo-porosity is shown in Fig. 24.
The convergence of the functional is shown in Fig. 25. It can be seen a different behaviour when compared to previous cases (Figs. 12 and 19), where the functional converges to lower values as the penalisation is increased. In previous cases, the low initial penalisation (p = 1) allows the occurrence of materials presenting intermediate properties, which would be a good solution if they were feasible. The main concept of the material model based on offset hyperbolic tangents is here verified, as it makes intermediate materials less interesting as the penalisation increases. Given that the optimisation of the cylindrical case starts from a moderate penalisation (p = 20), optimised layouts exploring artificial intermediate materials are not interesting from the beginning and the functional value increases along the iterations.
Non-adsorbent islands posses a radius of influence around them, resulting in an overall decrease of the adsorbent bed temperature, as seen in Fig. 26. Thereafter, the gas density field depicted in Fig. 27 becomes much more flat. Comparing with the homogeneous case, the mass of  gas admitted during the filling cycle is 10.00% higher when compared with the domain filled with activated carbon only, as shown in Table 6.

Remarks on the initial guess and penalisation
Topology optimisation problems are non-convex and often lead to different solutions depending on the initial guess and Fig. 28 Pseudo-porosity fields evolution for homogeneous (left) and arbitrary (right) initial guesses for prismatic tank optimisation starting with p ini = 1 on the strategy adopted for the penalisation continuation on the material model. The optimisation of the prismatic tank presented in Section 5.1 has been shown independent of the initial guess, as seen in Fig. 28, where the evolution of the optimisation starting from an homogeneous pseudo-porosity distribution e = 0.5 is shown on the left column while one starting from an arbitrary initial guess is shown on the right. In both cases, the distribution of pseudo-porosities goes towards a full-carbon domain during the first optimisation iterations, following an identical evolution after that towards the same final layout.
The cylindrical tank optimisation presented in Section 5.3 is ran from a material model which is already significantly penalised (p ini = 20). In this case, different initial guesses result in different final layouts as seen in Fig. 29, where the distributions on the top of each block show the evolution of the homogeneous initial guess, while the bottom ones refers to an arbitrary initial distribution of pseudo-porosity. Totally different layouts are obtained for each case, although the final objective function values are very similar: 2.31 kg of total gas intake for the homogeneous initial guess and 2.28 kg for the arbitrary.
The dependence of the final layouts on the initial guesses can be mitigated by following a continuation approach starting from a material model presenting low penalisation. Figure 30 shows the evolution of the pseudo-porosity distribution for an optimisation starting with penalisation p ini = 1 . Similar to the prismatic case, shown in Fig. 28, the initial optimisation steps take the pseudo-porosities to a full-carbon distribution, following the same evolution from then on. The final layout for this case results in a total gas intake of 2.30 kg, which is very similar to the 2.31 kg observed for the case starting the penalisation from p ini = 20. This particular case also highlights the nonconvexity of such optimisation problem, where a particular distribution of pseudo-porosities in a moderately penalised material model complies the objective function better than a less penalised one following a smooth continuation.

Conclusions
In this work, a topology optimisation method for the design of adsorption systems is formulated and tested for ANG tanks. Accordingly, a comprehensive model for ANG tanks is presented, including the kinetics of adsorption and the variation of adsorbent properties in function to the porosities. The solution of the resulting set of PDEs governing the adsorption problem is obtained via FEM. The procedure to assemble the FEM matrices from the PDEs weak form is described in detail. The implementation of this method is based on the calculation of the sensitivities via the adjoint method for transient systems, which is comprehensively presented.
Examples of ANG tanks of different shapes subjected to different boundary conditions are studied. Tanks with optimised distribution of non-adsorbent materials presented significant gains in the total mass of gas admitted during the filling cycle in comparison to tanks constituted solely by an adsorbent material. The improvements in total gas mass intake observed are 6.03% in a prismatic tank subjected to forced convection, 9.86% in a complex axisymmetric tank subjected to natural convection and a cold wall, and 10.00% in a cylindrical tank subjected to natural convection with an optimised distribution of a high heat capacity phase.
Topology optimisation has been proven to be a suitable approach for optimising the distribution of different phases across adsorbent beds to improve gas storage. In particular, good results were achieved for biphasic adsorbent beds constituted by materials presenting high heat capacity, encouraging further research on the use of phase-change materials (PCM) in ANG tanks.
Funding information This paper received the support of the RCGI Research Centre for Gas Innovation, sponsored by FAPESP (2014/50279-4) and Shell. Authors also acknowledge the financial support of CNPq (National Council for Research and Development) under grants 140664/2013-0, 304121/2013-4 and 200404/2014 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: Finite element formulation
FEM is based on a weak formulation to approximate solutions for a set of PDEs by satisfying an integral problem for a family of test functions, v. Consider the domain Ω ⊂ R n and the Hilbert space H 1 (Ω). The weak formulation of L(w) is given by its actuation on the test function v as (Karniadakis et al. 2005): For axisymmetric domains, it is multiplied by 2πr and given by: where r is the radial coordinate. On Γ N ⊂ ∂Ω, where Robin boundary condition applies, the weak term for axisymmetric domains is given by Karniadakis et al. (2005): Considering the Dirichlet boundary condition given by: the unknown w belongs to the space W defined as: in which w D ∈ H 1/2 (Ω).

A.1 Adsorption problem
From (17), (23) and (13), the following weak forms are obtained for mass and energy balances and for the adsorption rate: Given the convective heat transfer on the walls: The weak forms of the adsorption problem are approximated over a finite element, ω, by the following expressions (Hess et al. 2016): where P and T are, respectively, the vectors with the nodal values of pressure and temperature, N and M denotes the number of nodes used for each unknown and φ and ψ are the shape functions, which are assembled in matrices and . Analogously, the mass of adsorbed gas per volume of adsorbent (Q) is interpolated over finite elements by the shape function ξ : The interpolation of the total porosity (ε) and all other material properties is also done by using ξ : Making v M = φ, v E = ψ and v D = ξ , the FEM system is assembled based on the following coefficients:

A.2 Helmholtz filtering
The weak form of the modified Helmholtz-type PDE is given by: The interpolation of the pseudo-porosities is done over a finite element by using the shape function ξ : Making v H = ξ , the FEM system is assembled based on the following coefficients:

A.3 Derivation of the adjoint problem
In the optimisation problem given in (30), the functional of interest J a (w(e), e) is constrained to the PDE system F(w, e) = 0 given in (33). The difference in the total gas mass between the beginning and the end of the adsorption cycle can be expressed in terms of Dirac's delta inside a time-integral, resulting in the following augmented Lagrangian form: (74) where B denotes the total gas mass, given by: and, for simplicity: where δ denotes Dirac's delta operator. The vector of Lagrange multipliers, λ, is given by: The derivative of the reduced functional J in respect to a design variable e j is given by: Integrating the term involving dẇ/de j by parts, results in: where: Regarding the derivatives of w at the initial time t ia , although P ia and T ia are imposed and do not depend on the design variables, the initial density of gas adsorbed, Q ia , is dependent on the volume of adsorbent material inside the domain. Therefore, the derivative of the state variables in respect to e j is given by: where the equilibrium adsorbed density Q eq is given in (9) and the material model interpolation function m(e j ) is given in (28).
Since the term dw/de j is unknown for t > t ia , the terms b,f a and adj cannot be evaluated explicitly. Given that λ is arbitrary, it can be chosen in order to nullify these terms (Dahl et al. 2008). To ensure b,f a = 0 it is necessary to choose λ such that λ(t f a ) = 0. The term adj is null when the expression between brackets is equal to zero, which can be rewritten as: resulting in: The Lagrange multiplier λ is then given by the solution of the following adjoint problem: where, for simplicity: Such system poses a terminal value problem, which can be converted to an initial value problem by substituting the time variable t per τ = t f a −t. A new system is then defined in terms ofλ(τ ) = λ(t f a − t): (F t (w)) Tλ (τ ) + (F u (w)) Tλ (τ ) = δ τ f a − δ τ ia ∂B ∂w λ(τ ia ) = 0 (89) Letting the Laplace transform ofλ(τ ) be expressed as: The transformation of first two terms in the (89) gives: From the definition of the Laplace transform, the transformation of the last term is given by: which results in: Following, the inverse Laplace transform must be performed to go back to the reverse time domain τ , which for the last term is given by: where H denotes the heaviside step function, which is applied to the function B shifted by τ f a and τ ia . The adjoint problem given in (89) becomes:

Appendix C: Material properties
Material properties adopted in the examples presented in Section 5 are listed here. The properties of methane and activated carbon are given in Tables 8 and 9, respectively. Parameters related to the adsorption of methane on activated   carbon are given in Table 10. Optimised topologies are based on the distribution of steel across the adsorbent bed, which properties are given in Table 11.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.