Multiscale Continuum Model Research Articles

Modeling of Electric Double Layer to Analyze Reaction Kinetics

Electric double layers (EDLs) play a fundamental role in various electrochemical processes such as colloid dispersion, surface charging, and charge-transfer reactions. Increasingly, the role of EDLs on reaction kinetics is being studied (1). Despite the extensive history of EDL modeling, there remain challenges in accurate prediction of its structure. The characteristic length of EDL (nanometer scale) exceeds the grasp of typical ab-initio molecular-dynamic simulations in order to capture a whole picture of the statistical equilibrium. While continuum models offer a means to estimate the thermodynamic equilibrium of EDLs with substantially lower computational cost than molecular dynamics, state–of-the-art continuum models require parameter fitting(2) due to the lack of appropriate expressions for microscopic interactions.Herein, we propose a predictive multiscale continuum model of EDL that eliminates the empirical fitting and relies more on physical relations. This model computes the distribution of the electrostatic potential, electron density, and species’ concentrations by minimizing the total grand potential of the system. The grand potential is calculated from entropic energy, electrostatic energy, electron energy, solute-solute interactions, and electrode-solute interactions. These energy expressions include not only the terms from previous work(2), but also the microscopic interactions that are newly introduced in this work: polarization of hydration shells, electrostatic interaction in parallel plane toward the electrode, and ion-size-dependent entropy. The parameters that identify the electrode and electrolyte materials are obtained from independent experiments such as vacuum work function, bulk modulus, X-ray diffraction, and concentration-dependent dielectric constant. The developed model reproduces the trends in the experimental differential capacitance with multiple electrode and electrolyte materials (Ag (111)-NaF, Ag (111)-NaClO4, and Hg-NaF), which verifies the accuracy and predictiveness of the model in a system without significant specific electrode-electrolyte interactions. The difference in the capacitance between mercury and silver electrode was attributed to the difference in the electron stability in the metal that is parameterized from the material’s bulk modulus. Finally, extension of the model to systems with specific electrode-electrolyte interactions (e.g., Pt (111)-NaClO4) and the analysis of reaction kinetics will be introduced. Overall, the model framework and findings will provide insights into the EDL structures and guidance on how to tailor electrode-electrolyte interface for improved reaction kinetics. Acknowledgements: This work was partially supported by the Hydrogen and Fuel Cell Technologies Office (HFTO), Office of Energy Efficiency and Renewable Energy, U.S. Department of Energy (DOE) through the Million Mile Fuel Cell Truck (M2FCT) consortia, technology managers G. Kleen and D. Papageorgopoulos, under contract number DE-AC02- 05CH11231, as well as by a CRADA with Toyota Central R&D Labs., Inc. and by the Center for Ionomer-based Water Electrolysis (CIWE), a DOE Energy Earthshot Research Center. Reference: S. J. Shin, D. H. Kim, G. Bae, S. Ringe, H. Choi, H. K. Lim, C. H. Choi and H. Kim, Nat Commun, 13, 174 (2022).J. Huang, J Chem Theory Comput, 19, 1003 (2023).

Parameter-Fitting-Free Continuum Modeling of Electric Double Layer in Aqueous Electrolyte.

Electric double layers (EDLs) play fundamental roles in various electrochemical processes. Despite the extensive history of EDL modeling, there remain challenges in the accurate prediction of its structure without expensive computation. Herein, we propose a predictive multiscale continuum model of EDL that eliminates the need for parameter fitting. This model computes the distribution of the electrostatic potential, electron density, and species' concentrations by taking the extremum of the total grand potential of the system. The grand potential includes the microscopic interactions that are newly introduced in this work: polarization of solvation shells, electrostatic interaction in parallel plane toward the electrode, and ion-size-dependent entropy. The parameters that identify the electrode and electrolyte materials are obtained from independent experiments in the literature. The model reproduces the trends in the experimental differential capacitance with multiple electrode and nonadsorbing electrolyte materials (Ag(110) in NaF, Ag(110) in NaClO4, and Hg in NaF), which verifies the accuracy and predictiveness of the model and rationalizes the observed values to be due to changes in electron stability. However, our calculation on Pt(111) in KClO4 suggests the need for the incorporation of electrode/ion-specific interactions. Sensitivity analyses confirmed that effective ion radius, ion valence, the electrode's Wigner-Seitz radius, and the bulk modulus of the electrode are significant material properties that control the EDL structure. Overall, the model framework and findings provide insights into EDL structures and predictive capability at low computational cost.

Modeling the Environment-Dependent Kinetics of Oxygen Reduction Reaction – a Continuum Model for Electric Double Layer

For proton-exchange-membrane fuel cells (PEMFCs) to achieve broad commercialization, improved energy-conversion efficiency with minimal Pt-based electrocatalyst is required [1]. Because the sluggish rate of oxygen reduction reaction (ORR) limits the efficiency of PEMFCs [2], the efficiency improvement requires a better understanding of ORR kinetics and mechanism to design better catalyst. To understand the ORR mechanism, theoretical and experimental analyses have been conducted. While previous studies reasonably explained the catalyst-dependent activity on single crystal catalysts in 0.1 M perchloric acid solution [3, 4], the explicit effect of electrolyte and related microenvironments [5, 6] is not thoroughly understood. The change in the electrolyte alters the electric-double-layer (EDL) structure and thus the local microenvironment at the electrode/electrolyte interface. Thus, the structure of the EDL should be carefully analyzed to uncover the electrolyte-dependent reaction kinetics. In this talk, we propose a multiscale continuum model to predict the EDL structure and examine the effect of perchloric acid concentration on ORR activity on Pt (111). The model includes Density Potential Functional Theory (DPFT) for electron density and Modified Poisson Boltzmann equation for species’ density and electric potential [7]. Also, the interaction between adsorbents and electric field is taken into account by minimizing the grand potential. After model validation with experimentally measured double-layer capacity data as a function of applied potential and concentration, the effect of the perchloric acid concentration (0.02 M – 0.2 M) on ORR activity is analyzed and discussed. It is shown that the model reproduces the specific activity obtained in the experiments when assuming the oxygen adsorption is limiting the rate, which can be attributed to the large energetic barrier for solvent reorganization [8]. Then, extension of the model to PEMFC ionomer electrolytes will be introduced. Overall, the model framework and findings provide insights into the ORR mechanism and guidance on how to tailor catalyst materials for increased PEMFC performance. Acknowledgements: This work was partially supported by the Hydrogen and Fuel Cell Technologies Office (HFTO), Office of Energy Efficiency and Renewable Energy, U.S. Department of Energy (DOE) through the Million Mile Fuel Cell Truck (M2FCT) consortia, technology managers G. Kleen and D. Papageorgopoulos, under contract number DE-AC02- 05CH11231, as well as by a CRADA with Toyota Central R&D Labs., Inc.. Reference [1] D. A. Cullen, K. C. Neyerlin, R. K. Ahluwalia, R. Mukundan, K. L. More, R. L. Borup, A. Z. Weber, D. J. Myers and A. Kusoglu, Nature Energy, 6, 462 (2021).[2] I. E. L. Stephens, A. S. Bondarenko, U. Grønbjerg, J. Rossmeisl and I. Chorkendorff, Energy & Environmental Science, 5 (2012).[3] M. Shibata, M. Inaba, K. Shinozaki, K. Kodama and R. Jinnouchi, Journal of The Electrochemical Society, 169 (2022).[4] H. A. Hansen, V. Viswanathan and J. K. Nørskov, The Journal of Physical Chemistry C, 118, 6706 (2014).[5] K. Kodama, K. Motobayashi, A. Shinohara, N. Hasegawa, K. Kudo, R. Jinnouchi, M. Osawa and Y. Morimoto, ACS Catalysis, 8, 694 (2017).[6] M. Luo and M. T. M. Koper, Nature Catalysis, 5, 615 (2022).[7] J. Huang, J Chem Theory Comput, 19, 1003 (2023).[8] D. T. Limmer, A. P. Willard, P. Madden and D. Chandler, Proceedings of the National Academy of Sciences, 110, 4200 (2013).

A multiscale continuum model for the mechanics of hyperelastic composite reinforced with nanofibers

A multiscale continuum model for the mechanics of hyperelastic composite reinforced with nanofibers

Modeling the Environment-Dependent Kinetics of Oxygen Reduction Reaction – Effect of Relative Humidity

For proton exchange membrane fuel cells (PEMFCs) to achieve broad commercialization, the kinetics of oxygen reduction reaction (ORR) need to be better understood to improve the efficiency with minimized use of expensive Pt-based electrocatalyst(1). The increasing demands of PEMFC-powered heavy-duty vehicles make this issue critical, especially to meet the efficiency target (ultimately 72%)(2). Previous first principle-based theoretical analyses reasonably explained catalyst-dependent activity in ideal conditions (on single crystal catalysts in HClO4 solution under room temperature)(3). To design a practical catalyst, further analysis is required to bridge the gap between the ideal and actual PEMFC conditions (typically, covered with perfluorosulfonic acids under 60 to 90 °C, and 30 to 100% RH). In this talk, we examine the effect of relative humidity (RH) on ORR activity on a platinum surface. The effect of RH on ORR activity has been a controversial topic: an experimental study suggested lower RH increases the ORR activity(4), while another study showed the opposite trend(5). By using a multiscale continuum model, we validate our hypothesis that the controversial experimental results stem from the difference in the properties of the electrode/electrolyte interface. In the mathematical calculation, we assume that the polycrystalline platinum catalysts consist of multiple single crystalline surfaces(6) and that their activity is described by the summation of the activities on the single crystalline surfaces. Based on transition state theory, the reaction rate of each elementary step is calculated from the reaction free energy, reactant concentration, and activation energy with explicit accounting of the double-layer structure and interactions. In the presentation, the calculation results will be compared with experimental solid-state data using microelectrodes, and the mechanism of RH-dependent ORR activity will be discussed. Reference I. E. L. Stephens, A. S. Bondarenko, U. Grønbjerg, J. Rossmeisl and I. Chorkendorff, Energy & Environmental Science, 5 (2012).D. A. Cullen, K. C. Neyerlin, R. K. Ahluwalia, R. Mukundan, K. L. More, R. L. Borup, A. Z. Weber, D. J. Myers and A. Kusoglu, Nature Energy, 6, 462 (2021).H. A. Hansen, V. Viswanathan and J. K. Nørskov, The Journal of Physical Chemistry C, 118, 6706 (2014).M. Shibata, M. Inaba, K. Shinozaki, K. Kodama and R. Jinnouchi, Journal of The Electrochemical Society, 169, 044507 (2022).H. Xu, Y. Song, H. R. Kunz and J. M. Fenton, Journal of The Electrochemical Society, 152 (2005).G. A. Tritsaris, J. Greeley, J. Rossmeisl and J. K. Nørskov, Catalysis Letters, 141, 909 (2011).

A multiscale micromorphic model with strain rate relationship between MD simulations and macroscale experimental tests and dynamic heterogeneity for glassy polymers

This paper proposes a micromorphic simulation model with novel strain rate relationships between MD simulations and macroscale experiments at the ductile-brittle (DB) transition of glassy polymer and dynamic heterogeneity of deformation. Novelties of this paper are on the strain rate relationship between two different scale regions and investigations on the effects of molecular weight on the relationship. From molecular dynamics model having realistically entangled polymer networks and experiments, the strain rate relationships are derived and then realized through a two-scale multiscale continuum model based on the micromorphic theory. The investigation on the effects of molecular weights on the relationship yields that the relationship is maintained when the molecular weight is higher than the entanglement length. The relationship confirmed that the dynamic heterogeneity alleviates as the polymer enters into the plastic flow above the glass transition temperature and disappears above the melting temperature. These findings will be a basis for developing the constitutive relationship of glassy polymers through a multiscale approach that can incorporate the heterogeneous dynamics of molecules.

A four-compartment multiscale model of fluid and drug distribution in vascular tumours.

The subtle relationship between vascular network structure and mass transport is vital to predict and improve the efficacy of anticancer treatments. Here, mathematical homogenisation is used to derive a new multiscale continuum model of blood and chemotherapy transport in the vasculature and interstitium of a vascular tumour. This framework enables information at a range of vascular hierarchies to be fed into an effective description on the length scale of the tumour. The model behaviour is explored through a demonstrative case study of a simplified representation of a dorsal skinfold chamber, to examine the role of vascular network architecture in influencing fluid and drug perfusion on the length scale of the chamber. A single parameter, P, is identified that relates tumour‐scale fluid perfusion to the permeability and density of the capillary bed. By fixing the topological and physiological properties of the arteriole and venule networks, an optimal value for P is identified, which maximises tumour fluid transport and is thus hypothesised to benefit chemotherapy delivery. We calculate the values for P for eight explicit network structures; in each case, vascular intervention by either decreasing the permeability or increasing the density of the capillary network would increase fluid perfusion through the cancerous tissue. Chemotherapeutic strategies are compared and indicate that single injection is consistently more successful compared with constant perfusion, and the model predicts optimal timing of a second dose. These results highlight the potential of computational modelling to elucidate the link between vascular architecture and fluid, drug distribution in tumours.

A multiscale continuum model for inelastic behavior of woven composite

A multiscale continuum model to predict the anisotropic non-linear behavior and strength of 2D and 3D woven composite has been proposed. The proposed approach is based on representative volume element (RVE), which serves as a two-level approximation, first at the microscale via fibers enclosed with cohesive interphase embedded in a matrix and second at the mesoscale via tows interlaced with each other and reinforced inside a matrix. The constitutive property fields of the RVE’s at two scales are evaluated under periodic boundary conditions and quasi-static loading in six different directions. The mechanism’s such as built-in imperfections at the surface, residual stresses and glitches that are generally neglected in such models have been incorporated combinedly in the form of cohesive interphase at lower scale. The accuracy of the proposed approach is assessed by a comparison of the literature data and the predicted results. The fundamental failure mechanism at both the scales have been addressed using preeminent failure criterions. From the various techniques used, homogenization in conjunction with finite element analysis is a precise and consistent way to analyze the nonlinear mechanical response and prediction of failure in woven composites and their tows.

Quantifying Reaction and Rate Heterogeneity in Battery Electrodes in 3D through Operando X-ray Diffraction Computed Tomography.

In composite battery electrode architectures, local limitations in ionic and electronic transport can result in nonuniform energy storage reactions. Understanding such reaction heterogeneity is important to optimizing battery performance, including rate capability and mitigating degradation and failure. Here, we use spatially resolved X-ray diffraction computed tomography to map the reaction in a composite electrode based on the LiFePO4 active material as it undergoes charge and discharge. Accelerated reactions at the electrode faces in contact with either the separator or the current collector demonstrate that both ionic and electronic transport limit the reaction progress. The data quantify how nonuniformity of the electrode reaction leads to variability in the charge/discharge rate, both as a function of time and position within the electrode architecture. Importantly, this local variation in the reaction rate means that the maximum rate that individual cathode particles experience can be substantially higher than the average, control charge/discharge rate, by a factor of at least 2-5 times. This rate heterogeneity may accelerate rate-dependent degradation pathways in regions of the composite electrode experiencing faster-than-average reaction and has important implications for understanding and optimizing rate-dependent battery performance. Benchmarking multiscale continuum model parameters against the observed reaction heterogeneity permits extension of these models to other electrode geometries.

Tissue loading and microstructure regulate the deformation of embedded nerve fibres: predictions from single-scale and multiscale simulations

Excessive deformation of nerve fibres (axons) in the spinal facet capsular ligaments (FCLs) can be a cause of pain. The axons are embedded in the fibrous extracellular matrix (ECM) of FCLs, so understanding how local fibre organization and micromechanics modulate their mechanical behaviour is essential. We constructed a computational discrete-fibre model of an axon embedded in a collagen fibre network attached to the axon by distinct fibre-axon connections. This model was used to relate the axonal deformation to the fibre alignment and collagen volume concentration of the surrounding network during transverse, axial and shear deformations. Our results showed that fibre alignment affects axonal deformation only during transverse and axial loading, but higher collagen volume concentration results in larger overall axonal strains for all loading cases. Furthermore, axial loading leads to the largest stretch of axonal microtubules and induces the largest forces on axon's surface in most cases. Comparison between this model and a multiscale continuum model for a representative case showed that although both models predicted similar averaged axonal strains, strain was more heterogeneous in the discrete-fibre model.

A multi-scale continuum model of skeletal muscle mechanics predicting force enhancement based on actin-titin interaction.

Although recent research emphasises the possible role of titin in skeletal muscle force enhancement, this property is commonly ignored in current computational models. This work presents the first biophysically based continuum-mechanical model of skeletal muscle that considers, in addition to actin-myosin interactions, force enhancement based on actin-titin interactions. During activation, titin attaches to actin filaments, which results in a significant reduction in titin's free molecular spring length and therefore results in increased titin forces during a subsequent stretch. The mechanical behaviour of titin is included on the microscopic half-sarcomere level of a multi-scale chemo-electro-mechanical muscle model, which is based on the classic sliding-filament and cross-bridge theories. In addition to titin stress contributions in the muscle fibre direction, the continuum-mechanical constitutive relation accounts for geometrically motivated, titin-induced stresses acting in the muscle's cross-fibre directions. Representative simulations of active stretches under maximal and submaximal activation levels predict realistic magnitudes of force enhancement in fibre direction. For example, stretching the model by 20% from optimal length increased the isometric force at the target length by about 30%. Predicted titin-induced stresses in the muscle's cross-fibre directions are rather insignificant. Including the presented development in future continuum-mechanical models of muscle function in dynamic situations will lead to more accurate model predictions during and after lengthening contractions.

Gating of a mechanosensitive channel due to cellular flows

A multiscale continuum model is constructed for a mechanosensitive (MS) channel gated by tension in a lipid bilayer membrane under stresses due to fluid flows. We illustrate that for typical physiological conditions vesicle hydrodynamics driven by a fluid flow may render the membrane tension sufficiently large to gate a MS channel open. In particular, we focus on the dynamic opening/closing of a MS channel in a vesicle membrane under a planar shear flow and a pressure-driven flow across a constriction channel. Our modeling and numerical simulation results quantify the critical flow strength or flow channel geometry for intracellular transport through a MS channel. In particular, we determine the percentage of MS channels that are open or closed as a function of the relevant measure of flow strength. The modeling and simulation results imply that for fluid flows that are physiologically relevant and realizable in microfluidic configurations stress-induced intracellular transport across the lipid membrane can be achieved by the gating of reconstituted MS channels, which can be useful for designing drug delivery in medical therapy and understanding complicated mechanotransduction.

An Axisymmetric Boundary Element Model for Determination of Articular Cartilage Pericellular Matrix Properties In Situ via Inverse Analysis of Chondron Deformation

The pericellular matrix (PCM) is the narrow tissue region surrounding all chondrocytes in articular cartilage and, together, the chondrocyte(s) and surrounding PCM have been termed the chondron. Previous theoretical and experimental studies suggest that the structure and properties of the PCM significantly influence the biomechanical environment at the microscopic scale of the chondrocytes within cartilage. In the present study, an axisymmetric boundary element method (BEM) was developed for linear elastic domains with internal interfaces. The new BEM was employed in a multiscale continuum model to determine linear elastic properties of the PCM in situ, via inverse analysis of previously reported experimental data for the three-dimensional morphological changes of chondrons within a cartilage explant in equilibrium unconfined compression (Choi, et al., 2007, "Zonal Changes in the Three-Dimensional Morphology of the Chondron Under Compression: The Relationship Among Cellular, Pericellular, and Extracellular Deformation in Articular Cartilage," J. Biomech., 40, pp. 2596-2603). The microscale geometry of the chondron (cell and PCM) within the cartilage extracellular matrix (ECM) was represented as a three-zone equilibrated biphasic region comprised of an ellipsoidal chondrocyte with encapsulating PCM that was embedded within a spherical ECM subjected to boundary conditions for unconfined compression at its outer boundary. Accuracy of the three-zone BEM model was evaluated and compared with analytical finite element solutions. The model was then integrated with a nonlinear optimization technique (Nelder-Mead) to determine PCM elastic properties within the cartilage explant by solving an inverse problem associated with the in situ experimental data for chondron deformation. Depending on the assumed material properties of the ECM and the choice of cost function in the optimization, estimates of the PCM Young's modulus ranged from approximately 24 kPa to 59 kPa, consistent with previous measurements of PCM properties on extracted chondrons using micropipette aspiration. Taken together with previous experimental and theoretical studies of cell-matrix interactions in cartilage, these findings suggest an important role for the PCM in modulating the mechanical environment of the chondrocyte.

A multiscale continuum model of the grain-size dependence of the stress hysteresis in shape memory alloy polycrystals

In this paper, a multiscale continuum model is proposed to study the effect of grain size on the macroscopic dissipative response of shape memory alloy polycrystals during isothermal thermoelastic phase transition. In the simplest one-dimensional (1D) heterogeneous structural hierarchy, a series of non-local and non-convex double-well continuum elements are employed to model the micro-instability and the macroscopic stress hysteresis of the material under uniaxial quasi-static stretching. Three characteristic length scales (specimen size L, grain size l and intrinsic material length g) of a bulk polycrystal are imbedded in the 1D chain model and their important roles in the macroscopic dissipation are quantified. It is shown that the specific energy dissipation or the width of the stress hysteresis is governed by two non-dimensional ratios, N(= L/ l) and l ¯ = ( l / g ) . For a given specimen of size L, the hysteresis decreases rapidly at either very large or small values of l. In particular, it vanishes when the grain size is reduced to the nano-scale where the grain size and the interface thickness become comparable. The above results of the 1D model are reproduced in a two-dimensional (2D) non-local numerical experiment on the energy dissipation during multiple domain evolution in heterogeneous strips. The predictions of the two models agree well qualitatively with the recent experimental observations of the stress hysteresis in nano-grained superelastic NiTi polycrystals.

Upscaling pore‐scale reactive transport equations using a multiscale continuum formulation

Reactive transport equations are solved at the pore scale using the lattice Boltzmann (LB) method, and the results are upscaled using volume averaging over a representative elemental volume and are fit to a multiscale continuum model. The multiscale continuum model accounts for local concentration gradients within diffusion‐dominated matrix domains that are coupled to the primary continuum fluid. In general it is found that a multiscale continuum formulation is required to fit the upscaled pore‐scale results. However, it is also demonstrated that in some cases, multiscale processes may be represented by a single‐continuum model employing effective parameters that are not directly measurable. Provided that sufficient resolution of the pore‐scale geometry can be obtained, the pore‐scale model can be used to determine the most appropriate form of continuum formulation (single, dual, or multiple continua) that best fits the upscaled pore‐scale simulation and, simultaneously, to provide parameters needed for constitutive relations appearing in the multiscale continuum formulation. It is suggested that a multiscale continuum approach may help explain the observed discrepancy between laboratory and field‐derived reaction rates by explicitly representing distinct transport domains through separate interacting continua which could be responsible for the formation of preferential pathways. An example is presented on the basis of a multiscale, synthetic, structured porous medium describing transport of a tracer and linear reaction kinetics.

The application of the multiscale models for description of the dispersed composites

We studied composite materials reinforced by micro-particles. It is assumed that the following physical objects and appropriate local effects determine the specific properties of composites: cohesion fields and corresponding zones of internal interaction; adhesion interactions, which determine the peculiarity of interaction of contacting bodies. These effects are important and can define the mechanical properties of the mediums with developed interface surfaces. In this paper, we use the multiscale continuum model of solids to explain the specific properties of composite materials with thin structures associated with local interactions of special type between inclusions and a matrix.

Nanomechanical modeling of the nanostructures and dispersed composites

Nanoparticles, hyperthin films, nanotubes and composite materials obtained on the base of such nanostructures exhibit very attractive mechanical properties and are great interest to researchers from continuum mechanics. In this paper, we intend to develop the multiscale continuum model of solids to explain the uncommon properties of the thin structures and the composite materials with thin structures associated with special type local interactions between nanoparticles and matrix. We have used the variation approach assuming that the internal interactions are determined by general character of kinematical connections. This approach allows introducing the system of internal interactions of various types consequently, which correspond to various types of kinematic restrictions in the considered mediums. The model of Cosserat’s pseudocontinuum with both non-free deformations, which depend on the other generalized coordinates and free deformations and rotations (continuous field of defects) but with symmetric stress tensor was proposed. The particular model is considered, which according to the authors view may serve as base for investigations of the scale effects due to cohesion interactions. Description of the cohesion field near top of cracks and composite materials with inclusions are considered in this research.