A metafluid with multistable density and internal energy states

Investigating and tailoring the thermodynamic properties of different fluids is crucial to many fields. For example, the efficiency, operation range, and environmental safety of applications in energy and refrigeration cycles are highly affected by the properties of the respective available fluids. Here, we suggest combining gas, liquid and multistable elastic capsules to create an artificial fluid with a multitude of stable states. We study, theoretically and experimentally, the suspension’s internal energy, equilibrium pressure-density relations, and their stability for both adiabatic and isothermal processes. We show that the elastic multistability of the capsules endows the fluid with multistable thermodynamic properties, including the ability of capturing and storing energy at standard atmospheric conditions, not found in naturally available fluids.

T he thermodynamic properties of fluids are crucial to many fields, and specifically to energy and refrigeration cycles 1,2 . These widespread and essential cycles are leading causes for global warming [3][4][5][6] , and in particularly, refrigeration was recently listed as the single most polluting technology 3 , due to fluids used within such cycles. Creating a fluid with exceptional properties can thus be of great practical importance, yielding possibilities of advancement. In this work, we study the properties of a suspension composed of a multitude of lubricated multistable capsules (structures capable of transforming between different equilibrium deformation patterns 7,8 ) enclosing a compressible gas and immersed within another fluid. The thermodynamic properties of the suspension are determined by an interaction between the external liquid, the encapsulated gas, and the characteristic elastic energy profile of the capsules. By leveraging the elastic multistability of the suspended capsules, we can produce a fluid with multiple stable density points for a given pressure and temperature states as well as unstable regions.
Multistable structures are a special class of mechanical metamaterials. Meta-materials are architected structures made from assemblies of micro-level building blocks, usually arranged in a repeated pattern. While made from standard materials, the geometry of their building blocks can give rise to unique behaviors 9 . For example, meta-materials have been designed to have "negative mass" [10][11][12] , to manipulate light 13,14 and stress waves 15-17 , or to provide ultra-high stiffness to weight ratio 18 . Multistable metamaterials feature a multitude of possible equilibrium configurations for a prescribed load. These meta-stable configurations differ in the state of each mechanical unit (building block), which leads to a highly inhomogeneous and non-affine deformation field for the whole body 19 . Hence, by careful design of their building blocks, these multi-stable structures can be "programmed" to undergo morphological changes in response to external stimuli. The ability to manufacture multi-stable metamaterials at all scales has opened exciting possibilities in a range of applications such as vibration isolation, super-elastic behavior, sensing and actuation, soft robotics, deployable structures and more [20][21][22][23][24][25][26][27][28][29] .
In this work, we adopt a similar approach for the creation of multistable metafluids. Metafluids are artificial suspensions that include a mixture of a standard fluid and metamaterial particles engineered to give rise to properties/behaviors that are not found in natural fluids. We study the pressure-density relations, and stability, of a "metafluid", composed of a multitude of multistable gas-filled capsules suspended in a liquid. A simple example of a multistable capsule, used in the current study, is presented in Fig. 1, showing multiple connected bi-stable elements, produced by combining two non-identical conical frusta in opposing orientations. Bi-stability is therefore achieved by switching between extended and collapsed states through inversion of one frustum which creates a significant change in volume. Below, we fully define our model problem, discuss our approach for energy calculation, present our results for stable and unstable equilibria, present experimental data of pressure cycles, examine the use of metafluid in a refrigeration cycle and provide concluding remarks.

Results
We study the thermodynamic pressure-density relations and their stability for a fluidic suspension composed of multistable elastic capsules, enclosing a compressible gas and immersed in a liquid, as seen in Fig. 1. Each capsule is composed of n identical bi-stable elastic elements creating a single encapsulated volume, filled in with a compressible gas. For the sake of simplicity, body forces are hereafter neglected (see detailed analysis in Supplementary Note 4 and Supplementary Movie 3). We examine a control volume which includes a large number of capsules, so that homogenization may be applied, while keeping all physical properties spatially uniform. The density of the external liquid is denoted by ρ l . The number of capsules per-unit-mass is denoted by Φ. The total mass of the capsules per suspension mass m can be calculated via Φ, as m = Φm c , where m c is the mass of a single capsule. We define Ψ as the external fluid mass ratio. Thus, the mass of the external liquid is m l = Ψm and the total mass of the capsules is m(1 − Ψ). Thus we obtain a Φ-Ψ relation of Since the external fluid is assumed to be incompressible and thus its volume is constant, the equivalent density ρ of the suspension can be expressed by where v i is the gas-filled volume contained within bi-stable element i, and v a is all the additional gas and solid volume of the capsule. The volume v i is a multi-valued function of the difference between the pressure of the internal gas, denoted as p gas and the pressure in the external fluid, denoted simply as p. We thus define this pressure difference by p el = p gas − p. We approximate v i by the tri-linear model (see Fig. 1b) which is a common and useful approximation that enables a simplified analytic treatment of bistability. More specifically, according to this model, when considering the energy relation to the pressure p el and element These coefficients are commonly evaluated by fitting to an experiment 24,35,36 . Thus, the pressure-volume function for each of the bi-stable elements is given by, where k s is determined by requiring continuity between the spinodal region and phases I and II, yielding Þ=ðv II s À v I s Þ < 0. Next, we address the total gas-filled volume of a capsule, which is composed of the sum of all volumes of the n bi-stable elements (given by ∑ n i¼1 v i and additional gas-filled volume v a,gas ). Note that the pressure acting on all elements is identical and depends only on the sum of all volumes, so that the overall pressurevolume relation of the capsule is dictated by the number of elements in phase I (n I ), in phase II (n II ), and in the spinodal (n s ). For convenience we denote the vector containing the number of elements in each multistable state as the permutation, In other words, different permutations lead to different pressure-volume relations which are independent of the specific arrangement within a permutation. Thus, a capsule with n bistable elements has (n 2 + 3n + 2)/2 possible permutations. In the experiment, each capsule contains n = 18 elements which yields 190 possible permutations. Considering that the pressure acting on all elements is identical and equal to p el , we can calculate the volume of an element in each phase according to Eq. (2). Denoting the volume in each phase by v I (p el ), v s (p el ), v II (p el ), we obtain the total volume of gas within a capsule as per ÀÀ! Áfv I ðp el Þ; v s ðp el Þ; v II ðp el Þg þ v a;gas .
Internal energy. We now calculate the internal energy for all of the possible permutations {n I , n s , n II } of the suspension. For the current analysis, we neglect the compressibility of the external fluid and assume that all capsules are identical. Thus, the internal energy of the suspension is composed of the total energy of the ideal gas, denoted by U gas , and the energy of the multistable elastic capsules, denoted by U el .
The elastic energy U el of a bi-stable element i for the three possible phases is given by integration of (2) over v i , yielding and the total elastic potential energy per-unit-mass is thus To calculate U gas we assume an ideal gas model, p gas = ρ gas RT. For isothermal processes, we thus obtain the relation p gas ¼ m g RT=ðv a;gas þ per ÀÀ! Áfv I ; v s ; v II gÞ, where R is the specific gas constant, T is the constant temperature, m g is the mass of the gas within the capsule. For adiabatic processes, assuming ideal gas and isentropic compression/expansion, the gas pressure is given by p gas =p atm ¼ ½v atm =ðv a;gas þ per ÀÀ! Áfv I ; v s ; v II gÞ γ , where p atm is atmospheric pressure, and v atm is the volume taken by the gas contained within a single capsule at atmospheric pressure. In the adiabatic case, the temperature of the gas within the capsules will change during pressure variations according to T ¼ T 0 ðp=p 0 Þ 1À1=γ , but not the temperature of the external incompressible fluid. These spatial temperature variations may yield heat transfer between the gas, solid and liquid, and more complex dynamics. To allow simple treatment of this case, we assume in the current analysis that an adiabatic process corresponds to the more strict requirement of no heat transfer from the encapsulated gas.
The gas energy per-unit-mass of the suspension is defined as and thus the total energy of the system for an adiabatic process is given by U tot = U gas + U el + U ex , where the last RHS term is the external work energy U ex and p is the dictated pressure in the external fluid. For an isothermal process the first RHS term is replaced with Φp atm ln Â ðv a;gas þ per ÀÀ! Áfv I ; v s ; v II gÞ=v atm Ã . Equilibrium of the system occurs when the elastic force in all bi-stable elements is balanced by the difference between the external pressure, and the pressure of the internal gas. In the current analysis, all bi-stable elements are identical, and can be in one of three possible states (phase I, spinodal, and phase II). Thus, the equilibrium values in each of the three possible states can be calculated from the following relations, Equation (8) describes an adiabatic process, and for an isothermal process the first RHS term is replaced with ðm g RTÞ= ðv a;gas þ per ÀÀ! Áfv I ; v s ; v II gÞ.
Solving Eq. (8) for the equilibrium values v I , v s and v II , and substituting these into the energy equation (7) allows to calculate the energy of the suspension per given permutation per ÀÀ! ¼ fn I ; n s ; n II g and external pressure p. Each permutation is only valid in a specific range of external pressures, as defined in Eq. (4). In Fig. 2 we present all of the available energy equilibria states, for multistable capsules composed of 18 bi-stable elements, with 190 possible permutations. Both isothermal (orange lines) and adiabatic (blue lines) processes are presented. The lines are determined by changing the external pressure, for different permutations of the phases of the bi-stable elements. The results show that the multi-stability extends the possible state of an ideal gas from a single line to a region containing multitude of possible energies. Different permutations provide different energypressure curves, corresponding to a fluid with different properties. Changing the properties of the fluid can be achieved by increasing/reducing the pressure beyond the range of a specific permutation. In addition to changing fluid properties, such transitions between different equilibrium states allow the fluid to store and release energy, as discussed below.
Stability and the equation of state. We examine the stability of equilibrium states by considering the Hessian of the total energy U tot , given by H ij = ∂ 2 U tot /∂v i ∂v j . An equilibrium configuration is stable if it corresponds to a local minimum of the energy, requiring positive definite Hessian. In order to obtain the Hessian in terms of all n degrees of freedom, v i , we note that per ÀÀ! Áfv I ; v s ; v II g ¼ ∑ n j¼1 v j . Substituting this relation into Eq. (7) and differentiating the total energy with respect to the degrees of freedom, v i , yields the following set of n equations where α(i) denotes the phase of bi-stable element i, α(i) ∈ {I, s, II}.
The set of n equations in (9) is differentiated again in order to obtain the Hessian matrix, where k α(i) (i = 1, 2, …,n) is the stiffness of the i-th bistable element and k g is the (instantaneous) stiffness of the gas at the equilibrium configuration. For an adiabatic process and for iso- An equilibrium state is stable if and only if the Hessian is positive definite, which occurs if and only if all leading principal minors A j (the determinant of the submatrix of first j rows and j columns) are positive. The principal minors of the current configuration in (10) are In order to obtain conditions for positive A j , we examine separately three possible permutation states. In the first state the permutation is of the form (n I , 0, n II ) and thus there are no elements in the spinodal region. In this case, since all k α(i) and k g are positive, all minors A j in (11) are necessarily positive for all p. Another possible permutation state is the case of two or more elements in the spinodal region. The symmetry in terms of the degrees of freedom allows us to arbitrarily rearrange the degrees of freedom. We re-order the Hessian so that the first two elements are in the spinodal region, and thus k α(1) = k α(2) = k s < 0. The first minor in this case is given by A 1 = k g k s (1/k g + 1/k s ) and is positive only if k g > |k s |. The second minor A 2 ¼ k g ðk s Þ 2 ð1=k g þ 2=k s Þ) is positive only if k g < |k s |/2. Thus, since the two conditions cannot hold simultaneously, if there are at least two elements in the spinodal region, the configuration is unstable for all k g . Finally, we examine the case of permutations of the form (n I , 1, n II ), meaning that only a single element is in the spinodal region. We re-order the Hessian so that the n-th element is in the spinodal region. Thus all minors A j , j ∈ {1, …, n − 1}, are positive. The requirement that the minor A n is positive, which is needed for the stability condition to be satisfied, may be expressed for an adiabatic process as where {v I , v s , v II } are functions of the external pressure p, and are calculated from Eq. (8). For isothermal process the first LHS term is replaced with v a;gas þ per ÀÀ! Áfv I ; v s ; v II g 2 =ðv atm p atm Þ. Obtaining the stability condition allows us to present all possible densitypressure relations which correspond to all permutations of the form per ÀÀ! ¼ ðn I ; 0; n II Þ or permutations of the form per ÀÀ! ¼ ðn I ; 1; n II Þ along with the condition (12). The pressure-density relations of the fluid are presented in Fig. 3, where unstable states are marked by dashed lines and stable states are marked by solid lines. In Fig. 3a we present the external pressure p vs density ρ for adiabatic (blue lines) and isothermal (orange lines) processes. For a metafluid, the adiabatic process is associated with rapid actuation due to a limited heat transfer from the encapsulated gas to the surrounding liquid. For longer time periods, the gas temperature, surrounding liquid temperature, as well as the outside environment will balance, with an isothermal limit being appropriate.
In addition, pressure-density relations for gas only are presented by dotted lines with corresponding colors. The distance between stable points of the fluid is minimal near the permutation per ÀÀ! ¼ ð0; 0; nÞ and increases as we approach the reversed order permutation per ÀÀ! ¼ ðn; 0; 0Þ. In Fig. 3a, we show adiabatic processes, which were obtained for three different initial gas volumes within the capsules at atmospheric conditions (blue line is the reference curve, purple line has a decreased amount of gas, and the green line have an increased amount of gas). The effect of the amount of gas becomes more pronounced as the density increases and the average slope of the multi-stable curve increases with the addition of gas into the capsules. In addition, it is possible to see that permutations of the form per ÀÀ! ¼ ðn I ; 1; n II Þ become stable as the density of the fluid increases (see insets within Fig. 3b). The red curve in Fig. 3b, presents modified k II . In this case, a permutation with unstable element (n s = 1) results in a stable state of the system. The insets show that decreasing v atm or increasing k II creates regions where the fluid is stable although a bi-stable element is at the spinodal region. In both panels, several stable densities are possible for a given pressure, as well as several pressures for a single density. This feature appears for unstable curves with positive slopes.
Pressure cycles and experimental demonstration. In Fig. 4 we present analytic and experimental results for cyclical pressurization of the fluid. In Fig. 4a, each of the black curves denotes a stable equilibrium state corresponding to a different permutation, while the blue and red curves follow the actual loading and unloading paths, respectively, involving switching between different permutations. In order to access one of the multiple possible densities per a given pressure, we need to reach the specific permutation on which the desired density-pressure point is   (9,0,9) per (0,0, 18) per (12,0,6) per (    located. This can be achieved by increasing the pressure beyond the initial snap-close value associated with the current permutation (blue line in Fig. 4a), until reaching a stable permutation for the given pressure value. At this point, reducing the pressure will not change the permutation as long as the pressure is above the snap-open value of the current permutation (red line in Fig. 4a) and the system will remain at the same permutation (following the black line in Fig. 4a associated with that permutation). Thus, multi-stable pressure cycles are closed curves which involve intervals on constant permutation curves (black lines), intervals on the snap-close curves (blue lines), and intervals on the snapopen curves (red lines). This unique ability to choose a specific operation curve is demonstrated experimentally in what follows.
The experimental setup (Fig. 4c) consists of a cuboid 35 × 30 × 3.2 cm fluid-filled tank, a pressure flow controller, and a fluid reservoir. The tank contains 16 capsules, and is connected to a fluidic reservoir in which the pressure is regulated using a pressure controller (Elveflow OB1). Using this setup, we experimentally measured the parameters values of the a single capsule by eliminating the encapsulated gas effects (for more details, see Supplementary Notes 1). The following physical parameters were obtained: v I s ¼ 0:26 ml, v II s ¼ 1:04 ml, k I = 2300 kPa ml −1 , k II = 152 kPa ml −1 and each capsule is composed of n = 18 bistable elements. In addition, we experimentally measured the snapping pressures (see definitions in The surrounding fluid is water (ρ l = 0.997 g ml −1 ) and the internal fluid is air (γ =1.4, R = 0.286 J g −1 K −1 ). To obtain the change in density we measured the height of the fluid in the reservoir, which is equivalent to the change in the volume of the capsules (since the surrounding fluid is assumed incompressible). Throughout the experiment, the height of the liquid in the fluid reservoir and the pressure in the tank were continuously monitored. By using the density of the surrounding fluid and the volume of the metafluid tank, we calculated the equivalent density as a function of the pressure. The increase of tank volume due to elasticity of the walls was accounted for by pressurizing the tank without capsules, and measuring the volume-pressure relation.
In Fig. 4a, the orange and green curves present experimental results of two different pressure-density cycles of a multi-stable fluid. Each experimental curve presents three identical pressure cycles, repeated in each experiment. In Experiment #1, the initial state of the capsules at atmospheric conditions yields equivalent density of ρ 0 = 0.719 g ml −1 and permutation of per Fig. 4a1). We pressurized the tank linearly from p 0 = 101 kPa to 171 kPa over a period of 30 s. The fluid density increased and the permutations changed until reaching a steady state permutation of per ÀÀ! ¼ f14; 0; 4g and density of ρ = 1.07ρ 0 (see Fig. 4a2). At this stage, we reduced the pressure back to p 0 over a period of 30 seconds. In agreement with the model, the permutation of the fluid remained nearly unchanged (per ÀÀ! ¼ f13; 0; 5g) while the density decreased to ρ = 1.037ρ 0 (see Fig. 4a3). To return to the initial state, we applied a pressure of 30 kPa, which allowed us to change the state back to the initial fluid permutation of per ÀÀ! ¼ f5; 0; 13g and the density of ρ = 0.95ρ 0 (see Fig. 4a4). Finally, increasing the pressure back to 101 kPa compressed the fluid along a constant permutation. Similarly, in Experiment #2 we applied a smaller cycle, in which points 1-4 have the densities of ρ = ρ 0 , ρ = 1.04ρ 0 , ρ = 0.97ρ 0 , ρ = 1.03ρ 0 and the permutations of {7, 0, 11}, {11, 0, 7}, {11, 0, 7}, {7, 0, 11}, respectively. As it is evident from Fig. 4, a good agreement is observed between the theoretic stability states of the system and the experimental results for both cycles (see Supplementary Movies 1-2). In Fig. 4b we present the internal energy associated with the different equilibrium states. The internal energy stored within the fluid at atmospheric conditions in the initial state, point 1, is released by adding sufficient energy to the fluid reaching a different permutation, point 2. From point 2, energy is released along a single permutation until reaching point 3, which is the minimal energy at the atmospheric conditions. Similarly, storing energy within the fluid can be achieved by temporarily adding energy to change permutations, and reaching point 4. Reducing the pressure back to atmospheric pressure will not change the fluid permutation, and some of the energy will thus remain stored within the fluid.
Refrigeration cycles. Transfer of a metafluid between different permutations is analogous to phase-change in standard fluids. When the proposed metafluid snaps from one permutation to another, the internal gas rapidly and significantly expands or contracts with constant pressure. When a capsule is snapped open, the encapsulated gas cools down, and this translates to snap in heat absorption of the capsule, which is equivalent to the latent heat associated with evaporation. Similarly, snap-close is equivalent to condensation. In refrigeration cycles, phase-change is essential to achieving efficient cooling cycles 37,38 . As mentioned in the introduction, this often necessitates the use of environmentally hazardous gases with appropriate thermodynamic properties 39,40 . In this section, we discuss and characterize the operation of refrigeration cycles which utilize metafluids, with the aim of enabling the use of standard safe gases such as air for effective refrigeration without requiring phase change.
Cooling cycles in which the fluid remains in a gaseous state, such as the Brayton cycle, can operate with a variety of fluids, with air being the most commonly used. Such cycles, known as air cycle refrigeration (ACR), were for many decades the predominant technology for air conditioning of jet enginedriven airplanes, due to their relatively simple working mechanism and compactness. However, their performance is quite limited in comparison to vapor-compression cycles 41 , in which the refrigerant undergoes phase change throughout the cycle presenting enhanced performance thanks to the ability to absorb additional heat in the form of latent heat.
The thermodynamic cycle of a standard, gas only, reversed Brayton refrigerator 42 is presented in Fig. 5 (solid gray lines). This cycle is used to remove heat from a conditioned reservoir, state 1, where the temperature is T 1 = T c and the pressure is P 1 = P atm . The heat is transferred to the surrounding reservoir, state 3, where the temperature is T 3 = T room and the pressure is P 3 = αP 1 (α > 1 is the compression ratio). Let us assume that the gas enters the compressor at state 1. The gas is then compressed to state 2, where the pressure is P 2 = αP 1 , in an adiabatic and reversible process which causes the gas to heat up, reaching T 2 , which is well above the surrounding temperature, T 2 > T 3 . During this process, the pressure-volume relation follows from the equation P ¼ P 1 V 1 =V À Á γ and the temperature satisfies T ¼ T 1 ðV=V 1 Þ 1Àγ .
After stage 2 is reached, the refrigerant comes in contact with the surrounding reservoir and is cooled down to state 3 in an isobaric process. Next, the gas enters a turbine, which causes the gas to expand back to P 4 = P 1 , reaching the temperature of T 4 which is below that of the conditioned region. After stage 4 is reached, the refrigerant comes in contact with the conditioned reservoir and is heated up to state 1 in an isobaric process. The diagrams of a metafluid cooling cycle are presented in Fig. 5, along with those of the standard Brayton cycle using air. In Fig. 5a, we show a pressure volume diagram and in Fig. 5b we show a temperature-entropy diagram. The blue, green, red, and orange lines in Fig. 5 correspond to metafluid transitions 1-2, 2-3, 3-4, and 4-1, respectively. Gray lines denote standard air based Brayton cycle. The properties of the metafluid are the same as those presented in the previous section, and are based upon experimentally measured values. The externally applied pressure is different from the pressure of the gas, and isentropic relations yield where p el depends on the permutation of the metafluid, per ÀÀ! ¼ ðn I ; n s ; n II Þ, and is calculated via Eq. (2).
During the entire cycle, multiple snaps of the metafluid occur between different permutations. In order to quantify the efficiency of a metafluid-based cooling cycle, the coefficient of performance (COP) is calculated. The COP is defined as the ratio of the refrigeration effect to the net work input. Net work is calculated as the work performed on the compressor less the work performed by the turbine, and refrigeration effect is measured as heat transferred between states 4 and 1. The COP is given by See Supplementary Notes 2 for a detailed description of the calculation. Based on the exact parameters used in the experiment (see Fig. 4), we found that the metafluid demonstrated COP of 2.29, which is significantly higher than the Brayton air based cooling cycle's COP of 1.4. The efficiency improvement of 61% was achieved without any optimization of the system. In Supplementary Notes 3 we present the effect of various parameters, such as the number of bi-stable elements in the capsule and the mass of the encapsulated gas. In the figures presented in the Supplementary Notes, all of the system's parameters except for one are held constant. The results show that the efficiency of the system (COP) is mainly determined by the number of snaps occurring during the heat exchanging segments of the refrigeration cycle. Further optimization can be achieved if we modify the properties of the bi-stable frusta, showing optimal value of about 0:3k exp I in which there is COP value of 8.28, demonstrating that metafluids have the potential to replace vapor compression cycles and enable safe and efficient refrigeration while using common and safe gases, such as air.

Discussion
In this study, we introduce a new concept for a multistable metafluid, and investigate theoretically and experimentally the relationship between pressure and density both for adiabatic and isothermal processes. As a result of the variety of stable densities and energy states, it is possible to capture and store energy at standard atmospheric conditions, which is not possible with naturally occurring fluids. Due to these properties, metafluids are relevant to the ubiquitous refrigeration and heat cycles, as well as many other applications such as the creation of flow fields with physical memory attributed, as demonstrated in recent studies in the context of metamaterials 43,44 .
The experiments presented in this work involved large O(cm) particles. While these experiments are aimed at demonstrating the concept, future research on miniaturization of the multistable particles is needed for creating metafluids for real applications. Minimization of the capsules may be achieved by utilizing a wide range of established manufacturing technologies 45 , such as lithography 46,47 , high-resolution additive manufacturing 48,49 , and propagating photo-polymer wave-guide prototyping 50 , which were proven to enable bi-stable or multistable elements.
We examined the application of such a metafluid to refrigeration cycles and presented significant improvements to the efficiency of such cycles. We note however, that in order to realize a refrigeration device based on a metafluid, various other requirements need to be taken into account, such as optimizing heat transfer rates, as well as creating a system that does not damage the capsules. These considerations may alter the COP values for an actual device, and will need to be addressed in a future research. The results of this work suggest that the concept of multistable metafluids paves the way for the creation of futuristic non-polluting fluids with enhanced properties for energy and cooling cycles, leveraging multi-stability to harvest, store and release energy and heat.

Methods
Details on the geometry, design, fabrication, testing, analytical model, and numerical solutions are provided in the Supplementary Information.