Performance optimization of bijels as a novel type of catalyst support structure

Abstract

Due to their high surface-to-volume ratio, porous media are very suitable catalyst support materials. However, the stochastic morphology of commercially available supports generally results in poor transport properties. Accordingly, the catalyst material contained in the support structures is used inefficiently. These issues can be alleviated by using catalyst supports acquired from spinodally-derived architectures due to their beneficial percolation properties. In particular, architectures derived from bicontinuous interfacially jammed emulsion gels (bijels) provide a viable route to the manufacture of stable catalyst supports that address the aforementioned issues. We investigate this type of porous support using a hybrid simulation method that combines a multicomponent lattice Boltzmann solver with the discrete element method. First, we simulate the formation of bijels and report on the improved properties of the resulting porous structures compared to stochastic equivalents for use in chemical reactors. Hereafter, we further validate the enhanced performance of bijel-derived geometries through miscible reactive flow simulations. Our findings suggest that bijel-derived catalyst support structures allow for an almost threefold increase in reactor effectiveness.

---

New concepts

For the first time, we demonstrate – purely based on computer simulations – how bijel-derived porous materials with distinct morphological properties can be obtained. Due to the bicontinuous morphology of bijels, the porous materials are percolating by definition, i.e., there are no dead ends or closed pores. These properties make such materials perfect candidates for porous materials with a wide spectrum of applications. Here, we demonstrate their possible impact on the development of new catalyst support materials. We quantify the superior morphology of the bijel-derived materials as compared to classical stochastic catalyst supports and show with the help of reactive flow simulations that bijel-derived catalyst support structures enable an almost threefold increase in reactor effectiveness compared to generally used stochastic supports. Our multiscale concept involving simulations for the material generation, the analysis of the material properties, as well as for the demonstration of their applicability in real devices, is unique and opens a possible new way for the development of new materials with well-specified properties required for distinct applications.

1. Introduction

Heterogeneous catalysis is used to produce almost every chemical at some point in its production process.¹ Given the large impact catalysis has on chemical and energy industries, modeling and optimizing catalysis processes can be very worthwhile from an economic and environmental point of view.² Porous structures are generally used in catalytic devices due to their very large surface-to-volume ratio. However, the stochastic nature of commercially available catalyst support structures has several disadvantages. Irregular pore sizes and arrangements can result in structures with nonhomogeneous pore size distributions and dead pores. This, in turn, can lead to a loss of surface area and impede mass transport, rendering the search for improved support structures viable.

To improve the catalytic efficiency of porous support structures, a high flow permeability is required. Hence, the pore space should ideally be fully connected and not have any constrictions or dead ends.³ Additionally, for the sake of the stability of the structure, full solid-state connectivity also has to be maintained. Porous structures realized from spinodally decomposing fluids generally possess these properties.⁴ In particular, adding particles that arrest the phase separation to obtain bicontinuous interfacially jammed emulsion gels (bijels) seems to be a viable route to obtain these kinds of structures. Solidification of one of the two fluid phases then leads to a structure that has by definition no dead ends, no isolated pores, and uniform domain sizes.

If particles with partial wettability are present between a spinodally decomposing mixture of two immiscible fluids, bijels can be formed. This is possible since the free energy of the total system is reduced if the particle remains in contact with both fluids.⁵ This energy well is generally several orders of magnitude higher than the thermal energy, leading to the particles becoming trapped at the interface. As the phase separation of the fluids continues, the amount of available surface area for the particles steadily decreases. This forces the particles into close contact to the point where they arrest further interfacial coarsening. At this point, they formed either droplets covered with particles, known as a Pickering emulsion,⁶ or a quasi two-dimensional sheet percolating through the entire domain, known as a bijel.

The existence of bijels was first predicted in 2005 by Stratford et al.⁷ by means of lattice Boltzmann method simulations. Two years later, their existence was experimentally confirmed by Herzig et al.⁸ Since this pioneering work, the formation of bijels has been achieved using multiple fabrication methods for a variety of fluid systems and colloidal particle types.^9–22 In addition, several computational studies have further explored the design space and properties of bijels. The influence of particle wetting, particle volume fraction, and fluid component ratio on the transition from bijels to Pickering emulsions was further investigated by Jansen and Harting.²³ Using the same pseudopotential lattice Boltzmann model coupled with the discrete element method, Frijters et al.²⁴ further studied the domain sizes of particle-stabilized emulsions. Additionally, there have been numerical studies on bijel formation with magnetic²⁵ and anisotropic particles.²⁶ Furthermore, Carmack and Millett studied both bijel formation in confined regions²⁷ as well as their tunability in thin films using external electric fields.²⁸ The shear resistance of bijels was investigated by Bonaccorso et al.²⁹ Besides, phase-field simulations have been used to simulate the formation of bijels more accurately.³⁰ Finally, it is worth noting that recent developments in simulation techniques have reduced computational effort for simulating large numbers of particles at interfaces.³¹

In recent years, there has been an increasing interest in applying bijels for catalysis purposes, sparked by the successful generation of porous media from bijels in their arrested state. The possibility for solidification was first experimentally demonstrated by Lee and Mohraz.³² Subsequently, Witt et al.¹² demonstrated the large tunable range of domain sizes of the solidified structure, and Reeves et al.³³ worked on quantifying the resulting bijel topology. This morphological characterization was further extended by McDevitt et al.³⁴

The usability of bijel-derived morphologies for chemical reactions was experimentally demonstrated by Witt et al.³⁵ The authors synthesized bijel-based electrodes and reported more than a 1.5× improvement in energy density over previously reported electrode morphologies. These findings were further substantiated by Gross et al.,³⁶ who investigated bijel-based electrodes for gas generation reactions. They found that bijel-based electrodes expel product gas faster and require less overpotential to drive electrolysis. Further investigations into the manufacturing of this type of morphology for chemical reactions were done by Allasia et al.³⁷ They demonstrated the robustness and high thermal stability of this structure type by investigating its usability for the photodegradation of methylene blue, further underscoring its possible applicability in industrial settings. While these experimental studies indicate the significant potential of these structures for electrodes and catalysis, to the authors' best knowledge, thorough numerical studies on the optimization of these structures for chemical reactions are still lacking.

The remainder of this paper is organized as follows: first, we introduce the numerical method for our simulations in Section 2. Simulation data of bijel formations is presented in Section 3. In Section 4, an algorithm to generate stochastic packings of overlapping spheres with a predefined porosity and surface area is introduced. The morphologies of the resulting solid geometries are analyzed in Section 5. Hereafter, the improvement of bijel-based support structures with respect to stochastic ones is further substantiated in Section 6 by means of reactive flow simulations. Finally, a conclusion is presented in Section 7.

2. Numerical method

2.1. Multicomponent lattice Boltzmann model

For the simulation of bijel formation and reactive flows, we make use of the lattice Boltzmann method. This method is well established and has proven to be very suitable for modeling systems involving complex fluids.^38,39 In comparison with traditional Navier-Stokes solvers, the presence of complex boundary conditions and interactions is handled without much additional complexity. Additionally, the high degree of locality of the algorithm allows for efficient parallelization on supercomputers.⁴⁰

The lattice Boltzmann method is based on discretizing the Boltzmann equation in physical space, velocity space, and time. Each component of the simulated fluids obeys the lattice Boltzmann equation given by

f_i(x + c_iΔt,t + Δt) = f_i(x,t) + Ω_i(x,t) (1)

where f_i(x,t) is the single-particle distribution function at position x and time t. The lattice vectors are given by c_i and have a weight w_i. Additionally, the lattice constant following from the physical discretization, as well as the time step obtained from the time discretization, follow the general convention Δx = Δt = 1. Hence, is the speed of sound in our simulations. Furthermore, we consider a three-dimensional lattice and a velocity space discretization in 19 distinct vectors that connect neighboring lattice sites (D3Q19). We employ the Bhatnagar-Gross-Krook (BGK) collision operatorwhich allows for reconstructing the behavior predicted by the Navier-Stokes equations. Here, the relaxation time τ defines the kinematic viscosity as . The populations are generally relaxed towards the second-order truncated equilibrium distributionwith ρ the local density and u the velocity vector, which can be computed from the zeroth and first moments of the populations

On the other hand, in the cases where particles are included in the simulations, the third-order truncated equilibrium distribution is used to improve the algorithm's stability.⁴⁰

Only single-component fluids and binary fluid mixtures will be considered in this paper. We will distinguish the fluids for the binary fluid mixtures by referring to them as red fluid (denoted by the superscript “r”) and blue fluid (denoted by the superscript “b”). Additionally, we define a color field as ζ(x,t) = ρ^r(x,t) − ρ^b(x,t). The surface tension between the red and the blue fluid is introduced by means of the pseudopotential model of Shan and Chen.^41,42 For the red fluid, this can be written as

where g_rb is a force coupling constant and Ψ^r(x,t) the effective mass given bywith ρ^r₀ the initial density for the red fluid. The equations for the blue fluid can be obtained by swapping the red and blue indices due to symmetry considerations. Since these interactions only involve the nearest neighbors, they allow for the inclusion of repulsive forces between the two fluids without sacrificing much of the locality of the algorithm.⁴²

2.2. Suspended particles

Fully resolved colloidal particles are modeled in the lattice Boltzmann method based on pioneering work by Ladd and Aidun.^43–45 The particles are discretized on the lattice by marking all the lattice sites inside a particle as solid sites. Their momentum and angular momentum are integrated with a leap-frog integrator and governed by three different types of forces.

First, there are forces due to the exchange of momentum between the particles and the surrounding fluids. These are generally implemented by means of a bounce-back boundary condition for the fluid populations at the particle surface. To this end, the streaming step is locally modified to

f_i(x,t + Δt) = f_ī(x,t) (7)

with f_ī denoting a population moving into the opposite direction of f_i and into a solid node. Furthermore, we use the same modifications to the streaming step, with the momentum exchange coupling for solid particles as described in Jansen and Harting.²³

Second, there are forces due to two different types of potentials. When the distance between two particles is closer than the lattice discretization, their hydrodynamic interactions can no longer be modeled correctly by the lattice Boltzmann method. Then, a lubrication correction

is applied to correct for increasing fluid forces arising from fluid squeezed out of the interparticle gap.⁴⁶ In this equation, μ₀ is the average dynamic viscosity at the beginning of the simulation, r_p the particle radius, Δ_gap = ||x₁₂|| − 2r_p the inter-particle gap, a cutoff distance, u₁₂ = u₁ − u₂ the difference in the particle's velocity vectors and x₁₂ the vector connecting the particle centers. When the lubrication correction cannot avoid the overlap of particles, a force based on the Hertz potential⁴⁷with coupling constant K_Hertz = 100 is added to ensure the overlap is minimized.

Third and finally, there are wetting-related forces. These are created by setting virtual fluid densities for both fluid components on lattice sites covered by particles. The virtual fluid densities can then be used for the computation of wetting-related forces²³ applied like in eqn (5). Setting the ratio between the two virtual fluid densities to a number different than one allows for setting an equilibrium particle-interface contact angle that is different from 90 degrees.²³

For a more extensive description of the method and its implementation details, we refer the reader to our previous papers.^23,24,48

2.3. Boundary conditions

In all the conducted simulations, periodic boundaries are applied to both the fluid and the particles. Additionally, we apply bounce-back boundary conditions as shown in eqn (7) at all solid boundaries. Furthermore, in simulations with reactive flows, we also make use of two other boundary conditions that modify the populations.

First, we model the conversion of red to blue fluid as a first-order kinetic boundary condition. To this end, the fluid densities near the solid boundaries are changed as

andby scaling the populations. In this equation, the reaction rate ξ is positive, ρ^r is the local reactant density, and ρ^b is the local product density. Furthermore, NN(x) is the set of nearest neighbors of lattice node x, and s is an indicator function, being 0 for void sites and 1 for solid sites. This allows for properly accounting for the available surface area by rescaling based on the local number of nearest-neighbor solid nodes. Additionally, to avoid unphysical negative populations in the simulations due to forcing,⁴⁹ the density of the reactant converges towards a positive non-zero equilibrium density ρ^r_eq. An equilibrium density ρ^r_eq of 10⁻⁴ proved to be sufficient to resolve the aforementioned issue.

Second, we employ so-called recoloring layers in the reactive flow simulations to remove the reaction products. These recoloring layers are placed at the cross-sections with the lowest and the highest x-coordinate values. At a recoloring layer, the ratio between the densities is set to the initial density ratio: ρ^r_init/ρ^b_init = 1/(10⁻⁴). Given that the fluid is considered to be only weakly compressible, this effectively sets the concentration of the fluids at these sites equal to the initial concentrations. In the lattice Boltzmann method, this is implemented by multiplying all the populations by numbers that keep the total density constant but change their ratio. For the red fluid, this is done by applying

where the equation for the blue fluid can be obtained by swapping the indices.

3. Bijel formation

The bijel formations studied in this paper are all simulated on a 1024³ lattice. This very large lattice size is selected for two separate reasons. Firstly, a large lattice is required to avoid finite-size effects and to create a representative elementary volume of the real structure with a statistically relevant number of pores. Secondly, this lattice size is the largest that is computationally feasible with the resources at our disposal. The fluid densities are initialized with uniformly distributed random numbers between 0 and 0.5 for the red fluid and between 0 and 0.9 for the blue fluid. This sets an average red fluid density [small rho, Greek, macron]^r = 0.25 and an average blue fluid density [small rho, Greek, macron]^b = 0.45. Subsequently, we have a red-to-blue ratio of 5/9 for all the simulations in this paper. The suspended particles are initialized at randomly distributed non-overlapping positions. With a particle radius r_p = 5Δx and volume fractions Ξ of 0.15, 0.20 and 0.25, our simulations contain between approximately 300 000 and 500 000 particles. The virtual fluid densities of the particles are initialized as ρ^r_p = [small rho, Greek, macron]^r + Δρ^r for the red fluid and ρ^b_p = [small rho, Greek, macron]^b for the blue fluid. We vary Δρ^r between 0.0, 0.015, 0.030, 0.045 and 0.060. This results in the particles having a neutral wetting or a preferred wetting towards the minority (red) fluid component. Furthermore, with the used g_rb = 0.08, this leads to an equilibrium interface-particle contact angle θ_p between 90° and approximately 117°.²³ Additionally, we set τ = 1 for both fluids.

At early simulation times, we observe a period of rapid initial domain coarsening.²³ This domain coarsening starts to slow down when more and more particles accumulate at the interface and inhibit the phase separation. To quantify this behavior, we make use of the three-dimensional structure factor, which can be defined as

S(k,t) = |ζ(x,t)e^ik·x|/V (13)

here, k is a wavevector in the reciprocal space, and V = L³ is the total number of lattice nodes inside the system. We further smoothen S(k,t) by averaging over shells in k space to obtain the spherically averaged structure factorwhere k is the length of the wavevector k and the sum runs over spherical shells defined as in k space. Hereafter, it is possible to compute the characteristic domain size as

The time development of the characteristic domain size is shown in Fig. 1 for the simulations with a particle volume fraction Ξ = 0.20. It can be clearly seen that there is a period of rapid initial domain coarsening, which starts to slow down between time steps 10 000Δt and 20 000Δt. At this point, many particles have accumulated at the interface, and their lubrication and Hertz forces start resisting further coarsening of the fluid–fluid interface. However, in the studied simulations, the particle accumulation at the interface does not fully impede the domain coarsening. This is due to the fact that bijels are not thermodynamically but kinetically stabilized. The detachment energy keeping the particles trapped at the interface

ΔE_det = πr_p²σ(1 − |cos θ_p|)² (16)

is the defining factor for the stability of these structures.^50,51 From eqn (16) it follows that both an increase in the particle radius r_p or surface tension σ allow for an increase in the detachment energy. However, increasing the particle radius increases the computational expense. Additionally, the surface tensions that can be simulated with the pseudopotential lattice Boltzmann method are limited in range and significantly lower than what can be achieved in experiments. Besides, it should be noted that an increase of the equilibrium particle contact angle θ_p also decreases the detachment energy. Moreover, an increase in equilibrium particle contact angle also increases the preferred curvature of the interface, leading to a further increase in pinch-off effects. This can be clearly seen in Fig. 1. When the equilibrium contact angle is increased, particles are more quickly expelled, leading to a more rapid growth of the characteristic domain size R(t).

Fig. 1 The characteristic domain sizes for the bijel simulations with a particle volume fraction Ξ = 0.20. Different equilibrium particle contact angles are shown as a color gradient from blue (θ_p = 90°) to red (θ_p ≈ 117°). After a period of rapid initial domain coarsening, the particles start stabilizing the interface between the two fluids, and the increase in characteristic domain size becomes more moderate. The inset zoom shows the gradient in characteristic domain size. It clearly displays that although the rate of domain coarsening becomes slower, it does not tend to absolute zero for any of the simulations.

For the sake of this study, the bijel simulations are stopped at time step 100 000Δt. At this time, the accumulation of the particles at the interface is finished and the gradients of the characteristic domain size are small for all the simulations. The fact that domain coarsening is still occurring at this time step can be regarded as irrelevant since we are interested in structure generation instead of bijel stability.

To generate porous structures from the simulations, we make use of the color fields. The conversion process is illustrated in Fig. 2. First, all lattice sites with a color larger than or equal to zero become solid sites, while the remainder is set to void sites. Since the color field for the particles is equal to zero, this leads to the solidification of the particles and all the fluid sites where red is the majority component. Second, lattice sites that are not part of the resulting porous structure are handled. Since expelled particles floating around in the void space would be flushed out of the porous structure in experiments, their solid sites should be removed from the simulation domain. To this end, we use the Hoshen–Kopelman algorithm.⁵² This algorithm makes it straightforward to detect clusters of solid sites that do not span the entire simulation domain and should be converted to void sites. Thereafter, the Hoshen–Kopelman algorithm is applied a second time to detect and solidify very small non-percolating void spaces inside the structure related to discretization artifacts of the particles.

Fig. 2 A schematic of the conversion of a color field to a porous geometry. Picture (a) showcases the color field obtained from a small (256³) bijel simulation. Here, the color value of the particles is equal to zero. Picture (b) shows the porous geometry obtained from the thresholded color field. Picture (c) is the final porous geometry obtained from applying the Hoshen–Kopelman algorithm twice on the thresholded color field: once to remove left-over particles in the void space and a second time to solidify small non-percolating void spaces inside the structure.

4. Stochastic geometries

We model the stochastic morphologies of commercially available catalyst supports by using packings of overlapping spheres. These packings can be seen as a model for porous structures based on sintered powders of spherical particles. Hence, they closely resemble the stochastic morphology of porous structures used in certain types of commercially available catalyst supports.

When packings of overlapping spheres are generated on a lattice, they allow for the optimization of different properties based on the number of spheres, the radius of the spheres, and the distribution defining their placement probability. Moreover, if we assume that the sphere centers are located exactly on the lattice nodes and the placement probability is uniform, both the porosity

and the surface areacan be optimized simultaneously.

When starting the algorithm, we initialize a lattice L with uniformly distributed random numbers between zero and one. These numbers determine the location of the sphere centers; lattice nodes with numbers larger than a scalar threshold t are used as center locations. This allows us to define a thresholded lattice as L_t = L > t. For a given threshold t and sphere radius r, a porous geometry can then be obtained by computing the distance to the nearest sphere center for every lattice site. To this end, an Euclidean distance transform (EDT)⁵³ can be employed, which returns the Euclidean distance to the nearest sphere center. Accordingly, all the lattice sites where the distance to the nearest sphere center is smaller than the radius can be converted to porous geometries as L_p = EDT(L_t) < r.

Optimizing t and r is done by an algorithm based on a nested bisection search. Therefore, we perform a search on t and r in a hierarchical manner. We first select lower bounds t_lower, and r_lower and upper bounds t_upper, and r_upper for t and r, respectively. The exact values of the bounds do not matter as long as the interval contains the values giving the required porosity and surface area. Hereafter, we set our initial guesses for t and r equal to the average of the bounds, after which we optimize the porosity by doing a bisection search for the threshold between t_lower and t_upper. This involves computing the porosity several times for the given value of r while adjusting the value of t by changing its lower and upper bounds based on previous under- or overestimates. When t has converged, we compute the surface area, update the bounds for r, and repeat the process until convergence is reached for both t and r.

To properly compare the bijel-based porous structures with stochastically generated ones, we generate one equivalent stochastic structure for every bijel. With equivalent stochastic structure, we denote a stochastic structure with matched porosities and surface area generated by employing the algorithm described above. A visualization of three bijels with their equivalent stochastic structures can be seen in Fig. 3. The deviations between the bijels and the equivalent stochastic structures are always below 0.01% for the porosity and below 0.5% for the surface area.

Fig. 3 A visualization of six large scale (1024³) post-processed geometries. The left column shows the geometries obtained from bijel simulations with a particle volume fraction Ξ = 0.20 (insets a, c, e). The right column shows the equivalent stochastic geometries with matched surface area A and porosity ϕ (insets b, d, f). The equilibrium particle contact angle in the bijel simulations is set to θ_p = 90° (a), θ_p ≈ 103° (c), and θ_p ≈ 117° (e).

5. Porous structure analysis

5.1. Porosity and surface area

The first property studied is the relation between the porosity and the surface-to-volume ratio Λ = A/V. Fig. 4 shows the relationship between these two properties for all the large-scale bijel-derived structures. There is a clear correlation between the available surface area and the equilibrium particle contact angle. This is due to increased particle detachment for larger contact angles and hence decreased surface area, as discussed in Section 3. On the other hand, the porosity does not differ much when the equilibrium particle contact angle is changed. This is because the red-to-blue fluid ratio remains the same for all the simulations. Consequently, the main contribution to changes in the porosity is the differing amounts of particles that are removed since they are not adsorbed to any interface and remain in the fluid phase turned to void. The number of particles is generally not very large since particles have a preferred wetting towards the fluid phase being solidified. Therefore, they are almost always expelled towards this fluid phase, which does not change the final porosity.

Fig. 4 The porosity ϕ plotted against the surface-to-volume ratio Λ for the large scale simulations. The color gradient indicates the equilibrium particle contact angle from blue (θ_p = 90°) to red (θ_p ≈ 117°).

The porosity and surface area values for the equivalent stochastic structures are not shown in Fig. 4. Given that both the porosity and the surface area have been matched to the values of the bijels, their data points lie exactly on top of the displayed ones.

5.2. Pore size distributions

The geometries' pore size distribution (PSD) is further investigated by computing, for every lattice node in the void space, the radius of the largest sphere that fits inside a pore and still covers the node. For this purpose, we use an algorithm as proposed by Münch and Holzer⁵⁴ to compute a pore radius for every void site. Hereafter, we can apply binning to the resulting radii to obtain a representation of the pore size distribution.

Fig. 5 shows the PSDs and the cumulative PSDs for three bijel-derived geometries with a particle volume fraction of Ξ = 0.20. Additionally, their stochastic equivalents are shown too. For the geometries not shown in this figure, it is also observed that the average pore size of the bijels is always larger than the average pore size of their equivalent stochastic geometries. Besides, broader pore size distributions for the stochastic geometries are observed in all cases. In other words, for the bijel-derived geometries, the standard deviation in the pore sizes is significantly smaller for all cases.

Fig. 5 The pore size distributions for three bijels with a particle volume fraction of Ξ = 0.20 and their stochastic equivalents. The plotted equilibrium particle contact angles for the bijels are θ_p = 90° (blue triangles), θ_p ≈ 103° (purple circles) and θ_p ≈ 117° (red squares). The inset shows the corresponding cumulative pore size distributions. Bijels have significantly larger pores than their stochastic equivalents for all the geometries. Additionally, the curves for the cumulative pore size distributions are steeper for the bijels. This indicates that, in these geometries, the standard deviation of pore sizes is smaller than in their stochastic equivalents.

5.3. Mean surface curvature

To estimate the mean surface curvature, we first compute the normal vectors at the fluid-solid interface. We do this by defining a field q(x) in which solid sites are set to zero and fluid sites are set to one. An estimation of the normal vectors is then given ashere, the gradients of the q(x) field are computed using the lattice-weighted average schemewhere χ can be an arbitrary scalar field. With the normal vectors computed, we can compute the mean curvature at every lattice site using

H = −[(I − n ⊗ n)·∇]·n, (21)

where we use the same scheme in eqn (20) to compute the gradients. In Fig. 6, we report on the distribution of the mean curvature for all the sites at the fluid–solid boundary. These sites are defined as fluid sites having at least one lattice vector pointing into a solid site, or solid sites having at least one lattice vector pointing into a fluid site.

Fig. 6 The mean curvature distributions for three bijels with a particle contact angle of θ_p = 90° and their stochastic equivalents. It can be seen that the particle volume fraction, denoted by different colors and markers, has only a minor influence on the distribution. For the bijel, we see an asymmetric distribution in curvatures. While the distribution peaks around H = −0.3, due to the shape of the distribution, its mean is still close to zero, as theoretically expected. By computing a radius from the mean curvature as r_H = −2/H, we observe that the peak at H = −0.3 corresponds to the radius r_p of the suspended particles trapped at the interface. Similarly, for the stochastic structures, the sharp peak in the distribution is directly related to the radius r of the overlapping spheres.

As can be observed in Fig. 6, the mean curvature distribution for the bijels is asymmetric. While there is a peak around H = −0.3, due to the shape of the distribution, the mean of the distribution is still close to zero. This is in line with what can be theoretically expected from a bijel.⁵⁵ The peak around H = −0.3 is directly related to the curvature of the suspended particles. We can retrieve the original radius of the particles by computing a radius from the mean curvature as r_H = −2/H. For the stochastic structures, a similar argument applies. A sharp peak for one particular mean curvature can be observed. This curvature is directly related to the radius of the overlapping spheres out of which these structures consist.

5.4. Tortuosity

To compute the flow tortuosity of the geometries, we again employ the lattice Boltzmann method as described in Section 2. Here, we force a single-component fluid through the porous geometries by applying a body force with an acceleration constant a_x,y,z = 10⁻⁶ in the Cartesian directions to obtain three velocity fields. It was shown by Duda et al.⁵⁶ that such a velocity field for incompressible flow allows for the direct computation of the tortuosity. This idea can be extended to the regular lattice used by the lattice Boltzmann method to approximate the tortuosity along a Cartesian axis as⁵⁷

In this equation, u is the magnitude of the velocity vector, and u_x,y,z is the velocity vector component for one of the three Cartesian directions. Eqn (22) should give accurate approximations if the tortuosity is close to one and the flow does not contain recirculation zones. Given the nature of the geometries studied in this paper, this can be safely assumed. In Fig. 7, the computed tortuosities averaged over the three Cartesian directions are plotted. The two main branches in the data points correspond to the bijel-derived geometries and their equivalent stochastic geometries. For the bijels, we observe a slight increase in the tortuosities when the equilibrium particle contact angle is increased. On the contrary, since the stochastic structures are based upon overlapping sphere packings, their tortuosities remain much more constant over the range of studied porosity and surface area.

Fig. 7 The flow tortuosity averaged along the Cartesian axes 〈T〉 plotted against the equilibrium particle contact angle θ_p. Bijels with small equilibrium particle contact angles tend to have slightly less tortuous flow paths than their equivalent stochastic structures. However, when the contact angle of the particles is increased, the computed tortuosities of the bijels tend to become higher than the tortuosities of the stochastic structures.

5.5. Permeability

To compute the permeabilities of the stochastic structures, we make use of the same methodologies as Narvaez et al.^58,59 First, we rewrite Darcy's lawinto a form that can be used straightforwardly in lattice Boltzmann simulations. We make use of the same simulations that are also used to compute the tortuosity. Therefore, we can substitute expressions for the kinematic viscosity ν and the pressure gradient ∇p in eqn (23). This givesfor the permeability in the three Cartesian directions.

To quantify these geometries as chemical reactors, we postulate that their efficiency depends on the permeability κ and the surface-to-volume ratio A/V. To this end, we define a dimensionless reactor efficiency number as

The averaged permeabilities are plotted against the surface-to-volume ratio in Fig. 8. Additionally, isolines for the reactor efficiency number are shown to guide the eye. It can clearly be seen that there is some spread for the reactor efficiencies of the bijel-derived structures. However, the reactor efficiency numbers for the stochastic structures appear to lie on a single isoline. This coincides with trends predicted by the Kozeny-Carman equation for the permeability

which can be used as a theoretical model to study sintered bead packs.⁶⁰ Given the equivalent porosity and tortuosity for all the stochastic structures, their similar reactor efficiency numbers align with predictions from eqn (26).

Fig. 8 The average permeabilities for the porous geometries along the Cartesian axes 〈κ〉 plotted against the surface-to-volume ratio Λ for the bijels (a) and their stochastic equivalents (b). The equilibrium particle contact angles are shown by a gradient in color from blue (θ_p = 90°) to red (θ_p ≈ 117°). To guide the eye, we define a dimensionless reactor efficiency number as η = κΛ², for which we plot several isolines. For both the bijels and the stochastic geometries, a clear trade-off between permeability and surface area is visible. Additionally, it can be seen that for the simulations where the equilibrium particle contact angle is small, the bijels show significantly better efficiencies.

Furthermore, we plot the permeabilities against the equilibrium particle contact angles in Fig. 9. The permeabilities of the geometries converge towards similar values for larger equilibrium particle contact angles of the bijel. Since the surface area of these geometries is relatively small, the overall domain sizes need to become larger to accommodate the decrease in surface area. This, in turn, leads to fewer isolated pores and fewer pore necks in the stochastic geometries.

Fig. 9 The average permeability for the three coordinate directions 〈κ〉 plotted against the equilibrium particle contact angle θ_p. Standard deviations are so small that they have not been included in the plot. It can clearly be seen that the difference in permeability between bijels and their equivalent stochastic structures becomes smaller when the contact angle increases. The surface areas are smaller for the larger contact angles due to increased particle detachment in the bijel simulations. This smaller surface-to-volume ratio results in fewer isolated pores, fewer pore necks, and larger permeabilities in the equivalent stochastic structures.

6. Reactive flows

The reactor efficiency number, as defined in Section 5, enables us to get initial quantitative estimates of the usefulness of a porous structure as a catalyst support. However, while this number gives us good initial estimates, it only captures effects related to the permeability and surface-to-volume ratio. To this end, we resort to reactive flow simulations with uniform catalyst loading to obtain more accurate estimates for the performance of the geometries as catalyst supports. We have opted for uniform catalyst loading due to the uniformity of the domain sizes of the bijel throughout the domain. Since the uniformity and pore connectivity that can be expected for the bijel-based structures are not guaranteed for stochastic structures, we expect the deviations between these structures to be the most pronounced with this type of loading.

To reduce the computational cost of these reactive flow simulations, the size of the porous geometries is reduced to a quarter of the size of the original geometries in every coordinate direction. A reduced-size geometry is constructed by sampling the original size solid geometry at every fourth lattice site in every Cartesian direction. As a result, there are minor deviations in the expected permeability and surface-to-volume ratio of the subsampled geometries, as can be seen in Fig. 10. While these changes will slightly alter the results quantitatively, they do not have a qualitative impact.

Fig. 10 The percentage deviation from the expected values in average permeability 〈κ〉 and surface-to-volume ratio Λ plotted against the particle contact angle θ_p. The values are obtained from subsampling geometries with a particle volume fraction of Ξ = 0.20 from 1024³ to 256³. We observe minor deviations from the expected values, with the exception being the surface-to-volume ratio of the bijel. In the case of bijels, subsampling leads to a smoothing of the particulate surface. Hence, it leads to a relatively strong decrease in the available surface area. However, this decrease in surface-to-volume ratio is compensated for by an increase in permeability. The values of the curves for different particle volume fractions follow identical trends.

To quantify the performance of the geometries, we make use of the effectiveness factor. This is the amount of reactant converted at steady state divided by the amount of reactant that could have been converted if there were no retardation. For our simulation setup, we can write the effectiveness factor as

6.1. Diffusion-dominated flows

To simulate diffusion-dominated flows, we first initialize the void space with reactant and a minority density of product, ρ^r_init/ρ^b_init = 1/(10⁻⁴). Additionally, we apply the boundary conditions as described in Section 2.3 to convert reactant into product inside the geometry and to fix the concentrations to the initial values at its domain boundaries.

One of the geometries' resulting steady-state concentration profiles is plotted in Fig. 11. As can be seen in this figure, for small reaction rates, the concentration profile is close to linear, and the reaction rate dominates the behavior. However, diffusion becomes the more dominant factor when the reaction rate increases. At some point, the reaction happens so fast that the reactant cannot diffuse to the center of the system before being completely consumed. In those cases, a dead zone with zero effectiveness factor is formed at the center of the geometry.

Fig. 11 Concentration profiles in the bijel with particle volume fraction Ξ = 0.20 and particle contact angle θ_p = 90° for chemical reactions without advection. On the horizontal axis, we have the fractional distance to the center q; on the vertical axis, the mass fraction of the reactant ζ^r. The reaction rate ξ is displayed as a logarithmically varying color from blue (ξ = 10⁻⁷) to red (ξ = 10⁻¹).

It is possible to model the studied systems as a quasi-1D problem with a reaction-diffusion equation⁶¹

for the reactant. In this equation, x is used to index different cross sections. This equation can be rewritten in terms of the density ratio f and the Thiele modulus θ, describing the relation between diffusion and reaction rates in catalysts,⁶² aswithbeing the fractional distance from the center of the simulation domain in the x-direction. Additionally, the density ratio at a cross-section is computed in terms of the densities aswhere a difference due to the equilibrium concentration is taken into account. Solving eqn (29) with boundary conditions and f(−1) = f(1) = 1 gives

If we use eqn (32) to do least squares fits of the discrete density profiles produced by the reactive flow simulations, Thiele moduli can be approximated, giving us an indication of the influence of the rate of diffusion compared to the reaction rate. Since eqn (32) with fitted Thiele moduli gives us predictions of the concentration profiles, it also enables us to compute predictions for the effectiveness factors. These can then be compared with the effectiveness factors computed directly from the simulations to check the one-dimensional model's validity for various simulation parameters.

Fig. 12 shows ratios of the effectiveness factor from bijel simulations and the effectiveness factor obtained from simulations of their stochastic equivalents with the same simulation parameters. There is an improvement in the effectiveness factors of bijels with respect to stochastic structures for certain reaction rates and geometries. The predicted effectiveness factor ratios based on the computed Thiele moduli are also shown. There is a clear deviation between the values computed directly from the simulations and the theoretical predictions for large reaction rates. The 1D model description of the system loses its validity here since in these cases, the entirety of the reaction takes place in the first few cross-sections away from the recoloring layer. Hence, the structure of these first few layers has a major influence on the final steady-state profile of the simulations.

Fig. 12 The simulated effectiveness factor for bijels with Ξ= 0.20 divided by the simulated effectiveness factor of their equivalent stochastic geometry (closed symbols). The reaction rates are varied, and the simulations are in the diffusion-dominated flow regime. The symbol color denotes the equilibrium particle contact angle; varying from blue (θ_p = 90°) to red (θ_p ≈ 117°). Predictions for the effectiveness factor ratios based on the approximated Thiele moduli (open symbols) show a clear deviation from the simulated values for large reaction rates.

6.2. Advection-dominated flows

Advection-dominated flows are modeled in a similar manner as diffusion-dominated flows. However, we apply a strong acceleration a_x = 10⁻⁴ on the fluid to diminish diffusion-related effects. Additionally, we add 128 layers of void sites before and after the geometry in the x direction. This gives the flow room to equilibrate before entering the recoloring layer at the end of the domain, and ensures these simulations are in the high Péclet number regime.

Fig. 13 displays the concentration profiles of advection-dominated reactive flow simulations in one of the geometries. The linear concentration profiles in the low reaction rate regime and the strongly non-linear behavior in the high reaction rate regime are visible.

Fig. 13 The steady-state concentration profiles for advection-dominated flow in the bijel-derived geometry with particle volume fraction Ξ = 0.20 and particle contact angle θ_p = 90°. The mass fraction of the reactant ζ^r is plotted against the fractional distance to the center q. The reaction rate ξ is varied logarithmically from blue (ξ = 10⁻⁶) to red (ξ = 10⁻²). The separation between the inlet layer, geometry, and outlet layer is depicted with black vertical lines.

Employing the same methodologies for the diffusion-dominated flows, it is also possible to plot the ratio of the effectiveness factors for these simulations. Fig. 14 shows these ratios for the simulations with bijel-based geometries with Ξ = 0.20 and their stochastic equivalents. For these bijels, an improvement in catalyst effectiveness of approximately 100% can be realized. As can also be seen from this figure, a clear geometry-dependent optimum for the reaction rate exists. For these particular simulations, the ratio between the timescales for the advective mass transfer and the reaction rate is close to unity and, hence, optimal. When the particle volume fraction of the bijels is also varied, the ratio between the bijels and the stochastic geometries can be further improved. Fig. 15 shows the same ratios for the geometries with θ_p = 90° but different particle volume fractions. In this case, the simulations show that for a geometry with Ξ = 0.15, an improvement in the effectiveness factor of more than 170% can be realized.

Fig. 14 The effectiveness factor ratio α_bijel/α_stochastic for advection-dominated flow in the bijels with particle volume fraction Ξ = 0.20. The colors denote the equilibrium particle contact angle θ_p which is linearly varied from blue (θ_p = 90°) to red (θ_p ≈ 117°).

Fig. 15 The effectiveness factor ratio α_bijel/α_stochastic for advection-dominated flow in the geometries with varying particle volume fraction and a particle contact angle θ_p = 90°. A considerable improvement of the bijel over the stochastic structures is realized with the Ξ = 0.15 bijel, achieving an effectiveness factor improvement of more than 2.7 times.

7. Conclusion

This work presents an exploratory numerical study of the applicability of bijels as a novel type of catalyst support structure. The lattice Boltzmann method was used to generate bijel-derived porous structures over a range of pore sizes, porosities, surface areas, and permeabilities.

It was shown that the properties of bijel-based porous media are more beneficial than those of their stochastic equivalents (with the same porosity and surface area) for use in chemical reactors. This was further demonstrated by miscible reactive flow simulations, showing that optimization of the structures allows for an almost threefold increase in chemical reactor effectiveness. This is a significantly larger improvement over what has been demonstrated experimentally.

Experimental challenges related to the reproducibility and scalability of bijel-based support generation remain. Further, ensuring mechanical stability and structural heterogeneity for large porous structures is still challenging. However, several possible pathways for manufacturing stable bijel-based porous structures across a wide range of parameters are elaborated upon in the literature. To this end, we are convinced that it is possible to experimentally realize and study bijel-based porous reactor morphologies with the optimized particle contact angles and particle volume fractions that we report in this work.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data that support the findings of this study are openly available on Zenodo: https://doi.org/10.5281/zenodo.18861001.

Acknowledgements

We acknowledge the Helmholtz Association of German Research Centers (HGF) and the Federal Ministry of Education and Research (BMBF), Germany, for providing the necessary funding via the Innovation Pool project “Solar H2: Highly Pure and Compressed”. Additionally, we thank the Gauss Centre for Supercomputing e.V. (https://www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC).

References

1. G. Rothenberg, Catalysis: concepts and green applications, Wiley-VCH, 2008.
2. G. Centi, P. Ciambelli, S. Perathoner and P. Russo, Catal. Today, 2002, 75, 3–15.
3. M. Morbidelli, A. Gavriilidis and A. Varma, Catalyst design: optimal distribution of catalyst in pellets, reactors, and membranes, Cambridge University Press, 2001.
4. J. W. Cahn, J. Chem. Phys., 1965, 42, 93–99.
5. P. S. Clegg, Bijels: Bicontinuous Particle-stabilized Emulsions, Royal Society of Chemistry, 2020.
6. S. U. Pickering, J. Chem. Soc., Trans., 1907, 91, 2001–2021.
7. K. Stratford, R. Adhikari, I. Pagonabarraga, J. Desplat and M. E. Cates, Science, 2005, 309, 2198–2201.
8. E. M. Herzig, K. A. White, A. B. Schofield, W. C. K. Poon and P. S. Clegg, Nat. Mater., 2007, 6, 966–971.
9. M. E. Cates and P. S. Clegg, Soft Matter, 2008, 4, 2132.
10. K. A. White, A. B. Schofield, P. Wormald, J. W. Tavacoli, B. P. Binks and P. S. Clegg, J. Colloid Interface Sci., 2011, 359, 126–135.
11. J. W. Tavacoli, J. H. J. Thijssen, A. B. Schofield and P. S. Clegg, Adv. Funct. Mater., 2011, 21, 2020–2027.
12. J. A. Witt, D. R. Mumm and A. Mohraz, Soft Matter, 2013, 9, 6773.
13. H. Firoozmand and D. Rousseau, Food Hydrocolloids, 2015, 48, 208–212.
14. N. Hijnen, D. Cai and P. S. Clegg, Soft Matter, 2015, 11, 4351–4355.
15. M. F. Haase, N. Sharifi-Mood, D. Lee and K. J. Stebe, ACS Nano, 2016, 10, 6338–6344.
16. D. Cai, P. S. Clegg, T. Li, K. A. Rumble and J. W. Tavacoli, Soft Matter, 2017, 13, 4824–4829.
17. C. Huang, J. Forth, W. Wang, K. Hong, G. S. Smith, B. A. Helms and T. P. Russell, Nat. Nanotechnol., 2017, 12, 1060–1063.
18. G. Di Vitantonio, T. Wang, K. J. Stebe and D. Lee, Appl. Phys. Rev., 2021, 8, 021323.
19. F. Pizzetti, A. Rossetti, A. Marchetti, F. Castiglione, V. Vanoli, E. Coste, V. Veneruso, P. Veglianese, A. Sacchetti, A. Cingolani and F. Rossi, Adv. Mater. Interfaces, 2021, 8, 2100991.
20. H. Ching, T. J. Thorson, B. Paul and A. Mohraz, Mater. Adv., 2021, 2, 5067–5075.
21. A. J. Sprockel, M. A. Khan, M. De Ruiter, M. T. Alting, K. A. Macmillan and M. F. Haase, Colloids Surf., A, 2023, 666, 131306.
22. S. Amirfattahi, H. Honaryar and Z. Niroobakhsh, Small Struct., 2024, 5, 2300187.
23. F. Jansen and J. Harting, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 046707.
24. S. Frijters, F. Günther and J. Harting, Soft Matter, 2012, 8, 6542.
25. E. Kim, K. Stratford and M. E. Cates, Langmuir, 2010, 26, 7928–7936.
26. F. Günther, F. Janoschek, S. Frijters and J. Harting, Comput. Fluids, 2013, 80, 184–189.
27. J. M. Carmack and P. C. Millett, J. Chem. Phys., 2015, 143, 154701.
28. J. M. Carmack and P. C. Millett, Soft Matter, 2018, 14, 4344–4354.
29. F. Bonaccorso, S. Succi, M. Lauricella, A. Montessori, A. Tiribocchi and K. H. Luo, AIP Adv., 2020, 10, 095304.
30. J. M. Steenhoff and M. F. Haase, Phys. Chem. Chem. Phys., 2025, 27, 5117–5130.
31. C. Gu and L. Botto, Comput. Phys. Commun., 2020, 256, 107447.
32. M. N. Lee and A. Mohraz, Adv. Mater., 2010, 22, 4836–4841.
33. M. Reeves, K. Stratford and J. H. J. Thijssen, Soft Matter, 2016, 12, 4082–4092.
34. K. M. McDevitt, T. J. Thorson, E. L. Botvinick, D. R. Mumm and A. Mohraz, Materialia, 2019, 7, 100393.
35. J. A. Witt, D. R. Mumm and A. Mohraz, J. Mater. Chem. A, 2016, 4, 1000–1007.
36. S. J. Gross, K. M. McDevitt, D. R. Mumm and A. Mohraz, ACS Appl. Mater. Interfaces, 2021, 13, 8528–8537.
37. N. Allasia, O. Nevskyi, M. Marelli, I. Plazl, J. Albertazzi, V. Busini, F. Castiglione, F. Rossi and G. Vilé, Chem. Eng. J., 2024, 499, 155885.
38. R. Benzi, S. Succi and M. Vergassola, Phys. Rep., 1992, 222, 145–197.
39. S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond, Clarendon Press, Oxford University Press, 2001.
40. J. Harting, J. Chin, M. Venturoli and P. V. Coveney, Philos. Trans. R. Soc., A, 2005, 363, 1895–1915.
41. X. Shan and H. Chen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 47, 1815–1819.
42. X. Shan and G. Doolen, J. Stat. Phys., 1995, 81, 379–393.
43. A. J. C. Ladd, Phys. Rev. Lett., 1993, 70, 1339–1342.
44. C. K. Aidun and Y. Lu, J. Stat. Phys., 1995, 81, 49–61.
45. A. J. C. Ladd and R. Verberg, J. Stat. Phys., 2001, 104, 1191–1251.
46. A. J. C. Ladd, Phys. Fluids, 1997, 9, 491–499.
47. H. Hertz, J. Reine Angew. Math., 1882, 1882, 156–171.
48. Q. Xie, G. B. Davies, F. Günther and J. Harting, Soft Matter, 2015, 11, 3581–3588.
49. Z. Guo, C. Zheng and B. Shi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 046308.
50. Colloidal particles at liquid interfaces, ed B. P. Binks and T. S. Horozov, Cambridge University Press, 2006.
51. G. B. Davies, T. Krüger, P. V. Coveney and J. Harting, J. Chem. Phys., 2014, 141, 154902.
52. S. Frijters, T. Krüger and J. Harting, Comput. Phys. Commun., 2015, 189, 92–98.
53. Q. Z. Ye, [1988 Proceedings] 9th International Conference on Pattern Recognition, 1988, pp. 495–499.
54. B. Münch and L. Holzer, J. Am. Ceram. Soc., 2008, 91, 4059–4067.
55. H. Jinnai, Y. Nishikawa and T. Hashimoto, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, R2554–R2557.
56. A. Duda, Z. Koza and M. Matyka, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 036319.
57. M. Matyka and Z. Koza, AIP Conf. Proc., 2012, 4, pp. 17–22.
58. A. Narváez and J. Harting, Advances in Applied Mathematics and Mechanics, 2010, 2, 685–700.
59. A. Narvá, T. Zauner, F. Raischel, R. Hilfer and J. Harting, J. Stat. Mech.: Theory Exp., 2010, 2010, P11026.
60. I. Gueven, S. Frijters, J. Harting, S. Luding and H. Steeb, Granular Matter, 2017, 19, 28.
61. R. L. York, K. M. Bratlie, L. R. Hile and L. K. Jang, Catal. Today, 2011, 160, 204–212.
62. E. W. Thiele, Ind. Eng. Chem., 1939, 31, 916–920.

---