Stress, Modulus, and Strain Energy in Thin-Film Composite Battery Electrodes During Ion Insertion and Extraction

2.1. Estimation of the Apparent Young’s Modulus of Composite Electrodes Using Hashin-Shtrikman Bounds and Gibson-Ashby Porosity Correction

The porous composite electrode is considered as a two-level heterogeneous material. The first level is a dense solid skeleton composed of active material, conductive additive, and binder, and the second level introduces the porosity into the solid skeleton. The elastic modulus of the dense solid skeleton is estimated using Hashin-Shtrikman (HS) bounds [26], and then, the effect of porosity is delt with a Gibson-Ashby-type relative density correction [27].

The overall procedure can be expressed as

$\{m_i,{\rho }_i,E_i,{\nu }_i\}\to {\varphi }_i\to K_i,G_i\to K_{\mathrm{H}\mathrm{S}}^{\pm },G_{\mathrm{H}\mathrm{S}}^{\pm }\to E_s^{\pm }\to E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\pm },$ (2-1)

where $m_i$, ${\rho }_i$, $E_i$, and ${\nu }_i$ are the mass, density, Young’s modulus, and Poisson’s ratio of phase $i$, respectively. The subscripts $i=a,c,b$ denote the active material, conductive additive, and binder, respectively. ${\varphi }_i$ is the volume fraction of each phase. $K_i$ and $G_i$ mean the bulk modulus and shear modulus of phase $i$. $K_{\mathrm{H}\mathrm{S}}^{\pm }$ and $G_{\mathrm{H}\mathrm{S}}^{\pm }$ indicate the Hashin-Shtrikman-type estimates for the bulk modulus and shear modulus of the dense solid skeleton, respectively, where $\pm$ means the upper and lower bounds. The quantities $E_s^{\pm }$ denote the upper and lower estimates of the Young’s modulus of the dense solid skeleton, while $E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\pm }$ denote the corresponding apparent modulus of the whole porous electrode.

The true volume of each solid constituent is calculated from its mass and density as

$V_i={\displaystyle \frac{m_i}{{\rho }_i}},i=a,c,b.$ (2-2)

The total volume of the dense solid skeleton phase is therefore

$V_s=\sum _i{V}_i={\displaystyle \frac{m_a}{{\rho }_a}}+{\displaystyle \frac{m_c}{{\rho }_c}}+{\displaystyle \frac{m_b}{{\rho }_b}}.$ (2-3)

The volume fraction of each phase within the dense solid skeleton, excluding the pore phase, is defined as

${\varphi }_i={\displaystyle \frac{V_i}{V_s}}={\displaystyle \frac{m_i/{\rho }_i}{\sum _j{m}_j/{\rho }_j}},i,j=a,c,b,$ (2-4)

and $\sum _i{\varphi }_i=1.$The Hashin-Shtrikman formulation is used to estimate the bounds of the modulus of the dense solid skeleton. The Hashin-Shtrikman bounds provide rigorous variational upper and lower bounds for the effective elastic moduli of statistically isotropic multiphase composites. Compared with the Voigt and Reuss estimates, the Hashin-Shtrikman bounds are generally much narrower and more suitable for estimating the elastic properties of heterogeneous solid electrodes in the absence of detailed microstructural information.

Assuming each solid constituent to be isotropic and linearly elastic, the elastic constants $E_i$, bulk modulus $K_i$ and shear modulus $G_i$, of phase $i$ are converted according to

$K_i={\displaystyle \frac{E_i}{3(1-2{\nu }_i)}},$ (2-5a)

$G_i={\displaystyle \frac{E_i}{2(1+{\nu }_i)}}.$ (2-5b)

And then we have

$E_i={\displaystyle \frac{9K_i{G}_i}{3K_i+G_i}}$. (2-6)

If the Poisson’s ratios ${\nu }_i$ are not experimentally available, reasonable estimates should be adopted. For example, typical values may be 0.2~0.3 for carbonaceous active materials and conductive carbon, and 0.35~0.45 for PVDF-type polymer binders.

For a multiphase isotropic composite, the Hashin-Shtrikman-type estimate for the bulk modulus can be written as ²⁶

$K_{\mathrm{H}\mathrm{S}}\left(K_0,G_0\right)=K_0+{\left[\sum _i{\displaystyle \frac{{\varphi }_i}{K_i+\frac{4}{3}{G}_0}}\right]}^{-1}-(K_0+{\displaystyle \frac{4}{3}}{G}_0),$ (2-7)

where $K_0$ and $G_0$ are the reference bulk and shear moduli. It is clear that $K_{\mathrm{H}\mathrm{S}}\left(K_0,G_0\right)$ is factually only related to $G_0$.

The corresponding expression for the shear modulus is ²⁶

$G_{\mathrm{H}\mathrm{S}}(K_0,G_0)=G_0+{\left[\sum _i{\displaystyle \frac{{\varphi }_i}{G_i+{\beta }_0}}\right]}^{-1}-(G_0+{\beta }_0),$ (2-8)

where

${\beta }_0={\displaystyle \frac{G_0(9K_0+8G_0)}{6(K_0+2G_0)}}$. (2-9)

It can be seen that $G_{\mathrm{H}\mathrm{S}}(K_0,G_0)$ is related to both $K_0$ and $G_0$ through ${\beta }_0$.

For the upper estimate, the reference phase is selected as the stiffest phase:

$K_0=K_{\mathrm{m}\mathrm{a}\mathrm{x}},G_0=G_{\mathrm{m}\mathrm{a}\mathrm{x}}.$ (2-10)

Accordingly, the upper Hashin-Shtrikman estimate for the bulk modulus is

$K_{\mathrm{H}\mathrm{S}}^{+}={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{K_i+\frac{4}{3}{G}_{\mathrm{m}\mathrm{a}\mathrm{x}}}}\right]}^{-1}-{\displaystyle \frac{4}{3}}{G}_{\mathrm{m}\mathrm{a}\mathrm{x}},$ (2-11)

and the upper estimate for the shear modulus is

$G_{\mathrm{H}\mathrm{S}}^{+}={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{G_i+\frac{G_{\mathrm{m}\mathrm{a}\mathrm{x}}(9K_{\mathrm{m}\mathrm{a}\mathrm{x}}+8G_{\mathrm{m}\mathrm{a}\mathrm{x}})}{6(K_{\mathrm{m}\mathrm{a}\mathrm{x}}+2G_{\mathrm{m}\mathrm{a}\mathrm{x}})}}}\right]}^{-1}-{\displaystyle \frac{G_{\mathrm{m}\mathrm{a}\mathrm{x}}(9K_{\mathrm{m}\mathrm{a}\mathrm{x}}+8G_{\mathrm{m}\mathrm{a}\mathrm{x}})}{6(K_{\mathrm{m}\mathrm{a}\mathrm{x}}+2G_{\mathrm{m}\mathrm{a}\mathrm{x}})}}.$ (2-12)

Similarly, for the lower estimate, the reference phase is selected as the most compliant phase:

$K_0=K_{\mathrm{m}\mathrm{i}\mathrm{n}},G_0=G_{\mathrm{m}\mathrm{i}\mathrm{n}}.$ (2-13)

The lower Hashin-Shtrikman estimate for the bulk modulus is therefore

$K_{\mathrm{H}\mathrm{S}}^{-}={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{K_i+\frac{4}{3}{G}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right]}^{-1}-{\displaystyle \frac{4}{3}}{G}_{\mathrm{m}\mathrm{i}\mathrm{n}},$ (2-14)

and the lower estimate for the shear modulus is

$G_{\mathrm{H}\mathrm{S}}^{-}={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{G_i+\frac{G_{\mathrm{m}\mathrm{i}\mathrm{n}}(9K_{\mathrm{m}\mathrm{i}\mathrm{n}}+8G_{\mathrm{m}\mathrm{i}\mathrm{n}})}{6(K_{\mathrm{m}\mathrm{i}\mathrm{n}}+2G_{\mathrm{m}\mathrm{i}\mathrm{n}})}}}\right]}^{-1}-{\displaystyle \frac{G_{\mathrm{m}\mathrm{i}\mathrm{n}}(9K_{\mathrm{m}\mathrm{i}\mathrm{n}}+8G_{\mathrm{m}\mathrm{i}\mathrm{n}})}{6(K_{\mathrm{m}\mathrm{i}\mathrm{n}}+2G_{\mathrm{m}\mathrm{i}\mathrm{n}})}}.$ (2-15)

The corresponding upper and lower estimates of the Young’s modulus of the dense solid skeleton are calculated by converting $K_{\mathrm{H}\mathrm{S}}^{\pm }$ and $G_{\mathrm{H}\mathrm{S}}^{\pm }$ back to Young’s modulus through Equation (2-6), and we have

$E_s^{+}={\displaystyle \frac{9K_{\mathrm{H}\mathrm{S}}^{+}{G}_{\mathrm{H}\mathrm{S}}^{+}}{3K_{\mathrm{H}\mathrm{S}}^{+}+G_{\mathrm{H}\mathrm{S}}^{+}}},$ (2-16a)

$E_s^{-}={\displaystyle \frac{9K_{\mathrm{H}\mathrm{S}}^{-}{G}_{\mathrm{H}\mathrm{S}}^{-}}{3K_{\mathrm{H}\mathrm{S}}^{-}+G_{\mathrm{H}\mathrm{S}}^{-}}}.$ (2-16b)

Thus, the Young’s modulus of the dense solid skeleton is bounded by

$E_s^{-}\le E_s\le E_s^{+}.$ (2-17)

If a single representative value is required, the arithmetic mean of the two bounds may be used:

$E_s^{\mathrm{a}\mathrm{v}\mathrm{g}}={\displaystyle \frac{E_s^{-}+E_s^{+}}{2}}.$ (2-18)

After determining the modulus of the dense solid skeleton, the effect of porosity is introduced using a Gibson-Ashby-type scaling relation. The porosity of the porous electrode is defined as the ratio between the pore volume and the total porous electrode volume:

$\phi ={\displaystyle \frac{V_{\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{e}}}{V_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}}}.$ (2-19)

If the total electrode volume $V_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}$ is known from the electrode area $A$ and thickness $L$,

$V_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}=AL,$ (2-20)

then the porosity is calculated as

$\phi =1-{\displaystyle \frac{V_s}{V_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}}},$ (2-21)

or explicitly,

$\phi =1-{\displaystyle \frac{\frac{m_a}{{\rho }_a}+\frac{m_c}{{\rho }_c}+\frac{m_b}{{\rho }_b}}{AL}}.$ (2-22)

Alternatively, if the electrode density ${\rho }_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}$ is known, the electrode volume can be written as

$V_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}={\displaystyle \frac{m_a+m_c+m_b}{{\rho }_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}}},$ (2-23)

and the porosity becomes

$\phi =1-{\displaystyle \frac{\sum _i{m}_i/{\rho }_i}{\sum _i{m}_i/{\rho }_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}}}.$ (2-24)

According to Gibson-Ashby scaling relation, for a porous solid, the apparent Young’s modulus is related to the relative density by ²⁷

${\displaystyle \frac{E_{\mathrm{a}\mathrm{p}\mathrm{p}}}{E_s}}=C{\left({\displaystyle \frac{{\rho }_{\mathrm{a}\mathrm{p}\mathrm{p}}}{{\rho }_s}}\right)}^n,$ (2-25)

where $E_{\mathrm{a}\mathrm{p}\mathrm{p}}$ is the apparent modulus of the porous electrode, ${\rho }_{\mathrm{a}\mathrm{p}\mathrm{p}}$ is the apparent density of the porous electrode, ${\rho }_s$ is the density of the dense solid skeleton, $C$ is a structural coefficient, and $n$ is a porosity-dependent scaling exponent.

Since the relative density of the porous electrode is

${\displaystyle \frac{{\rho }_{\mathrm{a}\mathrm{p}\mathrm{p}}}{{\rho }_s}}=1-\phi ,$ (2-26)

the Gibson-Ashby relation becomes

$E_{\mathrm{a}\mathrm{p}\mathrm{p}}=C{E}_s(1-\phi {)}^n.$ (2-27)

Using the Hashin-Shtrikman upper and lower estimates for $E_s$, the corresponding bounds of the apparent electrode modulus are

$E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{+}=C{E}_s^{+}(1-\phi {)}^n.$ (2-28)

$E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{-}=C{E}_s^{-}(1-\phi {)}^n,$ (2-29)

Therefore, the apparent Young’s modulus of the porous composite electrode is estimated as

$C{E}_s^{-}(1-\phi {)}^n\le E_{\mathrm{a}\mathrm{p}\mathrm{p}}\le C{E}_s^{+}(1-\phi {)}^n.$ (2-30)

The parameters, $C$ and $n$ depend on pore morphology, pore connectivity, calendering history, particle packing, and binder distribution. For quantitative prediction, $C$ and $n$ should be calibrated using measured moduli of electrodes with different porosities or compaction densities.

When no experimental calibration is available, a common first approximation is to set $C=1,$ and to choose $n=2\text{–}4.$ For an initial estimate, one may take $n=3.$ Then the single-value estimate of the apparent modulus can then be approximated by

$E_{\mathrm{a}\mathrm{p}\mathrm{p}}\approx {\displaystyle \frac{E_s^{-}+E_s^{+}}{2}}(1-\phi {)}^3.$ (2-31)

More generally, one may use

$E_{\mathrm{a}\mathrm{p}\mathrm{p}}\approx C{\displaystyle \frac{E_s^{-}+E_s^{+}}{2}}(1-\phi {)}^n.$ (2-32)

Combining the Hashin–Shtrikman bounds with the Gibson–Ashby porosity correction gives

$E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\pm }=C{\displaystyle \frac{9K_{\mathrm{H}\mathrm{S}}^{\pm }{G}_{\mathrm{H}\mathrm{S}}^{\pm }}{3K_{\mathrm{H}\mathrm{S}}^{\pm }+G_{\mathrm{H}\mathrm{S}}^{\pm }}}(1-\phi {)}^n.$ (2-33)

Here,

$K_{\mathrm{H}\mathrm{S}}^{\pm }={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{K_i+\frac{4}{3}{G}_0^{\pm }}}\right]}^{-1}-{\displaystyle \frac{4}{3}}{G}_0^{\pm },$ (2-34)

and

$G_{\mathrm{H}\mathrm{S}}^{\pm }={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{G_i+\frac{G_0^{\pm }(9K_0^{\pm }+8G_0^{\pm })}{6(K_0^{\pm }+2G_0^{\pm })}}}\right]}^{-1}-{\displaystyle \frac{G_0^{\pm }(9K_0^{\pm }+8G_0^{\pm })}{6(K_0^{\pm }+2G_0^{\pm })}}.$ (2-35)

For the upper estimate,

$K_0^{+}=K_{\mathrm{m}\mathrm{a}\mathrm{x}},G_0^{+}=G_{\mathrm{m}\mathrm{a}\mathrm{x}},$ (2-36)

whereas for the lower estimate,

$K_0^{-}=K_{\mathrm{m}\mathrm{i}\mathrm{n}},G_0^{-}=G_{\mathrm{m}\mathrm{i}\mathrm{n}}.$ (2-37)

The above modeling provides a closed-form framework for estimating the apparent elastic modulus of a porous composite electrode using constituent masses, densities, elastic moduli, Poisson’s ratios, and electrode porosity.

2.2. Inverse Estimation of the Active Particle’s Young’s Modulus from the Apparent Modulus of Porous Composite Electrodes

When the apparent modulus of the porous electrode and the moduli of the conductive agent and binder are known, the Young’s modulus of the active particles can be estimated inversely. This is important when the modulus of the active particles varies with ion concentration and direct measurement is difficult.

Because the electrode is a heterogeneous porous medium and its microstructure is generally not explicitly resolved, the inverse problem is generally not unique. Here, the active particle Young’s modulus is therefore identified as an admissible interval rather than a single deterministic value. The forward model in Section 2.1 combines Hashin–Shtrikman bounds for the dense solid skeleton with a Gibson–Ashby-type porosity correction. The inverse problem is formulated by finding the range of active-particle moduli for which the measured apparent electrode modulus lies within the predicted upper and lower bounds.

The proposed inverse framework can be summarized as

$E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}\to E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}\to E_s^{-}\left(E_a\right)\le E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}\le E_s^{+}\left(E_a\right)\to E_a\in [E_a^{\mathrm{m}\mathrm{i}\mathrm{n}},E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}],$ (2-38)

where $E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}$ is the experimentally measured apparent Young’s modulus of the porous electrode, $E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}$ is the equivalent dense-skeleton modulus inferred from the porosity correction, and $E_s^{-}\left(E_a\right)$ and $E_s^{+}\left(E_a\right)$ are the Hashin–Shtrikman lower and upper estimates of the dense solid skeleton modulus as functions of the unknown active particle modulus $E_a$. The apparent Young’s modulus of the porous electrode is described by

$E_{\mathrm{a}\mathrm{p}\mathrm{p}}=C{E}_s(1-\phi {)}^n.$ (2-39)

The active particle modulus $E_a$ is treated as the unknown parameter, while the properties of the conductive additive and binder are assumed to be known.

For the active material, the bulk modulus $K_a$ and shear modulus $G_a$ are functions of the unknown active particle Young’s modulus $E_a$:

$K_a\left(E_a\right)={\displaystyle \frac{E_a}{3(1-2{\nu }_a)}},$ (2-40)

$G_a\left(E_a\right)={\displaystyle \frac{E_a}{2(1+{\nu }_a)}}.$ (2-41)

For the conductive additive and binder, the elastic constants are known:

$K_c={\displaystyle \frac{E_c}{3(1-2{\nu }_c)}},G_c={\displaystyle \frac{E_c}{2(1+{\nu }_c)}},$ (2-42)

$K_b={\displaystyle \frac{E_b}{3(1-2{\nu }_b)}},G_b={\displaystyle \frac{E_b}{2(1+{\nu }_b)}}.$ (2-43)

Therefore, the set of phase moduli entering the Hashin–Shtrikman model is

$\left\{K_i\right\}=\left\{K_a\right(E_a),K_c,K_b\},$ (2-44a)

$\left\{G_i\right\}=\left\{G_a\right(E_a),G_c,G_b\},$ (2-44b)

The Hashin–Shtrikman bounds provide upper and lower estimates of the effective bulk and shear moduli of a statistically isotropic multiphase composite.

For a given reference phase with bulk modulus $K_0$ and shear modulus $G_0$, the Hashin–Shtrikman estimate of the bulk modulus is given as

$K_{\mathrm{H}\mathrm{S}}(K_0,G_0;E_a)={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{K_i\left(E_a\right)+\frac{4}{3}{G}_0}}\right]}^{-1}-{\displaystyle \frac{4}{3}}{G}_0$, (2-45)

and the Hashin–Shtrikman estimate of the shear modulus is

$G_{\mathrm{H}\mathrm{S}}\left(K_0,G_0;E_a\right)={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{G_i\left(E_a\right)+{\beta }_0}}\right]}^{-1}-{\beta }_0.$ (2-46)

For the upper estimate, the stiffest reference phase is used:

$K_0^{+}=K_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(E_a\right),$ (2-47)

$G_0^{+}=G_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(E_a\right),$ (2-48)

where

$K_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(E_a\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left[K_a\left(E_a\right),K_c,K_b\right],$ (2-49)

$G_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(E_a\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left[G_a\left(E_a\right),G_c,G_b\right].$ (2-50)

The upper Hashin–Shtrikman estimate of the bulk modulus is then given by

$K_{\mathrm{H}\mathrm{S}}^{+}\left(E_a\right)={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{K_i\left(E_a\right)+\frac{4}{3}{G}_0^{+}}}\right]}^{-1}-{\displaystyle \frac{4}{3}}{G}_0^{+},$ (2-51)

and the upper Hashin–Shtrikman estimate of the shear modulus is

$G_{\mathrm{H}\mathrm{S}}^{+}\left(E_a\right)={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{G_i\left(E_a\right)+\frac{G_0^{+}(9K_0^{+}+8G_0^{+})}{6(K_0^{+}+2G_0^{+})}}}\right]}^{-1}-{\displaystyle \frac{G_0^{+}(9K_0^{+}+8G_0^{+})}{6(K_0^{+}+2G_0^{+})}}.$ (2-52)

For the lower estimate, the most compliant reference phase is used:

$K_0^{-}=K_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(E_a\right),$ (2-53)

$G_0^{-}=G_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(E_a\right),$ (2-54)

where

$K_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(E_a\right)=\mathrm{m}\mathrm{i}\mathrm{n}\left[K_a\left(E_a\right),K_c,K_b\right],$ (2-55)

$G_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(E_a\right)=\mathrm{m}\mathrm{i}\mathrm{n}\left[G_a\left(E_a\right),G_c,G_b\right].$ (2-56)

The lower Hashin–Shtrikman estimate of the bulk modulus is given by

$K_{\mathrm{H}\mathrm{S}}^{-}\left(E_a\right)={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{K_i\left(E_a\right)+\frac{4}{3}{G}_0^{-}}}\right]}^{-1}-{\displaystyle \frac{4}{3}}{G}_0^{-}.$ (2-57)

and the lower Hashin–Shtrikman estimate of the shear modulus is

$G_{\mathrm{H}\mathrm{S}}^{-}\left(E_a\right)={\left[\sum _i{\displaystyle \frac{{\varphi }_i}{G_i\left(E_a\right)+\frac{G_0^{-}(9K_0^{-}+8G_0^{-})}{6(K_0^{-}+2G_0^{-})}}}\right]}^{-1}-{\displaystyle \frac{G_0^{-}(9K_0^{-}+8G_0^{-})}{6(K_0^{-}+2G_0^{-})}}.$ (2-58)

The corresponding upper and lower Young’s moduli of the dense solid skeleton are obtained from

$E_s^{+}\left(E_a\right)={\displaystyle \frac{9K_{\mathrm{H}\mathrm{S}}^{+}\left(E_a\right)G_{\mathrm{H}\mathrm{S}}^{+}\left(E_a\right)}{3K_{\mathrm{H}\mathrm{S}}^{+}\left(E_a\right)+G_{\mathrm{H}\mathrm{S}}^{+}\left(E_a\right)}},$ (2-59)

$E_s^{-}\left(E_a\right)={\displaystyle \frac{9K_{\mathrm{H}\mathrm{S}}^{-}\left(E_a\right)G_{\mathrm{H}\mathrm{S}}^{-}\left(E_a\right)}{3K_{\mathrm{H}\mathrm{S}}^{-}\left(E_a\right)+G_{\mathrm{H}\mathrm{S}}^{-}\left(E_a\right)}}.$ (2-60)

Therefore, for a given active particle modulus $E_a$, the dense solid skeleton modulus satisfies

$E_s^{-}\left(E_a\right)\le E_s\left(E_a\right)\le E_s^{+}\left(E_a\right).$ (2-61)

According to the Gibson–Ashby relation, the measured apparent modulus of the porous electrode follows

$E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}=C{E}_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}(1-\phi {)}^n.$ (2-62)

Thus, the dense skeleton modulus required to reproduce the measured apparent modulus is

$E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}={\displaystyle \frac{E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{C(1-\phi {)}^n}}.$ (2-63)

This quantity represents the effective modulus that the dense solid skeleton would need to have in order to produce the measured porous electrode modulus under the assumed Gibson–Ashby porosity correction.

Instead of imposing a single equality between the predicted and measured modulus, the inverse problem is formulated using the Hashin–Shtrikman bounds. For a physically admissible active particle modulus $E_a$, the target dense skeleton modulus must lie between the lower and upper Hashin–Shtrikman estimates, that is

$E_s^{-}\left(E_a\right)\le E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}\le E_s^{+}\left(E_a\right).$ (2-64)

Substituting

$E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}={\displaystyle \frac{E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{C(1-\phi {)}^n}},$ (2-65)

one obtains the inverse admissibility condition:

$E_s^{-}\left(E_a\right)\le {\displaystyle \frac{E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{C(1-\phi {)}^n}}\le E_s^{+}\left(E_a\right).$ (2-66)

Equivalently, in terms of the apparent modulus,

$C{E}_s^{-}\left(E_a\right)(1-\phi {)}^n\le E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}\le C{E}_s^{+}(E_a\left)\right(1-\phi {)}^n.$ (2-67)

Therefore, the admissible active particle modulus set is defined as

$\mathcal{A}=\left\{E_a>0:C{E}_s^{-}\left(E_a\right)(1-\phi {)}^n\le E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}\le C{E}_s^{+}(E_a\left)\right(1-\phi {)}^n\right\}.$ (2-68)

If the admissible set is connected, it can be written as

$\mathcal{A}=[E_a^{\mathrm{m}\mathrm{i}\mathrm{n}},E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}].$ (2-69)

This interval represents the range of active particle Young’s moduli that are consistent with the measured apparent electrode modulus under the assumed Hashin–Shtrikman–Gibson–Ashby framework.

In most practical electrode systems, the predicted skeleton moduli $E_s^{-}\left(E_a\right)$ and $E_s^{+}\left(E_a\right)$ increase monotonically with the active-particle modulus $E_a$. Under this monotonicity assumption, the bounds of the admissible active particle modulus interval can be obtained by solving two scalar nonlinear equations.

The lower bound $E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}$ is obtained from the upper Hashin–Shtrikman estimate:

$C{E}_s^{+}\left(E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)(1-\phi {)}^n=E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}.$ (2-70)

Equivalently,

$E_s^{+}\left(E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)=E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}.$ (2-71)

The upper bound $E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}$ is obtained from the lower Hashin–Shtrikman estimate:

$C{E}_s^{-}\left(E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}\right)(1-\phi {)}^n=E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}.$ (2-72)

Equivalently,

$E_s^{-}\left(E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}\right)=E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}.$ (2-73)

Thus,

$E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}:E_s^{+}\left(E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)={\displaystyle \frac{E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{C(1-\phi {)}^n}}$ (2-74)

and

$E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}:E_s^{-}\left(E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}\right)={\displaystyle \frac{E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{C(1-\phi {)}^n}}.$ (2-75)

Because

$E_s^{-}\left(E_a\right)\le E_s^{+}\left(E_a\right),$ (2-76)

the modulus required to match the experimental value is reached at a smaller $E_a$ when the upper Hashin–Shtrikman estimate is used, and at a larger $E_a$ when the lower estimate is used. Therefore,

$E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}\le E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}.$ (2-77)

The final inverse estimate is

$E_a\in [E_a^{\mathrm{m}\mathrm{i}\mathrm{n}},E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}]$ (2-78)

with

$E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}={\left(E_s^{+}{\mid }_{\left[E_a^L,E_a^U\right]}\right)}^{-1}\left(E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}\right)$ (2-79)

and

$E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}={\left(E_s^{-}{\mid }_{\left[E_a^L,E_a^U\right]}\right)}^{-1}\left(E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}\right)$, (2-80)

provided that the inverse functions ${\left(\cdot \right)}^{-1}$ exist in the considered modulus range $\left[E_a^L,E_a^U\right]$.

The two boundary equations can be solved numerically by defining residual functions. For the lower bound of the active particle modulus, the residual is

$R^{+}\left(E_a\right)=C{E}_s^{+}\left(E_a\right)(1-\phi {)}^n-E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}.$ (2-81)

The root of this equation gives

$R^{+}\left(E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)=0.$ (2-82)

For the upper bound of the active particle modulus, the residual is

$R^{-}\left(E_a\right)=C{E}_s^{-}\left(E_a\right)(1-\phi {)}^n-E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}.$ (2-83)

The root of this equation gives

$R^{-}\left(E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}\right)=0.$ (2-84)

Alternatively, using the target dense skeleton modulus,

${\tilde{R}}^{+}\left(E_a\right)=E_s^{+}\left(E_a\right)-E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}},$ (2-85)

${\tilde{R}}^{-}\left(E_a\right)=E_s^{-}\left(E_a\right)-E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}.$ (2-86)

The roots are

${\tilde{R}}^{+}\left(E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)=0,$ (2-87)

${\tilde{R}}^{-}\left(E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}\right)=0.$ (2-88)

The use of the target dense skeleton modulus is often numerically convenient because the porosity correction is applied only once. Solving Equations (2-82) and (2-84) or (2-87) and (2-88) gives the $E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}$ and $E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}$, respectively. The lower bound $E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}$ corresponds to the case where the solid skeleton is assumed to achieve the upper Hashin–Shtrikman estimate. In this case, a smaller active particle modulus is sufficient to reproduce the measured apparent electrode modulus. In contrast, the upper bound $E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}$ corresponds to the case where the solid skeleton follows the lower Hashin–Shtrikman estimate. Because the lower bound predicts a more compliant skeleton, a larger active-particle modulus is required to match the same measured apparent electrode modulus.

A feasible solution exists only if the target dense skeleton modulus lies within the range spanned by the Hashin–Shtrikman bounds over the prescribed search interval. Let the search interval be

$E_a\in [E_a^L,E_a^U].$ (2-89)

A solution for $E_a^{\mathrm{m}\mathrm{i}\mathrm{n}}$ exists if

$\left[E_s^{+}\left(E_a^L\right)-E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}\right]\left[E_s^{+}\left(E_a^U\right)-E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}\right]\le 0.$ (2-90)

Similarly, a solution for $E_a^{\mathrm{m}\mathrm{a}\mathrm{x}}$exists if

$\left[E_s^{-}\left(E_a^L\right)-E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}\right]\left[E_s^{-}\left(E_a^U\right)-E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}\right]\le 0.$ (2-91)

If either condition is not satisfied, the corresponding root does not exist within the assumed search interval.

Although the interval-based inverse formulation is more physically appropriate, one may also define a representative skeleton modulus as the arithmetic mean of the Hashin–Shtrikman bounds:

$E_s^{\mathrm{a}\mathrm{v}\mathrm{g}}\left(E_a\right)={\displaystyle \frac{E_s^{-}\left(E_a\right)+E_s^{+}\left(E_a\right)}{2}}.$ (2-92)

A single equivalent active particle modulus can then be obtained by solving

$C{E}_s^{\mathrm{a}\mathrm{v}\mathrm{g}}\left(E_a\right)(1-\phi {)}^n=E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}.$ (2-93)

This value should not be interpreted as the unique intrinsic modulus of the active material. Rather, it should be regarded as a model-calibrated effective modulus.

The inverse estimate is highly sensitive to the porosity correction. Since $E_s^{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}}={\displaystyle \frac{E_{\mathrm{a}\mathrm{p}\mathrm{p}}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{C(1-\phi {)}^n}},$ small changes in $C$, $n$, or $\phi$ may produce large variations in the inferred active-particle modulus. Consequently, if possible, $C$ and $n$ should be calibrated using experimental modulus data from electrodes with different porosities or compaction densities.

3.1. Mechanical Description of the Constrained Electrode Layer

We first consider a freestanding electrode. The ion insertion induces a stress-free chemical expansion of the electrode film, and ions are treated as uniformly distributed in the active electrode film. The corresponding free linear chemical strain is denoted by

${\mathsf{\epsilon }}^{ch}\left(c\right),$where $c$ represents the ion concentration or state of charge.

Since the active electrode film was freestanding, it will expand freely by${\mathsf{\epsilon }}^{ch}\left(c\right)$ in all directions under isotropic chemical expansion. However, in a practical electrode, the active layer is bonded to a current collector, and its in-plane expansion is mechanically constrained, which is the situation we considered in the following discussion.

We denote the actual in-plane strain of the active layer of the substrate constrained electrode under equi-biaxial assumption by

${\epsilon }_{\parallel }^{tot,e}\left(c\right).$For isotropic active electrode layer, we have

${\epsilon }_x^{tot,e}\left(c\right)={\epsilon }_y^{tot,e}\left(c\right)={\epsilon }_{\parallel }^{tot,e}\left(c\right).$The in-plane mechanical strain of the active layer is therefore given by

${\epsilon }_{\parallel }^{m,e}\left(c\right)={\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)$, (3-1)

where ${\epsilon }_{\parallel }^{m,e}\left(c\right)<0$ indicates that the active layer is under in-plane compression.

For a sufficiently thick and stiff current collector, the in-plane deformation of the active layer may be approximated as fully constrained:

${\epsilon }_{\parallel }^{tot,e}\left(c\right)\approx 0.$ (3-2)

In this rigid-current-collector limit,

${\epsilon }_{\parallel }^{m,e}\left(c\right)\approx {\epsilon }_{\parallel }^{ch}\left(c\right)$. (3-3)

3.2. Concentration-Dependent Biaxial Modulus of the Active Layer and In-Plane Stress

The active layer is assumed to be isotropic and linearly elastic at each ion concentration. Both Young’s modulus and Poisson’s ratio may depend on concentration, i.e.,

$E_e=E_e\left(c\right),$ (3-4)

${\nu }_e={\nu }_e\left(c\right).$ (3-5)

Here, the $E_e$ and ${\nu }_e$ are equivalent to the apparent modulus $E_{app}$ and Poisson’s ratio ${\nu }_{app}$ of the composite electrode. In the following, we use the subscript e to replace the subscript app for simplicity. The concentration-dependent biaxial modulus of the active layer is defined as

$Y_e\left(c\right)={\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}$. (3-6)

Under plane-stress conditions,

${\sigma }_z^e=0,$ (3-7)

and under equi-biaxial in-plane deformation,

${\sigma }_x^e={\sigma }_y^e={\sigma }_{\parallel }^e.$ (3-8)

The in-plane stress in the active layer is therefore

${\sigma }_{\parallel }^e\left(c\right)=Y_e\left(c\right)\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$ (3-9)

or explicitly,

${\sigma }_{\parallel }^e\left(c\right)={\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$. (3-10)

In the rigid current collector limit, ${\epsilon }_{\parallel }^{tot,e}\left(c\right)\approx 0,$ and hence

${\sigma }_{\parallel }^e\left(c\right)\approx -{\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right)$. (3-11)

The negative sign denotes compressive stress in the active layer during ion insertion.

3.3. Out-of-Plane Strain and Thickness Evolution

Although the active layer is under plane stress, the out-of-plane strain is generally nonzero. The elastic out-of-plane Poisson strain is given by

${\epsilon }_z^{ela}\left(c\right)=-{\displaystyle \frac{{\nu }_e\left(c\right)}{E_e\left(c\right)}}\left[{\sigma }_x^e\left(c\right)+{\sigma }_y^e\left(c\right)\right].$ (3-12)

Since ${\sigma }_x^e={\sigma }_y^e={\sigma }_{\parallel }^e,$ we have

${\epsilon }_z^{ela}\left(c\right)=-{\displaystyle \frac{2{\nu }_e\left(c\right)}{E_e\left(c\right)}}{\sigma }_{\parallel }^e\left(c\right).$ (3-13)

Substituting Equation (3-10) into Equation (3-13) gives

${\epsilon }_z^{ela}\left(c\right)=-{\displaystyle \frac{2{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$. (3-14)

The total out-of-plane strain of the constrained active layer contains two contributions, i.e., the free chemical expansion strain, ${\epsilon }_z^{ch}\left(c\right)$, and the elastic Poisson strain, ${\epsilon }_z^{ela}\left(c\right)$. Therefore,

${\epsilon }_z^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)={\epsilon }_z^{ch}\left(c\right)-{\displaystyle \frac{2{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$. (3-15)

The active electrode layer thickness is then written as

$L_e\left(c\right)=L_{e0}\left[1+{\epsilon }_z^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)\right]$. (3-16)

where$L_{e0}$ is the initial active electrode layer thickness.

Substituting Equation (3-15) into Equation (3-16), we have

$L_e\left(c\right)=L_{e0}\left\{1+{\epsilon }_z^{ch}\left(c\right)-{\displaystyle \frac{2{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]\right\}$. (3-17a)

or equivalently,

$L_e\left(c\right)=L_{e0}\left\{1+{\epsilon }_{\parallel }^{ch}\left(c\right)-{\displaystyle \frac{2{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]\right\}$. (3-17b)

Equation (3-17) shows that the thickness evolution of the electrode layer depends directly on

${\epsilon }_{\parallel }^{ch}\left(c\right),{\nu }_e\left(c\right),{\epsilon }_{\parallel }^{tot,e}\left(c\right).$ (3-18)

Since ${\epsilon }_{\parallel }^{tot,e}\left(c\right)$ may itself depend on $E_e\left(c\right)$ through force balance with the current collector, the thickness $L_e\left(c\right)$ can also depend indirectly on $E_e\left(c\right)$.

3.4. Rigid Current Collector Limit

If the current collector is sufficiently thick and stiff, the active layer in-plane strain is approximately zero, i.e.,

${\epsilon }_{\parallel }^{tot,e}\left(c\right)=0.$ (3-19)

Equation (3-15) then becomes

${\epsilon }_z^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)={\epsilon }_{\parallel }^{ch}\left(c\right)+{\displaystyle \frac{2{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right).$ (3-20)

Thus,

${\epsilon }_z^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)={\displaystyle \frac{1+{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right)$, (3-21)

and the active layer thickness becomes

$L_e\left(c\right)=L_{e0}\left[1+{\displaystyle \frac{1+{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$. (3-22)

In this rigid substrate approximation, the thickness variation is governed by ${\epsilon }_{\parallel }^{ch}\left(c\right)$ and ${\nu }_e\left(c\right)$. The modulus $E_e\left(c\right)$ does not explicitly enter Equation (3-22), because the in-plane strain has already been prescribed as zero. However, the stress still depends directly on $E_e\left(c\right)$:

${\sigma }_{\parallel }^e\left(c\right)=-{\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right)$. (3-23)

3.5. Finite-Stiffness Current Collector: Coupled Dependence of $h_e\left(c\right)$

on$E_e\left(c\right)$When the current collector is not treated as perfectly rigid, the active layer and the current collector deform together. Perfect bonding requires in-plane strain compatibility:

${\epsilon }_{\parallel }^{tot,e}\left(c\right)={\epsilon }_{\parallel }^{cu}\left(c\right)={\epsilon }_{\parallel }^{tot}\left(c\right).$ (3-24)

The current collector is assumed to have no lithiation-induced chemical strain. Its stress is therefore

${\sigma }_{\parallel }^{cu}\left(c\right)=Y_{\mathrm{c}\mathrm{u}}{\epsilon }_{\parallel }^{tot}\left(c\right)$, (3-25)

where

$Y_{\mathrm{c}\mathrm{u}}={\displaystyle \frac{E_{\mathrm{c}\mathrm{u}}}{1-{\nu }_{\mathrm{c}\mathrm{u}}}}$ (3-26)

is the biaxial modulus of current collector.

The active electrode layer stress is

${\sigma }_{\parallel }^e\left(c\right)=Y_e\left(c\right)\left[{\epsilon }_{\parallel }^{tot}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$. (3-27)

Neglecting bending and assuming no external in-plane force, force balance per unit width gives

${\sigma }_{\parallel }^e\left(c\right)L_e\left(c\right)+{\sigma }_{\parallel }^{cu}\left(c\right)h_{\mathrm{c}\mathrm{u}}=0$, (3-28)

where $h_{\mathrm{c}\mathrm{u}}$ is the initial thickness of the current collector and is treated as constant. Noting that the variation in $h_{\mathrm{c}\mathrm{u}}$ induced by the elastic Poisson strain in the current collector is neglected.

Substituting Equations (3-25) and (3-27) into Equation (3-28) yields

$Y_e\left(c\right)L_e\left(c\right)\left[{\epsilon }_{\parallel }^{tot}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]+Y_{\mathrm{c}\mathrm{u}}{h}_{\mathrm{c}\mathrm{u}}{\epsilon }_{\parallel }^{tot}\left(c\right)=0.$ (3-29)

Solving Equation (3-29) for ${\epsilon }_{\parallel }^{tot}\left(c\right)$ gives

${\epsilon }_{\parallel }^{tot}\left(c\right)={\displaystyle \frac{Y_e\left(c\right)L_e\left(c\right)}{Y_e\left(c\right)L_e\left(c\right)+Y_{\mathrm{c}\mathrm{u}}{h}_{\mathrm{c}\mathrm{u}}}}{\epsilon }_{\parallel }^{ch}\left(c\right).$ (3-30)

Equation (3-30) shows explicitly that the actual in-plane strain depends on the concentration-dependent modulus $E_e\left(c\right)$.

Since $L_e\left(c\right)$ is given by Equation (3-17b), Equations (3-17b) and (3-30) form a coupled system:

$\left\{\begin{array}{c}{L}_e\left(c\right)=L_{e0}\left\{1+{\epsilon }_{\parallel }^{ch}\left(c\right)-{\displaystyle \frac{2{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}\left[{\epsilon }_{\parallel }^{tot}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]\right\},(3-31)\\ {\epsilon }_{\parallel }^{tot}\left(c\right)={\displaystyle \frac{Y_e\left(c\right)L_e\left(c\right)}{Y_e\left(c\right)L_e\left(c\right)+Y_{\mathrm{c}\mathrm{u}}{h}_{\mathrm{c}\mathrm{u}}}}{\epsilon }_{\parallel }^{ch}\left(c\right).(3-32)\end{array}\right.$Therefore, in the finite-stiffness current collector case, the thickness $L_e\left(c\right)$ depends on $E_e\left(c\right)$ indirectly through ${\epsilon }_{\parallel }^{tot}\left(c\right)$.

For compactness, we define in-plane extensional stiffness per unit width of the current collector

$S_{\mathrm{c}\mathrm{u}}=Y_{\mathrm{c}\mathrm{u}}{h}_{\mathrm{c}\mathrm{u}}$, (3-33)

and Poisson coupling factor

$\lambda \left(c\right)={\displaystyle \frac{2{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}$. (3-34)

Then, Equation (3-31) can be written as

$L_e\left(c\right)=L_{e0}\left[1+{\epsilon }_{\parallel }^{ch}\left(c\right)+\lambda \left(c\right){\epsilon }_{\parallel }^{ch}\left(c\right)-\lambda \left(c\right){\epsilon }_{\parallel }^{tot}\left(c\right)\right].$ (3-35)

Substituting Equation (3-30) into Equation (3-35) gives an implicit equation for $L_e\left(c\right)$. This equation can be rearranged into a quadratic form:

$Y_e\left(c\right)L_e^2\left(c\right)+\left[S_{\mathrm{c}\mathrm{u}}-L_{e0}{Y}_e\left(c\right)\left[1+{\epsilon }_{\parallel }^{ch}\left(c\right)\right]\right]L_e\left(c\right)-L_{e0}{S}_{\mathrm{c}\mathrm{u}}\left[1+\left(1+\lambda \left(c\right)\right){\epsilon }_{\parallel }^{ch}\left(c\right)\right]=0$. (3-36)

The physically meaningful positive solution is

$L_e\left(c\right)={\displaystyle \frac{L_{e0}{Y}_e\left(c\right)\left[1+{\epsilon }_{\parallel }^{ch}\left(c\right)\right]-S_{\mathrm{c}\mathrm{u}}+\sqrt{{\left[S_{\mathrm{c}\mathrm{u}}-L_{e0}{Y}_e\left(c\right)\left[1+{\epsilon }_{\parallel }^{ch}\left(c\right)\right]\right]}^2+4Y_e\left(c\right)L_{e0}{S}_{\mathrm{c}\mathrm{u}}\left[1+\left(1+\lambda \left(c\right)\right){\epsilon }_{\parallel }^{ch}\left(c\right)\right]}}{2Y_e\left(c\right)}}$. (3-36)

Equation (3-36) explicitly shows that, when the current collector has finite stiffness, the active layer thickness $L_e\left(c\right)$depends on the concentration-dependent modulus $E_e\left(c\right)$ through $Y_e\left(c\right)$. As long as the $E_e\left(c\right)$ and ${\nu }_e\left(c\right)$ are explicitly known, the $L_e\left(c\right)$ can be obtained by Equation (3-36).

After obtaining $L_e\left(c\right)$, the in-plane strain is calculated from

${\epsilon }_{\parallel }^{tot}\left(c\right)={\displaystyle \frac{Y_e\left(c\right)L_e\left(c\right)}{Y_e\left(c\right)L_e\left(c\right)+S_{\mathrm{c}\mathrm{u}}}}{\epsilon }_{\parallel }^{ch}\left(c\right),$ (3-37)

and the in-plane stress in the active electrode layer follows from

${\sigma }_{\parallel }^e\left(c\right)=Y_e\left(c\right)\left[{\epsilon }_{\parallel }^{tot}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$. (3-38)

On the other hand, if $E_e\left(c\right)$ is unknown while ${\nu }_e\left(c\right)$ can be treated as known and constant, and if ${\epsilon }_{\parallel }^{tot}\left(c\right)$ and ${\epsilon }_{\parallel }^{ch}\left(c\right)$ are measured experimentally, then $L_e\left(c\right)$ and $E_e\left(c\right)$ can be calculated by solving Equations (3-31) and (3-32) simultaneously. This is important for the practical application of the model, as ${\epsilon }_{\parallel }^{tot}\left(c\right)$ and ${\epsilon }_{\parallel }^{ch}\left(c\right)$ can be experimentally measured through carefully designed DIC-based in situ electrochemical tests.

3.6. Elastic Strain Energy Density

The elastic strain energy density in the active electrode layer under plane-stress conditions is

$u_e\left(c\right)={\displaystyle \frac{1}{2}}\left[{\sigma }_x^e\left(c\right){\epsilon }_x^{m,e}\left(c\right)+{\sigma }_y^e\left(c\right){\epsilon }_y^{m,e}\left(c\right)\right]$. (3-39)

For equi-biaxial deformation, ${\sigma }_x^e={\sigma }_y^e={\sigma }_{\parallel }^e,$ ${\epsilon }_x^{m,e}={\epsilon }_y^{m,e}={\epsilon }_{\parallel }^{m,e}$, and therefore

$u_e\left(c\right)={\sigma }_{\parallel }^e\left(c\right){\epsilon }_{\parallel }^{m,e}\left(c\right)$. (3-40)

Using Equation (3-9) or (3-10), we have

$u_e\left(c\right)=Y_e\left(c\right){\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]}^2$, (3-41)

or

$u_e\left(c\right)={\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]}^2$. (3-42)

The out-of-plane strain does not appear explicitly in Equation (3-42) because ${\sigma }_z=0$. However, it will affect the total strain energy through the current active layer thickness $L_e\left(c\right)$.

3.7. Elastic Strain Energy of the Active Electrode Layer

The elastic strain energy stored in the active electrode layer can be obtained by integrating the strain energy density over the current volume of the active electrode layer $V_e\left(c\right)$:

$U_e\left(c\right)={\int }_{V_e\left(c\right)}{u}_e\left(c\right)dV$. (3-43)

Assuming uniform stress and strain in the active layer, the energy per unit projected electrode area $A$ is

${\displaystyle \frac{U_e\left(c\right)}{A}}=L_e\left(c\right)Y_e\left(c\right){\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]}^2$, (3-44)

or

${\displaystyle \frac{U_e\left(c\right)}{A}}=L_e\left(c\right){\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]}^2$. (3-45)

For the rigid current collector case,

${\epsilon }_{\parallel }^{tot,e}\left(c\right)=0,$ and $L_e\left(c\right)=L_{e0}\left[1+{\displaystyle \frac{1+{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$,

thus, we have

${\displaystyle \frac{U_e\left(c\right)}{A}}=L_{e0}\left[1+{\displaystyle \frac{1+{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right)\right]{\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{{\epsilon }_{\parallel }^{ch}\left(c\right)}^2$, (3-46a)

or

$U_e\left(c\right)=A{L}_{e0}\left[1+{\displaystyle \frac{1+{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right)\right]{\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{{\epsilon }_{\parallel }^{ch}\left(c\right)}^2$. (3-46b)

For the finite-stiffness current collector case, $L_e\left(c\right)$ should be obtained from Equation (3-36), and then substituted into Equation (3-45).

When the elastic energy stored in the current collector is also included, the total strain energy becomes

$U_{\mathrm{t}\mathrm{o}\mathrm{t}}\left(c\right)=A{L}_e\left(c\right)Y_e\left(c\right){\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]}^2+A{h}_{\mathrm{c}\mathrm{u}}{Y}_{\mathrm{c}\mathrm{u}}{{\epsilon }_{\parallel }^{tot,e}\left(c\right)}^2$. (3-47)

Equation (3-47) accounts for both the compressive elastic energy stored in the active layer and the tensile elastic energy stored in the current collector.

3.8. Mechanical Work Generated During Ion Insertion

The incremental mechanical work associated with the constrained chemical expansion is

$dw={\sigma }_x^e\left(c\right)d{\epsilon }_x^{m,e}\left(c\right)+{\sigma }_y^e\left(c\right)d{\epsilon }_y^{m,e}\left(c\right)$. (3-48)

Under equi-biaxial in-plane deformation, ${\sigma }_x^e={\sigma }_y^e={\sigma }_{\parallel }^e,$ ${\epsilon }_x^{ch}={\epsilon }_y^{ch}={\epsilon }_{\parallel }^{ch}$, and ${\epsilon }_{\parallel }^{m,e}={\epsilon }_{\parallel }^{tot,e}-{\epsilon }_{\parallel }^{ch}$, then we have

$dw=2{\sigma }_{\parallel }^e\left(c\right)d\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]=2{\sigma }_{\parallel }^e\left(c\right)d{\epsilon }_{\parallel }^{tot,e}\left(c\right)-2{\sigma }_{\parallel }^e\left(c\right)d{\epsilon }_{\parallel }^{ch}\left(c\right)$. (3-49a)

For rigid current collector conditions, i.e., ${\epsilon }_{\parallel }^{tot,e}\left(c\right)\equiv 0$, then we have

$dw=-2{\sigma }_{\parallel }^e\left(c\right)d{\epsilon }_{\parallel }^{ch}\left(c\right)$. (3-49b)

The incremental work per unit projected area is therefore given by

$d\left({\displaystyle \frac{W}{A}}\right)=-2L_e\left(c\right){\sigma }_{\parallel }^e\left(c\right)d{\epsilon }_{\parallel }^{ch}\left(c\right)$. (3-50)

Integrating from the initial state (e.g., $c=0$) to $c=c_f$, the total mechanical work per unit projected area is

${\displaystyle \frac{W\left(c_f\right)}{A}}=-2{\int }_0^{c_f}{L}_e\left(c\right){\sigma }_{\parallel }^e\left(c\right){\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}dc$. (3-51)

Substituting Equation (3-27) into Equation (3-51) gives

$W\left(c_f\right)=2A{\int }_0^{c_f}{L}_e\left(c\right)Y_e\left(c\right){\epsilon }_{\parallel }^{ch}\left(c\right){\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}dc$ (3-52)

or explicitly,

$W\left(c_f\right)=2A{\int }_0^{c_f}{L}_e\left(c\right){\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right){\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}dc$, (3-53)

with

$L_e\left(c\right)=L_{e0}\left[1+{\displaystyle \frac{1+{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$. (3-54)

Equation (3-53) is the recommended general expression for the mechanical work when both $E_e\left(c\right)$ and $L_e\left(c\right)$ are concentration dependent. This expression should be used when the current collector is much stiffer than the active electrode layer and the active layer thickness variation is non-negligible.

For a current collector with finite stiffness, ${\epsilon }_{\parallel }^{tot,e}\left(c\right)$ is generally not zero or constant, thus the mechanical work must be calculated from the general incremental relation

$dw=2{\sigma }_{\parallel }^e\left(c\right)d\left[{\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$, (3-55)

and

$dW=2A{L}_e\left(c\right){\sigma }_{\parallel }^e\left(c\right)d\left[{\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$, (3-56)

using ${\sigma }_{\parallel }^e\left(c\right)={\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}\left[{\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]$.

The total mechanical work is then given by

$W\left(c_f\right)=2A{\int }_0^{c_f}{L}_e\left(c\right)Y_e\left(c\right)\left[{\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]{\displaystyle \frac{d}{dc}}\left[{\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]dc$, (3-57)

or equivalently,

$W\left(c_f\right)=2A{\int }_0^{c_f}{L}_e\left(c\right)Y_e\left(c\right)\left[{\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]\left[{\displaystyle \frac{d{\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)}{dc}}−{\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}\right]dc$, (3-58)

with ${\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)={\displaystyle \frac{Y_e\left(c\right)L_e\left(c\right)}{Y_e\left(c\right)L_e\left(c\right)+S_{cu}}}{\epsilon }_{\parallel }^{ch}\left(c\right)$. (3-59)

Defining $\eta \left(c\right)={\displaystyle \frac{Y_e\left(c\right)L_e\left(c\right)}{Y_e\left(c\right)L_e\left(c\right)+S_{cu}}}$, (3-60)being a dimensionless in-plane strain accommodation factor characterizing the fraction of the chemical expansion strain that is accommodated by the in-plane total deformation of the electrode, then the in-plane total strain can be written compactly as

${\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)=\eta \left(c\right){\epsilon }_{\parallel }^{ch}\left(c\right)$. (3-61)

Accordingly, the elastic mismatch strain becomes

${\epsilon }_{\parallel }^{m,e}\left(c\right)={\epsilon }_{\parallel }^{tot,e}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)=(\eta \left(c\right)-1){\epsilon }_{\parallel }^{ch}\left(c\right)$, (3-62)

and its derivative with respect to concentration is

${\displaystyle \frac{d}{dc}}\left[{\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]={\displaystyle \frac{d{\epsilon }_{\parallel }^{\mathrm{t}\mathrm{o}\mathrm{t},\mathrm{e}}\left(c\right)}{dc}}-{\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}={\displaystyle \frac{d\eta \left(c\right)}{dc}}{\epsilon }_{\parallel }^{ch}\left(c\right)+\left(\eta \right(c)-1){\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}$. (3-63)

Substituting Equations (3-61) and (3-63) into the general work formula (3-58) yields the more explicit form

$W\left(c_f\right)=2A{\int }_0^{c_f}{L}_e\left(c\right)Y_e\left(c\right)\left[\right(\eta \left(c\right)-1\left){\epsilon }_{\parallel }^{ch}\right(c\left)\right]\left[{\displaystyle \frac{d\eta \left(c\right)}{dc}}{\epsilon }_{\parallel }^{ch}\left(c\right)+\left(\eta \right(c)-1){\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}\right]dc$, (3-64)

where $L_e\left(c\right)$ and $Y_e\left(c\right)$ can be obtained from the coupled equations (3-31) and (3-32)

$\left\{\begin{array}{c}{L}_e\left(c\right)=L_{e0}\left\{1+{\epsilon }_{\parallel }^{ch}\left(c\right)-{\displaystyle \frac{2{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}\left[{\epsilon }_{\parallel }^{tot}\left(c\right)-{\epsilon }_{\parallel }^{ch}\left(c\right)\right]\right\},\\ {\epsilon }_{\parallel }^{tot}\left(c\right)={\displaystyle \frac{Y_e\left(c\right)L_e\left(c\right)}{Y_e\left(c\right)L_e\left(c\right)+Y_{\mathrm{c}\mathrm{u}}{h}_{\mathrm{c}\mathrm{u}}}}{\epsilon }_{\parallel }^{ch}\left(c\right).\end{array}\right.$The rigid current collector case can be recovered automatically as a limiting case of the finite-stiffness formulation. In the limit of infinitely large current collector stiffness,

$S_{cu}\to \mathrm{\infty }$, (3-65)

we have

$\eta \left(c\right)\to 0,{\epsilon }_{\parallel }^{tot,e}\left(c\right)\to 0$ (3-66)

Substituting Equations (3-65) and (3-66) into (3-64) gives

$W\left(c_f\right)=2A{\int }_0^{c_f}{L}_e\left(c\right)Y_e\left(c\right){\epsilon }_{\parallel }^{ch}\left(c\right){\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}dc$, (3-67)

which is exactly the same as the previously derived expression for the rigid current collector case Equation (3-53).

3.9. Difference Between Accumulated Mechanical Work and Final Elastic Strain Energy

When $E_e$, ${\nu }_e$, and $L_e$ are independent of concentration, the accumulated mechanical work may reduce to the final elastic strain energy. For example, for a rigid current collector, Equation (3-54)

$W\left(c_f\right)=2A{\int }_0^{c_f}{L}_e\left(c\right){\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right){\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}dc$will reduce to Equation (3-46b)

$U_e\left(c\right)=A{L}_{e0}\left[1+{\displaystyle \frac{1+{\nu }_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right)\right]{\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{{\epsilon }_{\parallel }^{ch}\left(c\right)}^2$.

However, when $E_e$, ${\nu }_e$, and $L_e$ are dependent of concentration

$E_e=E_e\left(c\right),$ ${\nu }_e={\nu }_e\left(c\right),$ and $L_e=L_e\left(c\right),$the process work and final stored elastic energy are generally different.

The accumulated work is a path integral,

$W\left(c_f\right)=2A{\int }_0^{c_f}{L}_e\left(c\right){\displaystyle \frac{E_e\left(c\right)}{1-{\nu }_e\left(c\right)}}{\epsilon }_{\parallel }^{ch}\left(c\right){\displaystyle \frac{d{\epsilon }_{\parallel }^{ch}\left(c\right)}{dc}}dc$,

whereas the final elastic energy is a state function,

$U_e\left(c_f\right)=A{L}_e\left(c_f\right){\displaystyle \frac{E_e\left(c_f\right)}{1-{\nu }_e\left(c_f\right)}}{\left[{\epsilon }_{\parallel }^{tot,e}\left(c_f\right)-{\epsilon }_{\parallel }^{ch}\left(c_f\right)\right]}^2$.

The difference arises because the elastic modulus, Poisson’s ratio, thickness, and constraint state may all evolve during ion insertion.

5.1. Physical System

We consider a planar film electrode of thickness $L$, perfectly bonded to a current collector/substrate of thickness $h$. A Cartesian coordinate $z$ is introduced along the thickness direction, with $z=0$ at the electrode/substrate interface, $z=L$ at the electrochemically active free surface of the electrode, and $-h\le z\le 0$ in the current collector. Since the in-plane dimensions are much larger than the thickness, edge effects are neglected and the problem is treated as one-dimensional in the thickness direction. The film electrode is assumed to deform within the regime of small strains and linear elasticity. For simplicity, the electrode material is taken as isotropic, with isotropic diffusion and isotropic chemical swelling. Meanwhile, the chemical diffusion coefficient $D$ and the electrode thickness $L$ are assumed to be constant and independent of concentration.

Because the electrode is bonded to a substrate, free in-plane expansion or contraction of the film is suppressed during ion insertion and extraction. As a result, concentration changes generate a mismatch between the natural chemical swelling strain of the electrode and the actual kinematically admissible deformation, which gives rise to in-plane stress. The top surface of the film is traction-free, so that the out-of-plane normal stress is negligible compared with the in-plane stress. Under these assumptions, the electrode develops an equi-biaxial in-plane stress state,

${\sigma }_x^e={\sigma }_y^e\equiv {\sigma }_{\parallel }^e(z,t),{\sigma }_z^e\approx 0.$ (5-1)

The ion concentration $c(z,t)$ satisfies the one-dimensional Fickian diffusion equation in the active electrode [29],

${\displaystyle \frac{\partial c}{\partial t}}=D{\displaystyle \frac{{\partial }^{2}c}{\partial z^2}},0<z<L,t>0,$ (5-2)

where $c(z,t)$ is the ion concentration and $D$ is the chemical diffusion coefficient.

In the present formulation, ion insertion and ion extraction are analyzed independently. For insertion experiments, all datasets are assumed to start from the same initial concentration $c_0^{\mathrm{i}n}$. For extraction experiments, all datasets are assumed to start from the same initial concentration $c_0^{\mathrm{e}x}$. These two initial concentrations may be different, depending on the experimental protocol.

For insertion,

$c(z,0)=c_0^{\mathrm{i}n},0\le z\le L.$ (5-3a)

For extraction,

$c(z,0)=c_0^{\mathrm{e}x},0\le z\le L.$ (5-3b)

Only the outer surface $z=L$ is assumed to be electrochemically active, while the substrate interface $z=0$ is impermeable. Accordingly,

${\displaystyle \frac{\partial c}{\partial z}}(0,t)=0.$ (5-4)

At the active surface, the imposed current density prescribes the ion molar flux. Let $i_0$ denote the applied current density, $\mathbb{n}$ the number of electrons transferred per ion, and $F$ Faraday’s constant. The interfacial molar flux is

$J_0={\displaystyle \frac{i_0}{\mathbb{n}F}}.$ (5-5)

In present case, we use lithium as an example, and thus $\mathbb{n}=1$. The ion insertion means lithiation and the ion extraction means delithiation.

The sign convention is chosen such that $J_0>0$ denotes inward lithium flux into the electrode during charging/lithiation. The electrochemical boundary condition at $z=L$ is therefore

$D{\displaystyle \frac{\partial c}{\partial z}}\left(L,t\right)=J_0.$ (5-6)

For lithiation,

$J_0>0,$ (5-7a)

whereas for delithiation,

$J_0<0.$ (5-7b)

Side reactions, interfacial kinetic limitations, concentration losses due to parasitic processes, and electrode thickness variation due to lithium insertion/extraction are neglected. Thus, the applied current is assumed to be fully converted into lithium insertion/extraction at the active surface. If both surfaces were electrochemically active, the second boundary condition would have to be modified accordingly; that case is not considered here.

The thickness-averaged concentration is defined as

$\bar{c}\left(t\right)={\displaystyle \frac{1}{L}}{\int }_0^{L}c(z,t)dz.$ (5-8)

Integrating the diffusion equation over the film thickness gives

${\displaystyle \frac{d\bar{c}}{dt}}={\displaystyle \frac{J_0}{L}}.$ (5-9)

Therefore, for a constant galvanostatic flux,

$\bar{c}\left(t\right)=c_0+{\displaystyle \frac{J_{0}t}{L}},$ (5-10)

where $c_0=c_0^{\mathrm{i}n}$ for lithiation and $c_0=c_0^{\mathrm{e}x}$ for delithiation.

It should be emphasized that Eq. (5-10) gives only the thickness-averaged concentration. Because lithium concentration is generally nonuniform through the electrode thickness, the local concentration $c(z,t)$ may be larger or smaller than the average value, especially at high current density or small diffusivity. Therefore, the concentration interval over which the local modulus function $Y_e\left(c\right)$ is identified should be determined from the local concentration field rather than from the averaged concentration alone.

Solving the diffusion equation under the above initial and boundary conditions gives

$c\left(z,t\right)=c_0+{\displaystyle \frac{J_{0}t}{L}}+{\displaystyle \frac{J_0{z}^2}{2DL}}-{\displaystyle \frac{J_{0}L}{6D}}-{\displaystyle \frac{2J_{0}L}{D{\pi }^2}}{\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \frac{{(-1)}^n}{n^2}}{e}^{-{\displaystyle \frac{n^2{\pi }^{2}Dt}{L^2}}}\cos \left({\displaystyle \frac{n\pi z}{L}}\right)$ (5-11)

Here $c_0$ should be interpreted as the initial concentration of the independently analyzed experiment group. Specifically, $c_0=c_0^{\mathrm{i}n}$ for lithiation datasets and $c_0=c_0^{\mathrm{e}x}$ for delithiation datasets.

For subsequent inverse identification, the concentration interval sampled by the experiments must be defined separately for lithiation and delithiation. In lithiation experiments with $J_0>0$, the local concentration generally increases with time, and the maximum concentration is typically attained near the end of the experiment and close to the electrochemically active surface. In delithiation experiments with $J_0<0$, the local concentration decreases with time, and the minimum concentration is typically attained near the end of the experiment and close to the electrochemically active surface.

5.2. Curvature and Constitutive Relation for Stress

Lithium insertion or extraction produces a local stress-free volumetric expansion or contraction in a freestanding electrode. This response is described by an isotropic chemical strain,

${\epsilon }_{ij}^{\mathrm{c}\mathrm{h}}=\beta \left(c−c_0\right){\delta }_{ij},$ (5-12)

where $\beta$ is the coefficient of compositional expansion, $c_0$ is the reference concentration for the considered experiment group, and ${\delta }_{ij}$ is the Kronecker delta. For lithiation, $c_0=c_0^{\mathrm{i}n}$; for delithiation, $c_0=c_0^{\mathrm{e}x}$. Physically, $\beta$ represents the linear strain per unit change in lithium concentration and can be estimated experimentally from strain measurements during low-rate lithiation/delithiation or after relaxation [30].

In the electrode/substrate bilayer structure, because the electrode film is perfectly bonded to the substrate, in-plane expansion or contraction is constrained. This constraint generates in-plane stress: generally compressive during lithiation and tensile during delithiation. If the bilayer is free to bend, this stress drives curvature evolution.

Within classical beam theory, the in-plane strain varies linearly across the thickness. The in-plane mechanical strains in the active electrode and current collector are written as [31]

${\epsilon }_{\parallel }^{m,e}\left(z,t\right)={\epsilon }_0\left(t\right)+\kappa \left(t\right)z-\beta \left[c\left(z,t\right)-c_0\right],0\le z\le L,$ (5-13a)

${\epsilon }_{\parallel }^{cu}\left(z,t\right)={\epsilon }_0\left(t\right)+\kappa \left(t\right)z,-h\le z\le 0,$ (5-13b)

where ${\epsilon }_0\left(t\right)$ is the in-plane strain at the electrode/current collector interface and $\kappa \left(t\right)$ is the curvature [32]. With the convention $\epsilon ={\epsilon }_0+\kappa z$, positive $\kappa$ means that the in-plane strain increases from the current collector side to the electrode surface side.

The biaxial modulus is defined as

$Y={\displaystyle \frac{E}{1-\nu }},$ (5-14)

where $E$ is Young’s modulus and $\nu$ is Poisson’s ratio.

The current collector is assumed to have a constant biaxial modulus $Y_{cu}$. In contrast, the biaxial modulus of the active electrode, which is factually the apparent biaxial modulus of the active electrode, generally depends on lithium concentration. Thus,

$Y_e=Y_e\left(c\right(z,t\left)\right)\equiv Y_e\left(c\right)$. (5-15)

The equi-biaxial stresses in the active electrode and current collector are therefore [33]

${\sigma }_{\parallel }^e\left(z,t\right)=Y_e\left(c\left(z,t\right)\right)\left[{\epsilon }_0\left(t\right)+\kappa \left(t\right)z-\beta \left(c\left(z,t\right)-c_0\right)\right],0\le z\le L,$ (5-16a)

${\sigma }_{\parallel }^{cu}(z,t)=Y_{cu}\left[{\epsilon }_0\right(t)+\kappa (t\left)z\right],-h\le z\le 0.$ (5-16b)

In the absence of external mechanical loading, force and moment balance require [28]

${\int }_{-h}^0{\sigma }_{\parallel }^{cu}(z,t)dz+{\int }_0^L{\sigma }_{\parallel }^e(z,t)dz=0,$ (5-17a)

${\int }_{-h}^0{\sigma }_{\parallel }^{cu}(z,t)zdz+{\int }_0^L{\sigma }_{\parallel }^e(z,t)zdz=0.$ (5-17b)

To express the governing equations compactly, the stiffness moments of the active layer are defined as

$A_m\left(t\right)={\int }_0^L{Y}_e\left(c\right(z,t\left)\right)z^{m}dz,m=\mathrm{0,1},2,$ (5-18)

and the chemical loading moments as

$B_m\left(t\right)={\int }_0^L{Y}_e\left(c\right(z,t\left)\right)\beta \left[c\right(z,t)-c_0]z^{m}dz,m=\mathrm{0,1}.$ (5-19)

Substitution into the force and moment balance equations gives

$\left[Y_{cu}h+A_0\left(t\right)\right]{\epsilon }_0\left(t\right)+\left[-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+A_1\left(t\right)\right]\kappa \left(t\right)=B_0\left(t\right),$ (5-20a)

$\left[-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+A_1\left(t\right)\right]{\epsilon }_0\left(t\right)+\left[{\displaystyle \frac{Y_{cu}{h}^3}{3}}+A_2\left(t\right)\right]\kappa \left(t\right)=B_1\left(t\right).$ (5-20b)

Equivalently,

$\left[\begin{array}{cc}{K}_{00}\left(t\right)& K_{01}\left(t\right)\\ K_{01}\left(t\right)& K_{11}\left(t\right)\end{array}\right]\left[\begin{array}{c}{\epsilon }_0\left(t\right)\\ \kappa \left(t\right)\end{array}\right]=\left[\begin{array}{c}{B}_0\left(t\right)\\ B_1\left(t\right)\end{array}\right],$ (5-21)

where

$K_{00}\left(t\right)=Y_{cu}h+A_0\left(t\right),$ (5-22a)

$K_{01}\left(t\right)=-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+A_1\left(t\right),$ (5-22b)

$K_{11}\left(t\right)={\displaystyle \frac{Y_{cu}{h}^3}{3}}+A_2\left(t\right).$ (5-22c)

Thus, for a prescribed constitutive relation $Y_e\left(c\right)$, the curvature predicted by the variable-modulus model is

${\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\left(t;Y_e\right)={\displaystyle \frac{K_{00}\left(t\right)B_1\left(t\right)-K_{01}\left(t\right)B_0\left(t\right)}{K_{00}\left(t\right)K_{11}\left(t\right)-K_{01}^2\left(t\right)}}.$ (5-23)

This equation defines the forward operator that maps the modulus-concentration relation to the curvature response,

$Y_e\left(c\right)\mapsto {\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\left(t;Y_e\right).$ (5-24)

Once ${\epsilon }_0\left(t\right)$ and $\kappa \left(t\right)$ are obtained, the local stress distribution follows directly from Eqs. (5-16a) and (5-16b). The average stress in the active electrode is

${\overline{\sigma }}_{\parallel }^e\left(t\right)={\displaystyle \frac{1}{L}}{\int }_0^L{Y}_e\left(c\left(z,t\right)\right)\left[{\epsilon }_0\left(t\right)+\kappa \left(t\right)z-\beta \left(c\left(z,t\right)-c_0\right)\right]dz,$ (5-25)

or equivalently,

${\overline{\sigma }}_{\parallel }^e\left(t\right)={\displaystyle \frac{1}{L}}\left[A_0\left(t\right){\epsilon }_0\left(t\right)+A_1\left(t\right)\kappa \left(t\right)-B_0\left(t\right)\right].$ (5-26)

The maximum absolute stress in the active electrode is evaluated as

${\sigma }_{max}^e\left(t\right)=\underset{0\le z\le L}{\mathrm{m}\mathrm{a}\mathrm{x}}\mid Y_e\left(c\right(z,t\left)\right)\left[{\epsilon }_0\left(t\right)+\kappa \left(t\right)z-\beta \left(c\right(z,t)-c_0)\right]\mid .$ (5-27)

Because $Y_e\left(c\right(z,t\left)\right)$ varies through the thickness, the maximum stress does not necessarily occur at $z=0$ or $z=L$.

5.3. Independent Inverse Identification of Modulus from Measured Curvature

In principle, if the curvature history ${\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(t\right)$ is measured in situ, Eq. (5-23) can be used to identify the unknown modulus-concentration relation $Y_e\left(c\right)$. However, direct pointwise recovery of $Y_e\left(c\right)$ from curvature data is generally ill-posed. The curvature at each time is a scalar quantity, whereas $Y_e\left(c\right)$ is an unknown continuous function. Moreover, the curvature response depends on $Y_e\left(c\right)$ through the thickness-weighted integrals $A_0\left(t\right)$, $A_1\left(t\right)$, $A_2\left(t\right)$, $B_0\left(t\right)$, and $B_1\left(t\right)$. Consequently, different functions $Y_e\left(c\right)$ may produce nearly indistinguishable curvature histories.

To address this non-uniqueness, a two-step inverse identification strategy is adopted. First, an apparent homogeneous modulus is extracted from the experimental curvature data using an effective-constant-modulus approximation. This reduced-order identification provides the magnitude and trend of the modulus evolution. Second, the inferred apparent modulus-concentration trend is used to guide the choice of candidate functional forms for $Y_e\left(c\right)$. Finally, the full variable-modulus model is used to identify a more accurate local modulus-concentration relation through parameterized and, when necessary, regularized inverse analysis.

In the present work, lithiation and delithiation experiments are analyzed independently. For lithiation, all experiments are assumed to start from the same initial concentration $c_0^{\mathrm{i}n}$, and the applied flux satisfies $J_0>0$. Multiple lithiation experiments at different current densities are jointly fitted using a common modulus function

$Y_e^{\mathrm{i}n}(c;{\theta }_{\mathrm{i}n}).$ (5-28a)

For delithiation, all experiments are assumed to start from the same initial concentration $c_0^{\mathrm{e}x}$, and the applied flux satisfies $J_0<0$. Multiple delithiation experiments are jointly fitted using a separate modulus function

$Y_e^{\mathrm{e}x}(c;{\theta }_{\mathrm{e}x}).$ (5-28b)

No constraint is imposed that $Y_e^{\mathrm{i}n}\left(c\right)$ and $Y_e^{\mathrm{e}x}\left(c\right)$ must be identical, unless fully reversible elastic behavior without path dependence, hysteresis, damage, or irreversible structural change is assumed.

Step I: Apparent modulus identification using an effective constant modulus

Strictly, because the lithium concentration varies through the electrode thickness, the active-layer modulus should be written as $Y_e\left(c\right(z,t\left)\right)$. A fully rigorous treatment therefore requires a specified constitutive relation $Y_e\left(c\right)$ and solution of the equilibrium equations with variable coefficients. To retain analytical tractability in the first stage, the spatially varying modulus field is approximated by an apparent homogeneous modulus at each time,

$Y_e\left(c\right(z,t\left)\right)\approx {\bar{Y}}_e\left(t\right).$ (5-29)

The quantity ${\bar{Y}}_e\left(t\right)$ should be interpreted as an apparent or effective biaxial modulus that reproduces the measured curvature within the homogeneous-modulus approximation. It is not, in general, identical to the local material function $Y_e\left(c\right)$.

With Eq. (5-29), the force and moment balance equations reduce to

$\left[Y_{cu}h+{\bar{Y}}_e\left(t\right)L\right]{\epsilon }_0\left(t\right)+\left[-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+{\displaystyle \frac{{\bar{Y}}_e\left(t\right)L^2}{2}}\right]\kappa \left(t\right)={\bar{Y}}_e\left(t\right)\beta J_{0}t,$ (5-30a)

$\left[-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+{\displaystyle \frac{{\bar{Y}}_e\left(t\right)L^2}{2}}\right]{\epsilon }_0\left(t\right)+\left[{\displaystyle \frac{Y_{cu}{h}^3}{3}}+{\displaystyle \frac{{\bar{Y}}_e\left(t\right)L^3}{3}}\right]\kappa \left(t\right)=\beta {\bar{Y}}_e\left(t\right)Q\left(t\right),$ (5-30b)

where

$Q\left(t\right)={\displaystyle \frac{J_{0}Lt}{2}}+{\displaystyle \frac{J_0{L}^3}{24D}}-{\displaystyle \frac{2J_0{L}^3}{D{\pi }^4}}{\displaystyle \sum _{n=1}^{\mathrm{\infty }}}{\displaystyle \frac{1-(-1{)}^n}{n^4}}\mathrm{e}\mathrm{x}\mathrm{p}(-{\displaystyle \frac{n^2{\pi }^{2}Dt}{L^2}})$. (5-31)

The sign of $J_0$ is retained in Eq. (5-31). Thus, the same expression applies to both lithiation and delithiation, provided that $J_0>0$ is used for lithiation and $J_0<0$ is used for delithiation.

Given the measured curvature ${\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(t_k\right)$, Eqs. (5-30a) and (5-30b) can be solved at each time $t_k$ for the two unknowns ${\epsilon }_0\left(t_k\right)$ and ${\bar{Y}}_e\left(t_k\right)$, under the physical constraint

${\bar{Y}}_e\left(t_k\right)>0.$ (5-32)

Data points at very early times may be excluded from Step I because the curvature signal is too small to identify the apparent modulus robustly.

For lithiation datasets, this procedure gives

$\left\{t_k,{\bar{Y}}_e^{\mathrm{i}n}\left(t_k\right)\right\},$ (5-33a)

and for delithiation datasets, it gives

$\left\{t_k,{\bar{Y}}_e^{\mathrm{e}x}\left(t_k\right)\right\}.$ (5-33b)

Using the thickness-averaged concentration,

${\bar{c}}_k=c_0+{\displaystyle \frac{J_0{t}_k}{L}},$ (5-34)

the apparent modulus history can be converted into an apparent modulus-concentration relation. For lithiation,

$\left\{{\bar{c}}_k^{\mathrm{i}n},{\bar{Y}}_e^{\mathrm{i}n}\left({\bar{c}}_k\right)\right\},{\bar{c}}_k^{\mathrm{i}n}=c_0^{\mathrm{i}n}+{\displaystyle \frac{J_0^{\mathrm{i}n}{t}_k}{L}},$ (5-35a)

where $J_0^{\mathrm{i}n}>0$.

For delithiation,

$\left\{{\bar{c}}_k^{\mathrm{e}x},{\bar{Y}}_e^{\mathrm{e}x}\left({\bar{c}}_k\right)\right\},{\bar{c}}_k^{\mathrm{e}x}=c_0^{\mathrm{e}x}+{\displaystyle \frac{J_0^{\mathrm{e}x}{t}_k}{L}},$ (5-35b)

where $J_0^{\mathrm{e}x}<0$.

These apparent relations provide useful reduced-order descriptors of the macroscopic modulus evolution. If the apparent modulus-concentration relation is approximately linear, a linear or weakly nonlinear model may be used as the initial candidate. If it decreases monotonically and tends to saturate, a saturating exponential model may be particularly appropriate. Logistic, Hill-type, or monotonic spline representations may also be considered when more complex transitions are observed and sufficient data are available.

The apparent relation ${\bar{Y}}_e\left(\bar{c}\right)$ is used only to guide the subsequent full inversion. It should not be interpreted as the exact local constitutive relation $Y_e\left(c\right)$, because the measured curvature is sensitive to thickness-weighted integrals of $Y_e\left(c\right(z,t\left)\right)$, rather than to its pointwise value.

Step II: Parameterized full inversion of${\mathit{Y}}_{\mathit{e}}\left(\mathit{c}\right)$

In the second step, the local concentration-dependent modulus is represented by a finite-dimensional parameterization. Since lithiation and delithiation are treated as independent inverse problems, two separate functions are introduced:

$Y_e^{\mathrm{i}n}(c;{\theta }_{\mathrm{i}n})$ (5-36a)

for lithiation, and

$Y_e^{\mathrm{e}x}(c;{\theta }_{\mathrm{e}x})$ (5-36b)

for delithiation.

The choice of parameterization is guided by the apparent trends obtained in Step I.

For example, if the apparent modulus is nearly linear, one may use

$Y_e^p(c;{\theta }_p)={\theta }_0^p+{\theta }_1^{p}c,{\theta }_p=\left[\begin{array}{c}{\theta }_0^p\\ {\theta }_1^p\end{array}\right],$ (5-37)

where $p=\mathrm{i}n$ or $p=\mathrm{e}x$.

If the modulus evolves in an exponential form, a quadratic exponential model may be adopted,

$Y_e^p(c;{\theta }_p)=\mathrm{e}\mathrm{x}\mathrm{p}\left({\eta }_0^p+{\theta }_1^p{s}_p+{\theta }_2^p{s}_p^2\right),$ (5-38)

where

${\theta }_p=\left[\begin{array}{c}{\eta }_0^p\\ {\theta }_1^p\\ {\theta }_2^p\end{array}\right].$ (5-39)

Here $s_p$ is the normalized concentration for the independently analyzed process $p$. The normalization range must be determined separately for lithiation and delithiation from the local concentration fields.

For lithiation, $p=\mathrm{i}n$, the normalized concentration is defined as

$s_{\mathrm{i}n}={\displaystyle \frac{c-c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{i}n}}{c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{i}n}-c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{i}n}}}.$ (5-40)

Because all lithiation experiments start from $c_0^{\mathrm{i}n}$, the lower bound is taken as

$c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{i}n}=c_0^{\mathrm{i}n}.$ (5-41)

For multiple lithiation experiments, $q=1,\dots ,Q_{\mathrm{i}n}$, the local concentration field of the $q$-th experiment is

$c_{\mathrm{i}n}^{\left(q\right)}(z,t)=c_0^{\mathrm{i}n}+{\displaystyle \frac{J_{0,\mathrm{i}n}^{\left(q\right)}t}{L}}+{\displaystyle \frac{J_{0,\mathrm{i}n}^{\left(q\right)}{z}^2}{2DL}}-{\displaystyle \frac{J_{0,\mathrm{i}n}^{\left(q\right)}L}{6D}}-{\displaystyle \frac{2J_{0,\mathrm{i}n}^{\left(q\right)}L}{D{\pi }^2}}{\displaystyle \sum _{n=1}^{\mathrm{\infty }}}{\displaystyle \frac{{\left(-1\right)}^n}{n^2}}\mathrm{e}\mathrm{x}\mathrm{p}(-{\displaystyle \frac{n^2{\pi }^{2}Dt}{L^2}})\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{n\pi z}{L}}\right),$ (5-42)

with

$J_{0,\mathrm{i}n}^{\left(q\right)}>0.$ (5-43)

The maximum local concentration sampled by all lithiation experiments is

$c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}^{\mathrm{i}n}=\underset{q=1,\dots ,Q_{\mathrm{i}n}}{\mathrm{m}\mathrm{a}\mathrm{x}}\left[\underset{0\le z\le L}{\mathrm{m}\mathrm{a}\mathrm{x}}{c}_{\mathrm{i}n}^{\left(q\right)}\left(z,t_{\mathrm{e}\mathrm{n}\mathrm{d}}^{\left(q\right)}\right)\right].$ (5-44)

To avoid numerical extrapolation caused by discretization or series truncation errors, a small safety factor $\alpha \ge 1$ is introduced:

$c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{i}n}=c_0^{\mathrm{i}n}+\alpha \left(c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}^{\mathrm{i}n}−c_0^{\mathrm{i}n}\right).$ (5-45)

Typically, $\alpha$ is chosen between 1.01 and 1.05.

For delithiation, $p=\mathrm{e}x$, the normalized concentration is defined as

$s_{\mathrm{e}x}={\displaystyle \frac{c-c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{e}x}}{c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{e}x}-c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{e}x}}}.$ (5-46)

If all delithiation experiments start from $c_0^{\mathrm{e}x}$, the upper bound can be taken as

$c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{e}x}=c_0^{\mathrm{e}x}.$ (5-47)

For multiple delithiation experiments, $q=1,\dots ,Q_{\mathrm{e}x}$, the local concentration field of the $q$-th experiment is

$c_{\mathrm{e}x}^{\left(q\right)}(z,t)=c_0^{\mathrm{e}x}+{\displaystyle \frac{J_{0,\mathrm{e}x}^{\left(q\right)}t}{L}}+{\displaystyle \frac{J_{0,\mathrm{e}x}^{\left(q\right)}{z}^2}{2DL}}-{\displaystyle \frac{J_{0,\mathrm{e}x}^{\left(q\right)}L}{6D}}-{\displaystyle \frac{2J_{0,\mathrm{e}x}^{\left(q\right)}L}{D{\pi }^2}}{\displaystyle \sum _{n=1}^{\mathrm{\infty }}}{\displaystyle \frac{{\left(-1\right)}^n}{n^2}}\mathrm{e}\mathrm{x}\mathrm{p}(-{\displaystyle \frac{n^2{\pi }^{2}Dt}{L^2}})\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{n\pi z}{L}}\right),$ (5-48)

with

$J_{0,\mathrm{e}x}^{\left(q\right)}<0.$ (5-49)

The minimum local concentration sampled by all delithiation experiments is

$c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}^{\mathrm{e}x}=\underset{q=1,\dots ,Q_{\mathrm{e}x}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left[\underset{0\le z\le L}{\mathrm{m}\mathrm{i}\mathrm{n}}{c}_{\mathrm{e}x}^{\left(q\right)}\left(z,t_{\mathrm{e}\mathrm{n}\mathrm{d}}^{\left(q\right)}\right)\right].$ (5-50)

A safety factor $\alpha \ge 1$ is introduced as

$c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^{\mathrm{e}x}=c_0^{\mathrm{e}x}-\alpha \left(c_0^{\mathrm{e}x}−c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}^{\mathrm{e}x}\right).$ (5-51)

Again, $\alpha$ is typically chosen between 1.01 and 1.05.

In addition to the quadratic exponential model, a saturating exponential model may be used when the apparent modulus changes rapidly at low lithium concentration and then approaches a plateau. For either lithiation or delithiation, this model is written as

$Y_e^p(c;{\theta }_p)=Y_{\mathrm{\infty }}^p+\left(Y_0^p−Y_{\mathrm{\infty }}^p\right)\mathrm{e}\mathrm{x}\mathrm{p}(-k^p{s}_p),$ (5-52)

where $Y_0^p$ is the modulus at $s_p=0$, $Y_{\mathrm{\infty }}^p$ is the saturated modulus at high $s_p$, and $k^p$ is the transition rate. The superscript $p$ denotes either lithiation or delithiation.

To guarantee positivity of the parameters, a logarithmic parameterization can be introduced:

$Y_0^p=\mathrm{e}\mathrm{x}\mathrm{p}\left(a^p\right),Y_{\mathrm{\infty }}^p=\mathrm{e}\mathrm{x}\mathrm{p}\left(b^p\right),k^p=\mathrm{e}\mathrm{x}\mathrm{p}\left(g^p\right).$ (5-53)

The saturating exponential model then becomes

$Y_e^p(c;{\theta }_p)=\mathrm{e}\mathrm{x}\mathrm{p}\left(b^p\right)+\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(a^p\right)-\mathrm{e}\mathrm{x}\mathrm{p}\left(b^p\right)\right]\mathrm{e}\mathrm{x}\mathrm{p}\left[-\mathrm{e}\mathrm{x}\mathrm{p}\left(g^p\right)s_p\right],$ (5-54)

with

${\theta }_p=\left[\begin{array}{c}{a}^p\\ b^p\\ g^p\end{array}\right].$ (5-55)

More flexible choices, such as B-splines or monotonic splines, may also be used when sufficient experimental data are available.

For a given parameter vector ${\theta }_p$, where $p=\mathrm{i}n$ or $p=\mathrm{e}x$, the local concentration field $c_p^{\left(q\right)}(z,t)$ is first calculated from the analytical diffusion solution. The stiffness and chemical loading moments are then evaluated as

$A_m^{(q,p)}(t;{\theta }_p)={\int }_0^L{Y}_e^p\left(c_p^{\left(q\right)}(z,t);{\theta }_p\right)z^{m}dz,m=\mathrm{0,1},2,$ (5-56)

$B_m^{(q,p)}(t;{\theta }_p)={\int }_0^L{Y}_e^p\left(c_p^{\left(q\right)}(z,t);{\theta }_p\right)\beta \left[c_p^{\left(q\right)}(z,t)-c_0^p\right]z^{m}dz,m=\mathrm{0,1}.$ (5-57)

Here

$c_0^p=\left\{\begin{array}{cc}{c}_0^{\mathrm{i}n},& p=\mathrm{i}n,\\ c_0^{\mathrm{e}x},& p=\mathrm{e}x.\end{array}\right.$ (5-58)

The corresponding stiffness matrix terms are

$K_{00}^{(q,p)}(t;{\theta }_p)=Y_{cu}h+A_0^{(q,p)}(t;{\theta }_p),$ (5-59a)

$K_{01}^{(q,p)}(t;{\theta }_p)=-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+A_1^{(q,p)}(t;{\theta }_p),$ (5-59b)

$K_{11}^{(q,p)}(t;{\theta }_p)={\displaystyle \frac{Y_{cu}{h}^3}{3}}+A_2^{(q,p)}(t;{\theta }_p).$ (5-59c)

The predicted curvature is then computed as

${\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(q,p)}(t;{\theta }_p)={\displaystyle \frac{K_{00}^{(q,p)}(t;{\theta }_p)B_1^{(q,p)}(t;{\theta }_p)-K_{01}^{(q,p)}(t;{\theta }_p)B_0^{(q,p)}(t;{\theta }_p)}{K_{00}^{(q,p)}(t;{\theta }_p)K_{11}^{(q,p)}(t;{\theta }_p)-{\left[K_{01}^{(q,p)}(t;{\theta }_p)\right]}^2}}.$ (5-60)

To reduce the influence of initial curvature, residual stress, or mounting errors, curvature increments are used:

$\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,p)}\left(t_k\right)={\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,p)}\left(t_k\right)-{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,p)}\left(0\right),$ (5-61)

$\mathsf{\Delta }{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(q,p)}(t_k;{\theta }_p)={\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(q,p)}(t_k;{\theta }_p)-{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(q,p)}(0;{\theta }_p).$ (5-62)

For lithiation, the objective function of the inverse problem is formulated as

$\underset{{\theta }_{\mathrm{i}n}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\mathsf{\Phi }}_{\mathrm{i}n}\left({\theta }_{\mathrm{i}n}\right)={\displaystyle \sum _{q=1}^{Q_{\mathrm{i}n}}}{\displaystyle \sum _{k=1}^{N_q}}{w}_{qk}^{\mathrm{i}n}{\left[\mathsf{\Delta }{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(\mathrm{q},\mathrm{i}\mathrm{n})}(t_k;{\theta }_{\mathrm{i}n})-\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,in)}\left(t_k\right)\right]}^2+{\mathsf{\Phi }}_{\mathrm{r}\mathrm{e}\mathrm{g}}^{\mathrm{i}n}\left({\theta }_{\mathrm{i}n}\right).$ (5-63)

For delithiation, the objective function of the inverse problem is formulated independently as

$\underset{{\theta }_{\mathrm{e}x}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\mathsf{\Phi }}_{\mathrm{e}x}\left({\theta }_{\mathrm{e}x}\right)={\displaystyle \sum _{q=1}^{Q_{\mathrm{e}x}}}{\displaystyle \sum _{k=1}^{N_q}}{w}_{qk}^{\mathrm{e}x}{\left[\mathsf{\Delta }{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(q,ex)}(t_k;{\theta }_{\mathrm{e}x})-\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,ex)}\left(t_k\right)\right]}^2+{\mathsf{\Phi }}_{\mathrm{r}\mathrm{e}\mathrm{g}}^{\mathrm{e}x}\left({\theta }_{\mathrm{e}x}\right).$ (5-64)

$w_{qk}^p$ is a weighting factor, and a useful weighting choice is

$w_{qk}^p={\displaystyle \frac{1}{N_q{\left[{\kappa }_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}}^{(q,p)}\right]}^2}},p=\mathrm{i}n,\mathrm{e}x,$ (5-65)

where

${\kappa }_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}}^{(q,p)}=\underset{k}{\mathrm{m}\mathrm{a}\mathrm{x}}\mid \mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,p)}\left(t_k\right)\mid .$ (5-66)

This weighting balances the contributions from experiments with different numbers of data points and different curvature magnitudes.

For low-dimensional models, such as the linear model, quadratic exponential model, and saturating exponential model, ${\mathsf{\Phi }}_{\mathrm{r}\mathrm{e}\mathrm{g}}^p$ may be omitted. For spline-based representations, regularization is required to suppress nonphysical oscillations. A typical smoothness penalty is

${\mathsf{\Phi }}_{\mathrm{r}\mathrm{e}\mathrm{g}}^p\left({\theta }_p\right)={\lambda }_p{\int }_{c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^p}^{c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}^p}{\left[{\displaystyle \frac{d^2{Y}_e^p(c;{\theta }_p)}{d{c}^2}}\right]}^{2}dc,p=\mathrm{i}n,\mathrm{e}x,$ (5-67)

where ${\lambda }_p$ controls the strength of regularization.

The inverse problem becomes the optimization of the objective functions. Physical constraints should be imposed during optimization. First, the modulus must remain positive over the fitting interval:

$Y_e^p(c;{\theta }_p)\ge Y_{\mathrm{m}\mathrm{i}\mathrm{n}}>0,c\in \left[c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^p,c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}^p\right],p=\mathrm{i}n,\mathrm{e}x.$ (5-68)

For the linear model,

$Y_e^p(c;{\theta }_p)={\theta }_0^p+{\theta }_1^{p}c,$ (5-69)

this positivity condition reduces to the two endpoint constraints:

${\theta }_0^p+{\theta }_1^p{c}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^p\ge Y_{\mathrm{m}\mathrm{i}\mathrm{n}},$ (5-70a)

${\theta }_0^p+{\theta }_1^p{c}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}^p\ge Y_{\mathrm{m}\mathrm{i}\mathrm{n}}.$ (5-70b)

If the electrode is known to soften with increasing lithium concentration, the monotonicity condition is

${\displaystyle \frac{d{Y}_e^p}{dc}}\le 0.$ (5-71)

For the linear model, this becomes

${\theta }_1^p\le 0.$ (5-72)

If the electrode is known to harden with increasing lithium concentration, the monotonicity condition is

${\displaystyle \frac{d{Y}_e^p}{dc}}\ge 0,$ (5-73)

which gives

${\theta }_1^p\ge 0.$ (5-74)

For the quadratic exponential model,

$Y_e^p(c;{\theta }_p)=\mathrm{e}\mathrm{x}\mathrm{p}\left({\eta }_0^p+{\theta }_1^p{s}_p+{\theta }_2^p{s}_p^2\right),$ (5-75)

positivity is automatically satisfied. Since

${\displaystyle \frac{d{Y}_e^p}{d{s}_p}}=Y_e^p\left(s_p\right)\left({\theta }_1^p+2{\theta }_2^p{s}_p\right),$ (5-76)

monotonic softening with increasing lithium concentration requires

${\theta }_1^p\le 0,{\theta }_1^p+2{\theta }_2^p\le 0.$ (5-77)

Monotonic hardening with increasing lithium concentration requires

${\theta }_1^p\ge 0,{\theta }_1^p+2{\theta }_2^p\ge 0.$ (5-78)

For the saturating exponential model,

$Y_e^p(c;{\theta }_p)=Y_{\mathrm{\infty }}^p+\left(Y_0^p−Y_{\mathrm{\infty }}^p\right)\mathrm{e}\mathrm{x}\mathrm{p}(-k^p{s}_p),$ (5-79)

the derivative with respect to $s_p$ is

${\displaystyle \frac{d{Y}_e^p}{d{s}_p}}=-k^p\left(Y_0^p−Y_{\mathrm{\infty }}^p\right)\mathrm{e}\mathrm{x}\mathrm{p}(-k^p{s}_p).$ (5-80)

Therefore, if

$Y_0^p>Y_{\mathrm{\infty }}^p>0,k^p>0,$ (5-81)

then

${\displaystyle \frac{d{Y}_e^p}{dc}}\le 0,$ (5-82)

and the modulus decreases monotonically with increasing lithium concentration. Under the logarithmic parameterization,

$Y_0^p=\mathrm{e}\mathrm{x}\mathrm{p}\left(a^p\right),Y_{\mathrm{\infty }}^p=\mathrm{e}\mathrm{x}\mathrm{p}\left(b^p\right),k^p=\mathrm{e}\mathrm{x}\mathrm{p}\left(g^p\right),$ (5-83)

this softening condition is imposed as

$a^p>b^p.$ (5-84)

For lithiation-induced hardening,

$Y_{\mathrm{\infty }}^p>Y_0^p>0,k^p>0,$ (5-85)

which corresponds to

$b^p>a^p.$ (5-86)

Additional bounds may also be imposed:

$Y_{\mathrm{m}\mathrm{i}\mathrm{n}}\le Y_e^p(c;{\theta }_p)\le Y_{\mathrm{m}\mathrm{a}\mathrm{x}},c\in \left[c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}^p,c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}^p\right].$ (5-87)

It should be noted that the monotonicity condition is imposed with respect to lithium concentration $c$, not with respect to time. Therefore, if the electrode softens with increasing lithium concentration, the condition $d{Y}_e/dc\le 0$ applies to both lithiation and delithiation inverse analyses, even though the concentration decreases with time during delithiation.

After optimization, the quality of the identified lithiation model should be evaluated using the residuals

$r_{qk}^{\mathrm{i}n}=\mathsf{\Delta }{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(q,in)}(t_k;{\widehat{\theta }}_{\mathrm{i}n})-\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,in)}\left(t_k\right),$ (5-88a)

where ${\widehat{\theta }}_{\mathrm{i}n}$ is the optimal lithiation parameter vector. Similarly, for delithiation,

$r_{qk}^{\mathrm{e}x}=\mathsf{\Delta }{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(q,ex)}(t_k;{\widehat{\theta }}_{\mathrm{e}x})-\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,ex)}\left(t_k\right),$ (5-88b)

where ${\widehat{\theta }}_{\mathrm{e}x}$ is the optimal delithiation parameter vector.

For a given process $p=\mathrm{i}n$ or $p=\mathrm{e}x$, the root-mean-square error for the $q$-th experiment is

${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}^{(q,p)}=\sqrt{{\displaystyle \frac{1}{N_q}}{\displaystyle \sum _{k=1}^{N_q}}{\left[r_{qk}^p\right]}^2}.$ (5-89)

The relative root-mean-square error is

${\mathrm{R}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}^{(q,p)}=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{k=1}^{N_q}}{\left[r_{qk}^p\right]}^2}{{\displaystyle \sum _{k=1}^{N_q}}{\left[\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,p)}\left(t_k\right)\right]}^2}}}.$ (5-90)

A suitable modulus model should reproduce the curvature histories for all experiments within the same process group with small and randomly distributed residuals. Systematic residual trends may indicate that the assumed functional form of $Y_e^p\left(c\right)$ is inadequate, or that other parameters such as $D$, $\beta$, or the mechanical assumptions of the bilayer model require refinement.

After the local constitutive relations $Y_e^{\mathrm{i}n}\left(c\right)$ and/or $Y_e^{\mathrm{e}x}\left(c\right)$ have been identified, the stress field can be calculated using Eqs. (5-16a) and (5-16b). The average active-layer stress and maximum absolute active-layer stress are then calculated using Eqs. (5-26) and (5-27), respectively. If only curvature increments are fitted, the resulting stress should be interpreted relative to the chosen initial state of each independently analyzed experiment group. If absolute stress is required, the chemical strain should be referred to a true stress-free reference concentration rather than only to the experimental initial concentration.

5.4. A Further Clarification

With the above formulation, the lithiation and delithiation inverse problems are independent:

${\theta }_{\mathrm{i}n}\ne {\theta }_{\mathrm{e}x}$ $\mathrm{in}\mathrm{general},$ (5-91)

and

$Y_e^{\mathrm{i}n}\left(c\right)\ne Y_e^{\mathrm{e}x}\left(c\right)$ $\mathrm{in}\mathrm{general}.$ (5-92)

The two functions may become identical only if one explicitly assumes fully reversible mechanical behavior without path dependence, hysteresis, damage, or irreversible microstructural evolution. Otherwise, the present framework allows lithiation and delithiation to produce distinct effective local modulus-concentration relations while still enabling multi-rate joint fitting within each individual process.

6.1. Physical System

The physical system is the same as Section V. We still use lithium as an example. The problem is still treated as one-dimensional in the thickness direction. The chemical diffusion coefficient $D$ and the thickness of the electrode $L$ are still taken to be constant and independent of concentration. Under the same assumptions as Section V, the active electrode still develops an equi-biaxial in-plane stress state, Equation (5-1), and the one-dimensional Fick’s diffusion equation keeps unchanged as

${\displaystyle \frac{\partial c}{\partial t}}=D{\displaystyle \frac{{\partial }^{2}c}{\partial z^2}},0<z<L,t>0.$ (6-1)

The initial concentration field at the beginning of the cycling experiment is prescribed as

$c\left(z,0\right)=c_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}\left(z\right),0\le z\le L.$ (6-2)

In many experiments, the electrode is initially equilibrated, in which case

$c\left(z,0\right)=c_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}},$ (6-3)

where $c_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$ is a spatially uniform constant. The framework, however, also permits a nonuniform initial concentration field if the cycling experiment starts from a partially relaxed or previously cycled state.

The boundary condition at $z=0$ remains to be

${\displaystyle \frac{\partial c}{\partial z}}(0,t)=0.$ (6-4)

The surface flux is prescribed by the applied galvanostatic current density. Since the protocol may contain several lithiation-delithiation cycles, we let $i\left(t\right)$ denote the signed current density, $\mathbb{n}$ the number of electrons transferred per lithium ion, and $F$ Faraday's constant. The molar lithium flux is then given by

$J\left(t\right)={\displaystyle \frac{i\left(t\right)}{\mathbb{n}F}}.$ (6-5)

The sign convention is chosen such that

$J\left(t\right)>0,$ for inward lithium flux into the electrode during lithiation, whereas

$J\left(t\right)<0$, for outward lithium flux during delithiation.

The electrochemical boundary condition at $z=L$ thus is

$D{\displaystyle \frac{\partial c}{\partial z}}\left(L,t\right)=J\left(t\right).$ (6-6)

This formulation applies to both lithiation-delithiation cycling and delithiation-lithiation cycling. For a lithiation-delithiation cycle, the flux sequence begins with a positive-flux segment followed by a negative-flux segment. For a delithiation-lithiation cycle, the sequence begins with a negative-flux segment followed by a positive-flux segment. The number of cycles may be one or greater than one.

The applied current is assumed to be fully converted into lithium insertion or extraction at the active surface. The thickness-averaged concentration is defined as

$\bar{c}\left(t\right)={\displaystyle \frac{1}{L}}{\int }_0^{L}c(z,t)dz.$ (6-7)

Integrating Eq. (6-1) through the film thickness gives

${\displaystyle \frac{d\bar{c}}{dt}}={\displaystyle \frac{J\left(t\right)}{L}}.$ (6-8)

Hence, for an arbitrary galvanostatic cycling protocol,

$\bar{c}\left(t\right)=\bar{c}\left(0\right)+{\displaystyle \frac{1}{L}}{\int }_0^{t}J\left(\tau \right)d\tau .$ (6-9)

For a piecewise-constant cycling protocol, if

$J\left(t\right)=J_r,t_r\le t<t_{r+1},$ (6-10)

then within the $r$-th segment,

$\bar{c}\left(t\right)=\bar{c}\left(t_r\right)+{\displaystyle \frac{J_r}{L}}\left(t-t_r\right).$ (6-11)

6.2. Concentration Field Under Galvanostatic Cycling

For a single-step lithiation or delithiation process beginning from a uniform concentration, the diffusion equation admits a closed-form analytical solution Equation (5-11). In cycling, however, the initial concentration field of each half-cycle is generally the nonuniform terminal concentration field inherited from the preceding half-cycle. Therefore, the single-step solution cannot be restarted from a uniform initial condition unless sufficient relaxation is imposed between segments.

To deal with arbitrary one-cycle or multi-cycle galvanostatic protocols, the cycling history is divided into $\mathbb{R}$ constant-flux segments. The $r$-th segment occupies the time interval

$t_r\le t<t_{r+1},$ (6-12)

with local time

$\tau =t-t_r.$ (6-13)

The imposed flux in this segment is

$J\left(t\right)=J_r.$ (6-14)

Here $J_r>0$ corresponds to lithiation and $J_r<0$ corresponds to delithiation.

The local concentration field during segment $r$ is denoted by

$c_r(z,\tau )=c(z,t_r+\tau ).$ (6-15)

It satisfies

${\displaystyle \frac{\partial c_r}{\partial \tau }}=D{\displaystyle \frac{{\partial }^2{c}_r}{\partial z^2}},0<z<L,$ (6-16)

with boundary conditions becoming

${\displaystyle \frac{\partial c_r}{\partial z}}(0,\tau )=0,$ (6-17)

$D{\displaystyle \frac{\partial c_r}{\partial z}}(L,\tau )=J_r.$ (6-18)

The initial condition for segment $r$ now is inherited from the previous segment:

$c_r\left(z,0\right)=c\left(z,t_r\right).$ (6-19)

For $r=1$,

$c_1(z,0)=c_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}\left(z\right)$. (6-20)

For $r>1$,

$c_r(z,0)=c_{r-1}(z,\mathsf{\Delta }{t}_{r-1}),$ (6-21)

where

$\mathsf{\Delta }{t}_{r-1}=t_r-t_{r-1}.$ (6-22)

For a constant flux $J_r$, the concentration field during the $r$-th segment can be expressed as

$c_r(z,\tau )={\displaystyle \frac{J_r\tau }{L}}+{\displaystyle \frac{J_r}{D}}\left({\displaystyle \frac{z^2}{2L}}−{\displaystyle \frac{L}{6}}\right)+a_0^{\left(r\right)}+{\displaystyle \sum _{n=1}^{\mathrm{\infty }}}{a}_n^{\left(r\right)}\mathrm{e}\mathrm{x}\mathrm{p}(-{\displaystyle \frac{n^2{\pi }^{2}D\tau }{L^2}})\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{n\pi z}{L}}\right).$ (6-23)

The coefficients are determined from the initial concentration field of the segment. Here, we define

$p_r\left(z\right)={\displaystyle \frac{J_r}{D}}\left({\displaystyle \frac{z^2}{2L}}−{\displaystyle \frac{L}{6}}\right).$ (6-24)

Then, we have

$a_0^{\left(r\right)}={\displaystyle \frac{1}{L}}{\int }_0^L\left[c_r(z,0)-p_r\left(z\right)\right]dz,$ (6-25)

and, for $n\ge 1$,

$a_n^{\left(r\right)}={\displaystyle \frac{2}{L}}{\int }_0^L\left[c_r(z,0)-p_r\left(z\right)\right]\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{n\pi z}{L}}\right)dz.$ (6-26)

Equations (6-23)-(6-26) provide a recursive analytical representation of the concentration field during galvanostatic cycling. They can be applied to lithiation-delithiation cycling, delithiation-lithiation cycling, one-cycle experiments, and multi-cycle experiments.

If each half-cycle is followed by sufficiently long relaxation such that the concentration field becomes uniform before the next galvanostatic segment begins, then the initial condition for each segment may be approximated as

$c_r\left(z,0\right)=\bar{c}\left(t_r\right).$ (6-27)

Under this additional relaxation approximation, the single-step uniform-initial-condition solution may be used repeatedly. Without sufficient relaxation, however, the recursive concentration solution in Equations (6-23)-(6-26) should be used.

6.3. Curvature and Constitutive Relation for Stress

For the cycling problem, a single reference concentration is introduced and denoted by $c_{\mathrm{r}\mathrm{e}\mathrm{f}}.$ The chemical strain is assumed isotropic:

${\epsilon }_{ij}^{\mathrm{c}\mathrm{h}}=\beta \left[c\right(z,t)-c_{\mathrm{r}\mathrm{e}\mathrm{f}}]{\delta }_{ij},$ (6-28)

The reference concentration $c_{\mathrm{r}\mathrm{e}\mathrm{f}}$ should correspond to the stress-free chemical reference state whenever possible. In practice, it may be chosen as the initial equilibrated concentration, the fully delithiated concentration, or another experimentally justified reference state. If the true stress-free reference state is unknown, $c_{\mathrm{r}\mathrm{e}\mathrm{f}}$ may be treated as an additional parameter or absorbed into the interpretation of stress increments.

The in-plane mechanical strain in the active electrode layer is given by

${\epsilon }_{\parallel }^{m,e}\left(z,t\right)={\epsilon }_0\left(t\right)+\kappa \left(t\right)z-\beta \left[c\left(z,t\right)-c_{\mathrm{r}\mathrm{e}\mathrm{f}}\right],0\le z\le L,$ (6-29)

The in-plane strain in the current collector remains

${\epsilon }_{\parallel }^{cu}(z,t)={\epsilon }_0\left(t\right)+\kappa \left(t\right)z,-h\le z\le 0.$ (6-30)

In the present cycling framework, lithiation and delithiation are assumed to share the same concentration-dependent modulus relation, i.e.,

$Y_e^{in}\left(c\right)=Y_e^{ex}\left(c\right)=Y_e\left(c\right).$ (6-31)

Equivalently, in the inverse problem, we have

$Y_e\left(c\right)=Y_e\left(c;\theta \right),$ (6-32)

where $\theta$ is a single parameter vector used for all lithiation and delithiation segments.

The equi-biaxial stresses are therefore given by

${\sigma }_{\parallel }^e(z,t)=Y_e\left(c\right(z,t);\theta )\left[{\epsilon }_0\left(t\right)+\kappa \left(t\right)z-\beta \left(c\right(z,t)-c_{\mathrm{r}\mathrm{e}\mathrm{f}})\right],0\le z\le L,$ (6-33)

and

${\sigma }_{\parallel }^{cu}(z,t)=Y_{cu}\left[{\epsilon }_0\right(t)+\kappa (t\left)z\right],-h\le z\le 0.$ (6-34)

In the absence of external mechanical loading, force and moment balance require

${\int }_{-h}^0{\sigma }_{\parallel }^{cu}(z,t)dz+{\int }_0^L{\sigma }_{\parallel }^e(z,t)dz=0,$ (6-35)

${\int }_{-h}^0{\sigma }_{\parallel }^{cu}(z,t)zdz+{\int }_0^L{\sigma }_{\parallel }^e(z,t)zdz=0.$ (6-36)

The active electrode layer stiffness moments are defined as

$A_m(t;\theta )={\int }_0^L{Y}_e\left(c\right(z,t);\theta )z^{m}dz,m=\mathrm{0,1},2,$ (6-37)

and the chemical loading moments are

$B_m(t;\theta )={\int }_0^L{Y}_e\left(c\right(z,t);\theta )\beta \left[c\right(z,t)-c_{\mathrm{r}\mathrm{e}\mathrm{f}}]z^{m}dz,m=\mathrm{0,1}.$ (6-38)

Substituting Equations (6-33), (6-34), (6-37) and (6-38) into the force and moment balance Equations (6-35) and (6-36) gives

$\left[Y_{cu}h+A_0(t;\theta )\right]{\epsilon }_0\left(t\right)+\left[-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+A_1(t;\theta )\right]\kappa \left(t\right)=B_0(t;\theta ),$ (6-39)

$\left[-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+A_1(t;\theta )\right]{\epsilon }_0\left(t\right)+\left[{\displaystyle \frac{Y_{cu}{h}^3}{3}}+A_2(t;\theta )\right]\kappa \left(t\right)=B_1(t;\theta ).$ (6-40)

Equivalently, Equations (6-39) and (6-40) can be rewritten as

$\left[\begin{array}{cc}{K}_{00}(t;\theta )& K_{01}(t;\theta )\\ K_{01}(t;\theta )& K_{11}(t;\theta )\end{array}\right]\left[\begin{array}{c}{\epsilon }_0\left(t\right)\\ \kappa \left(t\right)\end{array}\right]=\left[\begin{array}{c}{B}_0(t;\theta )\\ B_1(t;\theta )\end{array}\right],$ (6-41)

where

$K_{00}(t;\theta )=Y_{cu}h+A_0(t;\theta ),$ (6-42)

$K_{01}(t;\theta )=-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+A_1(t;\theta ),$ (6-43)

$K_{11}(t;\theta )={\displaystyle \frac{Y_{cu}{h}^3}{3}}+A_2(t;\theta ).$ (6-44)

For a prescribed modulus-concentration relation $Y_e(c;\theta )$, the predicted curvature can be given by

${\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}(t;\theta )={\displaystyle \frac{K_{00}(t;\theta )B_1(t;\theta )-K_{01}(t;\theta )B_0(t;\theta )}{K_{00}(t;\theta )K_{11}(t;\theta )-K_{01}^2(t;\theta )}}.$ (6-45)

This equation defines the forward chemo-mechanical operator for galvanostatic cycling, i.e.,

$Y_e(c;\theta )\mapsto {\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}(t;\theta ).$ (6-46)

Once ${\epsilon }_0\left(t\right)$ and $\kappa \left(t\right)$ are obtained, the local active electrode layer stress follows from Equation (6-33). The average active electrode layer stress can be given by

${\overline{\sigma }}_{\parallel }^e\left(t\right)={\displaystyle \frac{1}{L}}{\int }_0^L{Y}_e\left(c\right(z,t);\theta )\left[{\epsilon }_0\left(t\right)+\kappa \left(t\right)z-\beta \left(c\right(z,t)-c_{\mathrm{r}\mathrm{e}\mathrm{f}})\right]dz,$ (6-47)

or equivalently,

${\overline{\sigma }}_{\parallel }^e\left(t\right)={\displaystyle \frac{1}{L}}\left[A_0(t;\theta ){\epsilon }_0\left(t\right)+A_1(t;\theta )\kappa \left(t\right)-B_0(t;\theta )\right].$ (6-48)

The maximum absolute active-layer stress can be given by

${\sigma }_{max}^e\left(t\right)=\underset{0\le z\le L}{\mathrm{m}\mathrm{a}\mathrm{x}}\mid Y_e\left(c\right(z,t);\theta )\left[{\epsilon }_0\left(t\right)+\kappa \left(t\right)z-\beta \left(c\right(z,t)-c_{\mathrm{r}\mathrm{e}\mathrm{f}})\right]\mid .$ (6-49)

6.4. Inverse Identification of a Common Modulus Function During Cycling

The goal of the present inverse problem is to identify a single concentration-dependent biaxial modulus relation $Y_e\left(c\right)$ from measured curvature histories during galvanostatic lithiation-delithiation or delithiation-lithiation cycling.

Unlike an independent lithiation or delithiation analysis, the present framework imposes

$Y_e^{\mathrm{i}\mathrm{n}}\left(c\right)=Y_e^{\mathrm{e}\mathrm{x}}\left(c\right)=Y_e\left(c\right).$ (6-50)

Therefore, all lithiation segments, all delithiation segments, all cycles, and all current-rate experiments are fitted using one common parameter vector $\theta$.

This assumption corresponds to reversible elastic behavior with no modulus hysteresis between lithium insertion and extraction. It does not require the curvature response itself to be identical during lithiation and delithiation at the same average concentration, because the through-thickness concentration profile may differ between the two directions. In cycling, different concentration gradients can produce different curvature responses even when the local material function $Y_e\left(c\right)$ is unique.

To identify the local modulus function during cycling, a two-step inverse identification strategy similar to that used in Section V is adopted.

Step I: Apparent modulus identification using an effective constant modulus for cycling

Although the primary objective of the present framework is to identify a local concentration-dependent modulus $Y_e\left(c\right)$, it is useful to introduce an apparent homogeneous electrode modulus as a preliminary diagnostic quantity. This apparent modulus provides a reduced-order description of the overall chemo-mechanical response and can be used to guide the choice of a parametric form for $Y_e\left(c\right)$. However, its interpretation in galvanostatic cycling is more subtle than in a single-step lithiation or delithiation experiment, because the curvature depends not only on the thickness-averaged concentration but also on the through-thickness concentration profile inherited from previous cycling segments.

In this apparent-modulus approximation, the spatially varying modulus is replaced by a time-dependent homogeneous modulus,

$Y_e\left(c\right(z,t\left)\right)\approx {\bar{Y}}_e\left(t\right),$ （6-51）

where ${\bar{Y}}_e\left(t\right)$ is not interpreted as the true local elastic modulus at concentration $c$, but as an effective homogeneous modulus that reproduces the measured curvature at time $t$ under the assumed concentration field $c(z,t)$. In the cycling setting, this quantity is denoted by

${\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}\left(t\right).$ (6-52)

The superscript “cyc” emphasizes that the apparent modulus is extracted from cycling data and therefore contains information from both the local concentration distribution and the cycling history.

Under the homogeneous-modulus approximation, the active-layer modulus is independent of $z$ at a fixed time. The stiffness moments in the active electrode become

$A_0\left(t\right)={\int }_0^L{\bar{Y}}_e\left(t\right)dz={\bar{Y}}_e\left(t\right)L,$ (6-53)

$A_1\left(t\right)={\int }_0^L{\bar{Y}}_e\left(t\right)zdz={\bar{Y}}_e\left(t\right){\displaystyle \frac{L^2}{2}},$ (6-54)

$A_2\left(t\right)={\int }_0^L{\bar{Y}}_e\left(t\right)z^{2}dz={\bar{Y}}_e\left(t\right){\displaystyle \frac{L^3}{3}}.$ (6-55)

The chemical loading moments reduce to

$B_0\left(t\right)={\bar{Y}}_e\left(t\right)\beta {\int }_0^L\left[c\right(z,t)-c_{\mathrm{r}\mathrm{e}\mathrm{f}}]dz,$ (6-56)

$B_1\left(t\right)={\bar{Y}}_e\left(t\right)\beta {\int }_0^L\left[c\right(z,t)-c_{\mathrm{r}\mathrm{e}\mathrm{f}}]zdz.$ (6-57)

For compactness, we define the concentration moments as

$G_0\left(t\right)=\beta {\int }_0^L\left[c\right(z,t)-c_{\mathrm{r}\mathrm{e}\mathrm{f}}]dz,$ (6-58)

$G_1\left(t\right)=\beta {\int }_0^L\left[c\right(z,t)-c_{\mathrm{r}\mathrm{e}\mathrm{f}}]zdz.$ (6-59)

Then, we have

$B_0\left(t\right)={\bar{Y}}_e\left(t\right)G_0\left(t\right),$ (6-60)

$B_1\left(t\right)={\bar{Y}}_e\left(t\right)G_1\left(t\right).$ (6-61)

Similarly, we define the geometric moments of the active layer as

$E_0=L,$ (6-62)

$E_1={\displaystyle \frac{L^2}{2}},$ (6-63)

$E_2={\displaystyle \frac{L^3}{3}}.$ (6-64)

The substrate stiffness moments are defines as

$S_0=Y_{cu}h,$ (6-65)

$S_1=-{\displaystyle \frac{Y_{cu}{h}^2}{2}},$ (6-66)

$S_2={\displaystyle \frac{Y_{cu}{h}^3}{3}}.$ (6-67)

With these definitions, the bilayer stiffness terms become

$K_{00}\left(t\right)=S_0+{\bar{Y}}_e\left(t\right)E_0,$ (6-68)

$K_{01}\left(t\right)=S_1+{\bar{Y}}_e\left(t\right)E_1,$ (6-69)

$K_{11}\left(t\right)=S_2+{\bar{Y}}_e\left(t\right)E_2.$ (6-70)

At each experimental time $t_k$, the measured curvature may be used to infer ${\bar{Y}}_e\left(t_k\right)$. Let the curvature used for apparent-modulus extraction be denoted by

${\kappa }_{\mathrm{u}\mathrm{s}\mathrm{e}}\left(t_k\right).$ (6-71)

If absolute curvature is reliable and the stress-free reference state is known, one may set

${\kappa }_{\mathrm{u}\mathrm{s}\mathrm{e}}\left(t_k\right)={\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(t_k\right).$ (6-72)

If only curvature increments are meaningful because of residual stress, initial curvature, or mounting offsets, one may instead use

${\kappa }_{\mathrm{u}\mathrm{s}\mathrm{e}}\left(t_k\right)=\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(t_k\right)={\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(t_k\right)-{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(0\right).$ (6-73)

However, when curvature increments are used, the extracted apparent modulus should be interpreted only as an incremental diagnostic quantity rather than an absolute material property.

The curvature predicted by the homogeneous-modulus approximation is

${\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\left(t\right)={\displaystyle \frac{K_{00}\left(t\right)B_1\left(t\right)-K_{01}\left(t\right)B_0\left(t\right)}{K_{00}\left(t\right)K_{11}\left(t\right)-K_{01}^2\left(t\right)}}.$ (6-74)

Substituting Equations (6-60)-(6-70) into Eq. (6-74) gives

${\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\left(t\right)={\displaystyle \frac{[S_0+{\bar{Y}}_e{E}_0]{\bar{Y}}_e{G}_1-[S_1+{\bar{Y}}_e{E}_1]{\bar{Y}}_e{G}_0}{{\left[S_0+{\bar{Y}}_e{E}_0\left]\right[S_2+{\bar{Y}}_e{E}_2]-[S_1+{\bar{Y}}_e{E}_1\right]}^2}}.$ (6-75)

At a fixed time, all quantities except ${\bar{Y}}_e$ are known from the substrate properties, film geometry, measured curvature, and computed concentration field. Therefore, ${\bar{Y}}_e\left(t\right)$ can be obtained by solving

${\kappa }_{\mathrm{u}\mathrm{s}\mathrm{e}}\left(t\right)={\displaystyle \frac{[S_0+{\bar{Y}}_e{E}_0]{\bar{Y}}_e{G}_1-[S_1+{\bar{Y}}_e{E}_1]{\bar{Y}}_e{G}_0}{{\left[S_0+{\bar{Y}}_e{E}_0\left]\right[S_2+{\bar{Y}}_e{E}_2]-[S_1+{\bar{Y}}_e{E}_1\right]}^2}}.$ (6-76)

Equation (6-76) can be rearranged into a quadratic equation for ${\bar{Y}}_e\left(t\right)$. For clarity, define

$Y={\bar{Y}}_e\left(t\right),\kappa ={\kappa }_{\mathrm{u}\mathrm{s}\mathrm{e}}\left(t\right),$ (6-77)

and suppress the explicit time dependence of $G_0$ and $G_1$. Then Eq. (6-76) becomes

$\kappa \left[(S_0+Y{E}_0)(S_2+Y{E}_2)-(S_1+Y{E}_1{)}^2\right]=Y\left[(S_0+Y{E}_0)G_1-(S_1+Y{E}_1)G_0\right].$ (6-78)

Expanding both sides yields

$\mathcal{A}\left(t\right)Y^2+\mathcal{B}\left(t\right)Y+\mathcal{C}\left(t\right)=0$, (6-79)

where

$\mathcal{A}\left(t\right)=\kappa \left(E_0{E}_2−E_1^2\right)-\left(E_0{G}_1−E_1{G}_0\right),$ (6-80)

$\mathcal{B}\left(t\right)=\kappa \left(S_0{E}_2+E_0{S}_2-2S_1{E}_1\right)-\left(S_0{G}_1−S_1{G}_0\right),$ (6-81)

and

$\mathcal{C}\left(t\right)=\kappa \left(S_0{S}_2−S_1^2\right).$ (6-82)

Thus, the apparent homogeneous modulus at time $t$ is obtained from

${\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}\left(t\right)={\displaystyle \frac{-\mathcal{B}\left(t\right)\pm \sqrt{{\mathcal{B}}^2\left(t\right)-4\mathcal{A}\left(t\right)\mathcal{C}\left(t\right)}}{2\mathcal{A}\left(t\right)}}.$ (6-83)

Among the two mathematical roots, the physically admissible one should be selected by requiring

${\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}\left(t\right)>0.$ (6-84)

If both roots are positive, the root that lies within a physically reasonable modulus range should be chosen. If no positive real root exists, then the homogeneous-modulus approximation is inconsistent with the measured curvature, the assumed concentration field, or the adopted reference state at that time point. In practice, such points should be excluded from apparent-modulus plotting or treated as indicators of measurement noise, model mismatch, or insufficient identifiability.

When $\mathcal{A}\left(t\right)$ is close to zero, Equation (6-79) becomes nearly linear. In that case, the apparent modulus may be obtained from

${\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}\left(t\right)\approx -{\displaystyle \frac{\mathcal{C}\left(t\right)}{\mathcal{B}\left(t\right)}},\mid \mathcal{A}\left(t\right)\mid \ll 1.$ (6-85)

The discriminant

$\mathsf{\Delta }\left(t\right)={\mathcal{B}}^2\left(t\right)-4\mathcal{A}\left(t\right)\mathcal{C}\left(t\right)$ (6-86)

also provides a useful diagnostic. A negative discriminant indicates that, at that time point, no real homogeneous modulus can exactly reproduce the measured curvature under the assumed concentration field and reference strain. This may occur when experimental noise is large, when the curvature baseline is uncertain, when $D$, $\beta$, or $c_{\mathrm{r}\mathrm{e}\mathrm{f}}$ is inaccurate, or when the true electrode response cannot be represented by a homogeneous elastic modulus.

After evaluating ${\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}\left(t\right)$ at all admissible time points, the apparent modulus may be plotted against the average concentration,

$\bar{c}\left(t\right)={\displaystyle \frac{1}{L}}{\int }_0^{L}c(z,t)dz.$ (6-87)

This produces the apparent relation

${\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}={\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}\left(\bar{c}\right).$ (6-88)

Such a plot can be useful for identifying broad trends. For example, it may suggest whether the modulus decreases, increases, saturates, or varies approximately linearly with lithium content. It may also help determine whether a linear, exponential, saturating exponential, or spline representation should be used for the full local modulus function $Y_e\left(c\right)$.

However, the apparent relation in Eq. (6-88) should not be interpreted as the true material function $Y_e\left(c\right)$. In cycling, the same average concentration may correspond to different local concentration profiles. For example, during lithiation, a given $\bar{c}$ may be associated with a higher lithium concentration near the active surface and a lower concentration near the substrate. During delithiation, the same $\bar{c}$ may occur with the opposite type of concentration gradient. Therefore,

$\bar{c}\left(t_a\right)=\bar{c}\left(t_b\right)$ (6-89)

does not generally imply

$c(z,t_a)=c(z,t_b).$ (6-90)

Consequently, even if the true local modulus is a single reversible function $Y_e\left(c\right)$, the apparent modulus may exhibit different lithiation and delithiation branches when plotted against $\bar{c}$. This apparent hysteresis can arise solely from diffusion-induced concentration gradients and does not necessarily imply intrinsic modulus hysteresis.

This distinction is particularly important in high-rate cycling. At low current density or after sufficiently long relaxation, the concentration field becomes nearly uniform:

$c(z,t)\approx \bar{c}\left(t\right).$ (6-91)

In this limiting case, ${\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}\left(t\right)$ is more closely related to the local modulus $Y_e\left(\bar{c}\right)$. By contrast, at high current density or for short half-cycles, the concentration profile may be strongly nonuniform, and the apparent modulus becomes a weighted effective quantity influenced by the spatial distribution of lithium, the substrate constraint, and the bending moment generated by the chemical strain gradient.

Therefore, ${\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}\left(t\right)$ should be regarded as a diagnostic or visualization tool rather than the final identified modulus. It can be used to provide an initial estimate of the modulus magnitude, to reveal approximate softening or hardening trends with lithium content, to guide the choice of a parametric form for $Y_e(c;\theta )$, to identify time regions where the data are inconsistent with a simple homogeneous elastic response, and to compare different cycling rates qualitatively before performing the full inversion.

Step II: Parameterized full inversion of${\mathit{Y}}_{\mathit{e}}\left(\mathit{c}\right)$

The full inverse identification should be based on the local variable-modulus model,

$Y_e=Y_e(c;\theta ),$ (6-92)

using the complete through-thickness concentration field $c(z,t)$ generated by the entire cycling history. In the final inversion, all lithiation and delithiation segments should be fitted simultaneously using a single parameter vector $\theta$, rather than using ${\bar{Y}}_e^{\mathrm{c}\mathrm{y}\mathrm{c}}\left(\bar{c}\right)$ as the material law directly.

First, we need to define the concentration interval for fitting. The modulus function is identified only over the concentration interval sampled by the experiments. Since the local concentration may differ from the average concentration, the fitting interval must be determined from the local concentration fields over all experiments, all cycles, all times, and all thickness positions.

Suppose there are $Q$ cycling experiments. For the $q$-th experiment, the imposed flux history is $J^{\left(q\right)}\left(t\right)$, and the resulting concentration field is

$c^{\left(q\right)}(z,t).$ (6-93)

The minimum and maximum local concentrations sampled by all cycling experiments are

$c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}=\underset{q}{\mathrm{m}\mathrm{i}\mathrm{n}}\underset{t}{\mathrm{m}\mathrm{i}\mathrm{n}}\underset{0\le z\le L}{\mathrm{m}\mathrm{i}\mathrm{n}}{c}^{\left(q\right)}(z,t),$ (6-94)

$c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}=\underset{q}{\mathrm{m}\mathrm{a}\mathrm{x}}\underset{t}{\mathrm{m}\mathrm{a}\mathrm{x}}\underset{0\le z\le L}{\mathrm{m}\mathrm{a}\mathrm{x}}{c}^{\left(q\right)}(z,t).$ (6-95)

To avoid numerical extrapolation caused by discretization, series truncation, or quadrature error, a small safety margin may be introduced:

$c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}=c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}-{\delta }_c,$ (6-96)

$c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}=c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}+{\delta }_c,$ (6-97)

where ${\delta }_c\ge 0$ is a small concentration margin, which can be set as

${\delta }_c={\displaystyle \frac{\alpha -1}{2}}\left(c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}−c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}\right),$ (6-98)

where the multiplicative safety factor $\alpha \ge 1$.

A unified normalized concentration variable is then defined as

$s={\displaystyle \frac{c-c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}}{c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}-c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}}}.$ (6-99)

This same normalized variable is used for both lithiation and delithiation portions of the cycling experiment.

Second, we need to parameterize the modulus function. The local biaxial modulus is represented by a finite-dimensional function

$Y_e\left(c\right)=Y_e(c;\theta ).$ (6-100)

The choice of parameterization may be guided by the apparent modulus trends, prior mechanical measurements, or model comparison.

A linear model may be written as

$Y_e(c;\theta )={\theta }_0+{\theta }_{1}c,$ (6-101)

where

$\theta =\left[\begin{array}{c}{\theta }_0\\ {\theta }_1\end{array}\right].$ (6-102)

A quadratic exponential model may be written as

$Y_e(c;\theta )=\mathrm{e}\mathrm{x}\mathrm{p}({\eta }_0+{\theta }_{1}s+{\theta }_2{s}^2),$ (6-103)

where

$\theta =\left[\begin{array}{c}{\eta }_0\\ {\theta }_1\\ {\theta }_2\end{array}\right].$ (6-104)

This form automatically guarantees positivity.

If the modulus changes rapidly over part of the concentration range and then approaches a plateau, a saturating exponential model may be used:

$Y_e(c;\theta )=Y_{\mathrm{\infty }}+(Y_0-Y_{\mathrm{\infty }})\mathrm{e}\mathrm{x}\mathrm{p}(-ks),$ (6-105)

where $Y_0$ is the modulus at $s=0$, $Y_{\mathrm{\infty }}$ is the saturated modulus at large $s$, and $k$ is a transition-rate parameter.

To guarantee parameter positivity, a logarithmic parameterization can be introduced:

$Y_0=\mathrm{e}\mathrm{x}\mathrm{p}\left(a\right),Y_{\mathrm{\infty }}=\mathrm{e}\mathrm{x}\mathrm{p}\left(b\right),k=\mathrm{e}\mathrm{x}\mathrm{p}\left(g\right).$ (6-106)

Then, we have

$Y_e(c;\theta )=\mathrm{e}\mathrm{x}\mathrm{p}\left(b\right)+\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(a\right)-\mathrm{e}\mathrm{x}\mathrm{p}\left(b\right)\right]\mathrm{e}\mathrm{x}\mathrm{p}[-\mathrm{e}\mathrm{x}\mathrm{p}(g\left)s\right],$ (6-107)

with

$\theta =\left[\begin{array}{c}a\\ b\\ g\end{array}\right].$ (6-108)

More flexible representations, such as B-splines or monotonic splines, may also be used if sufficient cycling data are available.

After choosing the modulus function, the forward calculation is conducted. For a given parameter vector $\theta$, the forward calculation proceeds as follows. First, for each cycling experiment $q$, the concentration field $c^{\left(q\right)}(z,t)$ is calculated from the imposed flux history $J^{\left(q\right)}\left(t\right)$. For a piecewise-constant galvanostatic protocol, the recursive solution in Equations (6-23)-(6-26) is used segment by segment.

Second, at each measurement time $t_k$, the stiffness moments are evaluated by

$A_m^{\left(q\right)}(t_k;\theta )={\int }_0^L{Y}_e\left(c^{\left(q\right)}\right(z,t_k);\theta )z^{m}dz,m=\mathrm{0,1},2.$ (6-109)

The chemical loading moments are evaluated as

$B_m^{\left(q\right)}(t_k;\theta )={\int }_0^L{Y}_e\left(c^{\left(q\right)}\right(z,t_k);\theta )\beta \left[c^{\left(q\right)}\right(z,t_k)-c_{\mathrm{r}\mathrm{e}\mathrm{f}}]z^{m}dz,m=\mathrm{0,1}.$ (6-110)

The stiffness matrix terms are given by

$K_{00}^{\left(q\right)}(t_k;\theta )=Y_{cu}h+A_0^{\left(q\right)}(t_k;\theta ),$ (6-111)

$K_{01}^{\left(q\right)}(t_k;\theta )=-{\displaystyle \frac{Y_{cu}{h}^2}{2}}+A_1^{\left(q\right)}(t_k;\theta ),$ (6-112)

$K_{11}^{\left(q\right)}(t_k;\theta )={\displaystyle \frac{Y_{cu}{h}^3}{3}}+A_2^{\left(q\right)}(t_k;\theta ).$ (6-113)

The predicted curvature is then given by

${\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{\left(q\right)}(t_k;\theta )={\displaystyle \frac{K_{00}^{\left(q\right)}(t_k;\theta )B_1^{\left(q\right)}(t_k;\theta )-K_{01}^{\left(q\right)}(t_k;\theta )B_0^{\left(q\right)}(t_k;\theta )}{K_{00}^{\left(q\right)}(t_k;\theta )K_{11}^{\left(q\right)}(t_k;\theta )-\left[K_{01}^{\left(q\right)}\right(t_k;\theta ){]}^2}}.$ (6-114)

To reduce the influence of initial curvature, residual stress, and mounting error, curvature increments may be used:

$\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{\left(q\right)}\left(t_k\right)={\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{\left(q\right)}\left(t_k\right)-{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{\left(q\right)}\left(0\right),$ (6-115)

$\mathsf{\Delta }{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{\left(q\right)}(t_k;\theta )={\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{\left(q\right)}(t_k;\theta )-{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{\left(q\right)}(0;\theta ).$ (6-116)

If significant baseline drift occurs during long multi-cycle experiments, curvature increments may instead be defined relative to the beginning of each cycle or half-cycle. Let $t_{qm}^{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}$ denote the starting time of the $m$-th cycle segment or half-cycle in experiment $q$. Then

$\mathsf{\Delta }{\kappa }_{exp}^{\left(q,m\right)}\left(t_k\right)={\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{\left(q\right)}\left(t_k\right)-{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{\left(q\right)}\left(t_{qm}^{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}\right),$ (6-117)

$\mathsf{\Delta }{\kappa }_{pred}^{\left(q,m\right)}(t_k;\theta )={\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{\left(q\right)}(t_k;\theta )-{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{\left(q\right)}(t_{qm}^{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}};\theta ).$ (6-118)

Using segment-wise increments emphasizes the incremental chemo-mechanical response within each cycle or half-cycle, whereas using the global initial increment preserves more information about the accumulated stress evolution.

The inverse problem is formulated as a single joint optimization problem:

$\underset{\theta }{\mathrm{m}\mathrm{i}\mathrm{n}}{\mathsf{\Phi }}_{\mathrm{c}\mathrm{y}\mathrm{c}}\left(\theta \right).$ (6-119)

For $Q$ cycling experiments, the objective function may be written as

${\mathsf{\Phi }}_{\mathrm{c}\mathrm{y}\mathrm{c}}\left(\theta \right)={\displaystyle \sum _{q=1}^Q}{\displaystyle \sum _{k=1}^{N_q}}{w}_{qk}{\left[\mathsf{\Delta }{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{\left(q\right)}(t_k;\theta )-\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{\left(q\right)}\left(t_k\right)\right]}^2+{\mathsf{\Phi }}_{\mathrm{r}\mathrm{e}\mathrm{g}}\left(\theta \right),$ (6-120)

where $N_q$ is the number of curvature data points in experiment $q$, $w_{qk}$ are weights, and ${\mathsf{\Phi }}_{\mathrm{r}\mathrm{e}\mathrm{g}}$ is a regularization term.

A useful weighting choice is

$w_{qk}={\displaystyle \frac{1}{N_q[{\kappa }_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}}^{\left(q\right)}{]}^2}},$ (6-121)

with

${\kappa }_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}}^{\left(q\right)}=\underset{k}{\mathrm{m}\mathrm{a}\mathrm{x}}\mid \mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{\left(q\right)}\left(t_k\right)\mid .$ (6-122)

This weighting balances experiments with different numbers of data points and different curvature magnitudes.

For multi-cycle experiments, one may also balance the contributions from individual cycles or half-cycles. If $m=1,\dots ,M_q$ indexes the cycle segments or half-cycles in experiment $q$, the objective function may be written as

${\mathsf{\Phi }}_{\mathrm{c}\mathrm{y}\mathrm{c}}\left(\theta \right)={\displaystyle \sum _{q=1}^Q}{\displaystyle \sum _{m=1}^{M_q}}{\displaystyle \sum _{k=1}^{N_{qm}}}{w}_{qmk}{\left[\mathsf{\Delta }{\kappa }_{pred}^{\left(q,m\right)}(t_k;\theta )-\mathsf{\Delta }{\kappa }_{exp}^{\left(q,m\right)}\left(t_k\right)\right]}^2+{\mathsf{\Phi }}_{\mathrm{r}\mathrm{e}\mathrm{g}}\left(\theta \right).$ (6-123)

A corresponding segment-wise weight is

$w_{qmk}={\displaystyle \frac{1}{N_{qm}[{\kappa }_{scale}^{\left(q,m\right)}{]}^2}},$ (6-124)

where

${\kappa }_{scale}^{\left(q,m\right)}=\underset{k}{\mathrm{m}\mathrm{a}\mathrm{x}}\mid {∆\kappa }_{exp}^{\left(q,m\right)}\left(t_k\right)\mid .$ (6-125)

For low-dimensional models, such as the linear, quadratic exponential, or saturating exponential forms, regularization may be omitted. For spline-based representations, regularization is generally required to suppress nonphysical oscillations. A typical smoothness penalty is

${\mathsf{\Phi }}_{\mathrm{r}\mathrm{e}\mathrm{g}}\left(\theta \right)=\lambda {\int }_{c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}}^{c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}}{\left[{\displaystyle \frac{d^2{Y}_e(c;\theta )}{d{c}^2}}\right]}^{2}dc,$ (6-126)

where $\lambda$ controls the strength of regularization.

During the optimization process, several physical constraints should be applied. The modulus must remain positive over the entire fitting interval, i.e.,

$Y_e(c;\theta )\ge Y_{\mathrm{m}\mathrm{i}\mathrm{n}}>0,c\in [c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}},c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}].$ （6-127）

Additional bounds may also be imposed:

$Y_{\mathrm{m}\mathrm{i}\mathrm{n}}\le Y_e(c;\theta )\le Y_{\mathrm{m}\mathrm{a}\mathrm{x}},c\in [c_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}},c_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}].$ (6-128)

For the linear model,

$Y_e(c;\theta )={\theta }_0+{\theta }_{1}c,$ (6-129)

the positivity condition reduces to endpoint constraints:

${\theta }_0+{\theta }_1{c}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{f}\mathrm{i}\mathrm{t}}\ge Y_{\mathrm{m}\mathrm{i}\mathrm{n}},$ (6-130)

${\theta }_0+{\theta }_1{c}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{f}\mathrm{i}\mathrm{t}}\ge Y_{\mathrm{m}\mathrm{i}\mathrm{n}}.$ (6-131)

If the active electrode is known to soften with increasing lithium concentration, then

${\displaystyle \frac{d{Y}_e}{dc}}\le 0.$ (6-132)

For the linear model, this gives

${\theta }_1\le 0.$ (6-133)

If the electrode is known to harden with increasing lithium concentration, then

${\displaystyle \frac{d{Y}_e}{dc}}\ge 0,$ (6-134)

which gives

${\theta }_1\ge 0.$ (6-135)

For the quadratic exponential model,

$Y_e(c;\theta )=\mathrm{e}\mathrm{x}\mathrm{p}({\eta }_0+{\theta }_{1}s+{\theta }_2{s}^2),$ (6-136)

positivity is automatic. Since

${\displaystyle \frac{d{Y}_e}{ds}}=Y_e\left(s\right)({\theta }_1+2{\theta }_{2}s),$ (6-137)

monotonic softening with increasing lithium concentration requires

${\theta }_1\le 0,{\theta }_1+2{\theta }_2\le 0.$ (6-138)

Monotonic hardening requires

${\theta }_1\ge 0,{\theta }_1+2{\theta }_2\ge 0.$ (6-139)

For the saturating exponential model,

$Y_e(c;\theta )=Y_{\mathrm{\infty }}+(Y_0-Y_{\mathrm{\infty }})\mathrm{e}\mathrm{x}\mathrm{p}(-ks),$ (6-140)

the derivative with respect to $s$ is

${\displaystyle \frac{d{Y}_e}{ds}}=-k(Y_0-Y_{\mathrm{\infty }})\mathrm{e}\mathrm{x}\mathrm{p}(-ks).$ (6-141)

Therefore, if

$Y_0>Y_{\mathrm{\infty }}>0,k>0,$ (6-142)

Then

${\displaystyle \frac{d{Y}_e}{dc}}\le 0,$ (6-143)

and the modulus decreases monotonically with increasing lithium concentration. Under logarithmic parameterization, this softening condition is imposed as

$a>b.$ (6-144)

If instead the electrode hardens with increasing lithium concentration, then

$Y_{\mathrm{\infty }}>Y_0>0,k>0,$ (6-145)

which corresponds to

$b>a.$ (6-146)

The monotonicity condition is always imposed with respect to lithium concentration $c$, not with respect to time. Thus, if the electrode softens with increasing lithium concentration, the modulus decreases during lithiation and increases during delithiation, even though the same function $Y_e\left(c\right)$ is used in both directions.

After optimization, the residual for the $q$-th experiment is defined as

$r_{qk}=\mathsf{\Delta }{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{\left(q\right)}(t_k;\widehat{\theta })-\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{\left(q\right)}\left(t_k\right),$ (6-147)

where $\widehat{\theta }$ is the optimal parameter vector.

The root-mean-square error is

${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}^{\left(q\right)}=\sqrt{{\displaystyle \frac{1}{N_q}}{\displaystyle \sum _{k=1}^{N_q}}{r}_{qk}^2}.$ (6-148)

The relative root-mean-square error is

${\mathrm{R}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}^{\left(q\right)}=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{k=1}^{N_q}}{r}_{qk}^2}{{\displaystyle \sum _{k=1}^{N_q}}[\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{\left(q\right)}\left(t_k\right){]}^2}}}.$ (6-149)

For segment-wise or half-cycle-wise analysis,

$r_{qmk}=\mathsf{\Delta }{\kappa }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(q,m)}(t_k;\widehat{\theta })-\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,m)}\left(t_k\right),$ (6-150)

${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}^{(q,m)}=\sqrt{{\displaystyle \frac{1}{N_{qm}}}{\displaystyle \sum _{k=1}^{N_{qm}}}{r}_{qmk}^2},$ (6-151)

and

${\mathrm{R}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}^{(q,m)}=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{k=1}^{N_{qm}}}{r}_{qmk}^2}{{\displaystyle \sum _{k=1}^{N_{qm}}}[\mathsf{\Delta }{\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}^{(q,m)}\left(t_k\right){]}^2}}}.$ (6-152)

A suitable common modulus model should reproduce curvature histories over lithiation segments, delithiation segments, and multiple cycles with small and randomly distributed residuals. Systematic residuals that differ between lithiation and delithiation may indicate that the assumption

$Y_e^{\mathrm{i}\mathrm{n}}\left(c\right)=Y_e^{\mathrm{e}\mathrm{x}}\left(c\right)$ (6-153)

is inadequate, or that additional physics such as damage, plasticity, hysteresis, stress-dependent diffusion, side reactions, or evolving interfacial conditions must be included.

After the common local modulus relation $Y_e(c;\widehat{\theta })$ has been identified, the stress field is calculated using Equation (6-33). The average active-layer stress and maximum absolute active-layer stress are then calculated using Equations (6-47) and (6-49), respectively.

If curvature increments are fitted rather than absolute curvatures, the resulting stresses should be interpreted relative to the chosen reference state. If absolute stresses are required, the chemical strain should be referred to a true stress-free concentration $c_{\mathrm{r}\mathrm{e}\mathrm{f}}$, and any initial residual stress should be included or independently characterized.