Data-driven topology optimization of all-solid-
state batteries considering conductive additive
material informed by microstructure analysis

Naoyuki Ishida¹*, Kozo Furuta²†, Masashi Kishimoto¹†,
Takamitsu Sasaki¹†, Hiroshi Iwai¹†, Kazuhiro Izui²†,
Shinji Nishiwaki¹†

¹*Department of Mechanical Engineering and Science, Kyoto University,
Kyotodaigaku-katsura, Kyoto, 615-8246, Kyoto, Japan.

²Department of Micro Engineering, Kyoto University,
Kyotodaigaku-katsura, Kyoto, 615-8246, Kyoto, Japan.

*Corresponding author(s). E-mail(s):

ishida.naoyuki.23m@st.kyoto-u.ac.jp;

Contributing authors: furuta.kozo.5n@kyoto-u.ac.jp;

kishimoto.masashi.3m@kyoto-u.ac.jp;

sasaki.takamitsu.36c@st.kyoto-u.ac.jp; iwai.hiroshi.4x@kyoto-u.ac.jp;

izui@me.kyoto-u.ac.jp; shinji@prec.kyoto-u.ac.jp;

†These authors contributed equally to this work.

Abstract

All-solid-state batteries (ASSBs) have emerged as a promising alternative to con-
ventional lithium-ion batteries, offering enhanced safety and efficiency through
the utilization of solid electrolytes, which simultaneously improve effective energy
density and mitigate ignition risks. However, these batteries exhibit reduced con-
ductivity and electrochemical reaction rates compared to conventional lithium-ion
batteries due to insufficient contact area between constituent materials. Previ-
ous experimental analyses have demonstrated that battery performance can be
enhanced through the implementation of functionally graded materials (FGM),
which introduce gradients in the volume fractions of mixed materials within
the battery structure. Nevertheless, existing FGM approaches have been con-
strained to one-dimensional design modifications and empirical trial-and-error
methodologies. To address these limitations, this paper presents a novel topology
optimization (TO) method for determining the optimal spatial distribution of
material volume fractions, specifically electrolytes, active materials, and conduc-
tive additives, within the composite anode of ASSBs. In contrast to conventional
approaches that rely on the Bruggeman model, we propose a physiochemically
rigorous material interpolation function derived through microstructural analysis
and least-squares approximation. This methodology enables the precise evalu-
ation of effective properties at the appropriate length scales while accounting
for the inherent heterogeneity of composite anodes. We formulate the optimiza-
tion problem as the minimization of overpotential during discharge processes,
with sensitivities derived using continuous Lagrange adjoint methods. Subse-
quently, we develop a comprehensive optimization algorithm. Through numerical
examples, we demonstrate that the optimized material distribution obtained via
the proposed method achieves an 11% enhancement in capacity compared to
conventional uniform composite distributions.

**Keywords:** Topology optimization, All-solid-state batteries, Sensitivity analysis, Time-dependent problem, Material interpolation function, Bruggeman model

# 1 Introduction

All-solid-state batteries (ASSBs) have gained interest as a notable energy storage
device. Compared to conventional liquid-based lithium-ion batteries, ASSBs are less
flammable and have a lower risk of ignition or explosion, resulting in improved safety
(Janek and Zeier, 2016). Furthermore, the use of solid electrolytes instead of liquid
electrolytes increases effective energy density and enables more efficient energy stor-
age (Wang et al., 2020; Zhang et al., 2020; Nanda et al., 2021). This technological
advantage will allow for more effective use of renewable energy and increase the cruis-
ing range of electric vehicles, and as a result, this contributes to the achievement of a
carbon-neutral society (Cano et al., 2018). Therefore, the development and widespread
use of ASSBs is a key approach for laying the foundation for a sustainable energy
system.

However, ASSBs have a drawback compared to liquid-based lithium-ion batteries, i.e., the insufficient contact area among materials composed in the composite region, resulting in a decrease of conductivity, diffusion coefficient, and electrochemical reaction rate (Ren et al., 2023; Deysher et al., 2022; Minnmann et al., 2021). To mitigate this issue, many research groups are making significant efforts in various ways, such as material development (Yong-Gun et al., 2020) and the improvement of manufacturing processes (Giraldo et al., 2021). Recently, conductive additive materials have been increasingly incorporated into the components of ASSBs to enhance electronic conductivity (Giraldo et al., 2021; Reisacher et al., 2023). As mentioned above, performance improvements and material development have been carried out, but due to the complexity of coupled reaction-transport phenomena in the electrodes, a design method for the batteries has yet to be established.

In such a battery development situation, optimization methods have the potential
to improve the performance of all kinds of batteries. A typical application is the
material composition optimization of battery cathode and anode materials to enhance
performances including effective energy density, cycle life, and charge/discharge rate (Lv et al., 2022; Ong, 2019; Shahzad et al., 2024). Those researchers used a combined method of high-throughput experiments and numerical approaches to efficiently find optimized solutions in a short period.

Another application of the optimization method is obtaining the optimized con-
figuration of battery devices such as the particle sizes of the active material and solid
electrolyte. Many research groups have attempted to find optimal values for the prop-
erties of the materials present in the battery by combining experiments and numerical
simulations to increase the contact area among materials. For example, Shi et al.
(2020) performed parametric optimization for the battery performance by varying the
particle sizes of the active material and solid electrolyte in the cathode based on a
numerical model confirmed by experimental data. They found that an optimal parti-
cle size ratio improved the battery performance. Clausnitzer et al. (2023) developed
a method to find the optimal solution based on structure-resolved simulations tar-
geting the values of the active material fraction, the density after sintering, and the
particle size as objective functions. Moreover, they revealed the effective energy den-
sity could be improved using the functionally graded material (FGM) approach (Deng
et al., 2024), wherein the volume fraction of active material is increased towards the
current collector, and the volume fraction of electrolyte is increased towards the sep-
arator. That is, the increase of the contact area and appropriate material distribution
by these proposed methods can increase the transport of both Li ions and electrons
and charge transfer, resulting in improving power density of the battery.

However, although previously proposed optimization methods provide practical
solutions, they are limited in terms of the design variables' degree of freedom and the
range of feasible solutions. Thus, topology optimization (TO) is one of most promising
approaches to obtain further improved and innovative design configurations.

Development of a TO method originated from structural problems (Bendsøe and Kikuchi, 1988), and it is now extended to various physics problems. Borrvall and Petersson (2003) and Matsumori et al. (2013) have applied it to fluid and thermal fluid problems, respectively, and they obtained the optimized fluid channel layout of heat exchangers. Wadbro and Berggren (2006) and Lu et al. (2013) have constructed a TO method for the acoustic problem and obtained an optimized acoustic horn and metasurface design, respectively.

TO has also been applied to various types of electrochemical devices, including
batteries and fuel cells. Iwai et al. (2011) constructed the TO method to enhance the
efficiency of solid oxide fuel cells (SOFC). Yaji et al. (2018) applied TO to improve
the performance of redox flow batteries, which utilize liquid electrolytes at both the
anode and cathode. The TO methods were also constructed by Lin et al. (2022) and
Roy et al. (2022) to refine internal flow paths within the redox flow batteries, further
increasing operational efficiency. Additionally, Mo et al. (2021) and Guo et al. (2022)
proposed the TO method that can provide cooling plates design, addressing the ther-
mal management challenges inherent to battery systems. Pejman and Najafi (2023)
presented a TO method for providing the optimized configurations of solid and compli-
ant phases of structural battery composites (SBC) by modeling them with equivalent
physical networks. SBC embodies a novel class of multifunctional materials, integrat-
ing the load-bearing capabilities of carbon fiber composites with the energy storage
properties characteristic of lithium-ion batteries (Snyder et al., 2007, 2008; Wong et al.,
2007; Liu et al., 2009). Furthermore, Alizadeh et al. (2023) proposed a TO method for
optimizing porous reactors, which concerns the spatial distribution of active and elec-
trolyte material volume fractions. As for ASSBs, Ishida et al. (2024) proposed a TO
method that uses a homogenization method to optimize the material distributions of
active and electrolyte materials within a unit cell defined in a composite anode.

In this paper, we propose a TO method that provides the optimized macroscopic material distribution configuration for ASSBs, i.e., the distributions of the volume fractions of electrolytes, active materials, and conductive additives in the composite anode. While conventional topology optimization typically employs the solid isotropic material with penalization (SIMP) method to penalize intermediate densities (Bendsøe, 1989; Bendsøe and Sigmund, 1999), this approach is not suitable for battery optimization where intermediate material compositions play a crucial role. Instead, we require material interpolation functions that accurately reflect the physicochemical properties of mixed material states, as these intermediate compositions are fundamental to battery performance.

The Bruggeman model has been widely used in battery analysis to predict effective transport coefficients based on material volume fractions (Chung et al., 2013). This model interpolates the material’s macroscopic properties with a function corresponding to the penalty parameter *p* of 1.5 in the SIMP method, assuming that the material particles are spherical (Tjaden et al., 2016). The unit cell size used in the Bruggeman model, in which transport is evaluated, is assumed to be sufficiently large relative to the particle size, i.e., approximately the size of the composite anode region. However, in this paper, the material distributions within the composite anode vary on a scale smaller than the size of the composite anode region, and this fact indicates the use of the Bruggeman model is not appropriate in the proposed optimization method. In other words, the size of the unit cell should be significantly smaller than the composite anode region.

Thus, this paper develops the material interpolation function with microstructure
analysis, i.e., that randomly distributes materials in a pseudo unit cell, evaluates their
properties, and uses the least-squares method to obtain an approximated function
from the obtained datasets. Next, the optimization problem is formulated to minimize
the overvoltage during the discharge process, which is one of the measures of energy
loss in a battery system. The sensitivity is based on the sensitivity derivation proposed
by Ishida et al. (2024). After that, the optimization algorithm is developed. Finally,
numerical examples demonstrate the effectiveness of the proposed method.

The rest of paper is organized as follows. Section 2 describes the mathematical formulation of governing equations for ASSBs and the optimization problem for minimizing the overvoltage. Section 3 presents the approximation method for the material interpolation functions. In Sect. 4 the optimization algorithm is detailed. After that, Sect. 5 provides several numerical examples to confirm the validity of the proposed method. Finally, the conclusions of this paper are presented in Sect. 6.

The variables used in this paper are summarized in Table 1.
<table><caption>**Table 1:** List of variables used in the paper</caption><thead><tr><th>Variable</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>ε<sub>a,n</sub></td><td>Volume fraction of active material in the composite anode</td><td>—</td></tr><tr><td>ε<sub>c</sub></td><td>Volume fraction of conductive additive in the composite anode</td><td>—</td></tr><tr><td>ε<sub>e,n</sub></td><td>Volume fraction of electrolyte material in the composite anode</td><td>—</td></tr><tr><td>c<sub>n</sub></td><td>Li concentration in the composite anode</td><td>mol/m<sup>3</sup></td></tr><tr><td>c<sub>p</sub></td><td>Li concentration in the composite cathode</td><td>mol/m<sup>3</sup></td></tr><tr><td>φ<sub>n</sub></td><td>Electrical potential of active material in the composite anode</td><td>V</td></tr><tr><td>φ<sub>p</sub></td><td>Electrical potential of active material in the composite cathode</td><td>V</td></tr><tr><td>φ<sub>e</sub></td><td>Electrical potential of electrolyte material</td><td>V</td></tr><tr><td>D<sup>eff</sup><sub>n</sub></td><td>Li diffusion coefficient of active material in the composite anode</td><td>m<sup>2</sup>/s</td></tr><tr><td>D<sup>eff</sup><sub>p</sub></td><td>Li diffusion coefficient of active material in the composite cathode</td><td>m<sup>2</sup>/s</td></tr><tr><td>σ<sup>eff</sup><sub>n</sub></td><td>Electronic conductivity in the composite anode</td><td>S/m</td></tr><tr><td>σ<sup>eff</sup><sub>p</sub></td><td>Electronic conductivity in the composite cathode</td><td>S/m</td></tr><tr><td>κ<sup>eff</sup></td><td>Li-ion conductivity</td><td>S/m</td></tr><tr><td>j<sup>Li</sup></td><td>Charge-transfer current density</td><td>A/m<sup>3</sup></td></tr><tr><td>η</td><td>Overvoltage</td><td>V</td></tr><tr><td>a<sup>eff</sup><sub>s</sub></td><td>Active surface area per electrode unit volume for charge-transfer reaction</td><td>m<sup>2</sup>/m<sup>3</sup></td></tr><tr><td>I</td><td>Current</td><td>A</td></tr><tr><td>Q<sub>cell</sub></td><td>Capacity</td><td>C</td></tr><tr><td>J</td><td>Objective function</td><td>V</td></tr><tr><td>a<sub>1</sub> - a<sub>2</sub></td><td>Adjoint variables</td><td>V · m/mol</td></tr><tr><td>a<sub>3</sub> - a<sub>5</sub></td><td>Adjoint variables</td><td>m/(S · s)</td></tr><tr><td>J'<sub>1</sub></td><td>Sensitivity for ε<sub>a,n</sub></td><td>V</td></tr><tr><td>J'<sub>2</sub></td><td>Sensitivity for ε<sub>c</sub></td><td>V</td></tr><tr><td>ρ</td><td>Randomly generated material distribution in the unit cell [0 : 1]</td><td>—</td></tr><tr><td>C<sub>1</sub> - C<sub>3</sub></td><td>Dimensionless scalar fields defined in the unit cell</td><td>—</td></tr><tr><td>q<sub>1</sub> - q<sub>3</sub></td><td>Dimensionless fluxes of C<sub>1</sub>, C<sub>2</sub>, and C<sub>3</sub></td><td>—</td></tr><tr><td>K</td><td>Dimensionless transport coefficient</td><td>—</td></tr><tr><td>Q<sub>1</sub> - Q<sub>3</sub></td><td>Dimensionless fluxes of C<sub>1</sub>, C<sub>2</sub>, and C<sub>3</sub> on Γ<sub>1</sub>, Γ<sub>2</sub>, and Γ<sub>3</sub></td><td>—</td></tr><tr><td>W<sub>1</sub> - W<sub>3</sub></td><td>Approximated dimensionless transport coefficients</td><td>—</td></tr></tbody></table>

# 2 Formulation

This section describes the mathematical formulation for the optimization used in this study. First, Sect. 2.1 describes the battery configuration and design variables to be optimized. Section 2.2 then presents the analytical model of the ASSBs and clarifies the governing equations to be solved. Finally, Sect. 2.3 formulates an optimization problem aimed at minimizing the internal losses of the battery.

## 2.1 Design variables settings

Figure 1 shows the battery configuration that is the analysis domain of this research. The battery consists of three regions: composite anode (negative electrode) Ω<sub>n</sub>, electrolyte Ω<sub>e</sub>, and composite cathode (positive electrode) Ω<sub>p</sub>. The composite cathode is a composite region of active and electrolyte materials, while the composite anode is the composite region of active, electrolyte, and conductive additive materials. In the electrolyte, only the electrolyte material exists. The three materials composing the battery have different roles. The active material stores Li that moves between the cathode and anode during charging and discharging, and the electrolyte carries Li ions.
<figure><img src="image_1.png" /><figcaption>Fig. 1: Model of All-solid-state-batteries</figcaption></figure>

The conductive additive material is a highly conductive material added to increase the electronic conductivity of the composite anode.

The design goal of this study is to find the optimal distributions of the three materials in the composite anode $\Omega_n$. That is, the volume fraction of each material is assumed to be not constant in the composite anode, but depends on the position coordinates, and then the volume fractions of these materials are set to the design variables. Let $\epsilon_{a,n}$, $\epsilon_{e,n}$, and $\epsilon_c$ be the volume fractions of the active, electrolyte, and conductive additive materials in the composite anode, respectively. Here, this paper assumes that there are no voids in the composite anode.

## 2.2 Governing equations

<figure><img src="image_2.png" /><figcaption>Fig. 2: Domains and boundaries definitions</figcaption></figure>

The analysis in this study is conducted in a two-dimensional domain, as shown in Fig. 2. As mentioned in the previous section, the computational domain $\Omega$ consists of three regions: composite anode $\Omega_n$, composite cathode $\Omega_p$, and electrolyte $\Omega_e$, with parameters $L_n$, $L_p$, and $L_e$, respectively. The height is $L_h$. The left boundary of the composite anode is represented as $\Gamma_l$, and the right boundary of the composite cathode is represented as $\Gamma_r$.
as $\Gamma_r$. The interfaces between the composite anode and the electrolyte and between the composite cathode and the electrolyte are denoted by $\Gamma_{n-e}$ and $\Gamma_{p-e}$, respectively.

Next, time and spatial scales are defined within the computational domain, as shown in the following equations:

$$t \in \mathcal{T}[0, t_{\max}] \subseteq R_{\ge 0}, \qquad (1)$$

$$x = (x_1, x_2) \in \Omega \subseteq R^2. \qquad (2)$$

where $t_{\max}$ is the terminated time of the discharge analysis. The termination condition for the analysis of the discharge process to determine $t_{\max}$ is described in detail in the latter part of this section.

There are five variables required in this analysis model. Four of those are the Li concentrations, $c_n$ and $c_p$, and electrical potentials, $\phi_n$ and $\phi_p$, defined in each composite region. The fifth variable is the electrical potential $\phi_e$ in the electrolyte phase, defined in all regions where the electrolyte is present. In this study, as in Ishida et al. (2024), the transport number of Li ions is assumed to be unity (Kanno, 2023), and Li and electron diffusion follows Fick's law (Fuller and Harb, 2018). The governing equations, initial conditions and boundary conditions to be solved are described as follows:

$$\epsilon_{a,n} \frac{\partial c_n}{\partial t} = \nabla \cdot (D_n^{\text{eff}} \nabla c_n) - \frac{j^{\text{Li}}}{F} \quad \text{in } \mathcal{T} \times \Omega_n, \qquad (3)$$

$$\epsilon_{a,p} \frac{\partial c_p}{\partial t} = \nabla \cdot (D_p^{\text{eff}} \nabla c_p) - \frac{j^{\text{Li}}}{F} \quad \text{in } \mathcal{T} \times \Omega_p, \qquad (4)$$

$$\nabla \cdot (\sigma_n^{\text{eff}} \nabla \phi_n) = j^{\text{Li}} \quad \text{in } \mathcal{T} \times \Omega_n, \qquad (5)$$

$$\nabla \cdot (\sigma_p^{\text{eff}} \nabla \phi_p) = j^{\text{Li}} \quad \text{in } \mathcal{T} \times \Omega_p, \qquad (6)$$

$$\nabla \cdot (\kappa^{\text{eff}} \nabla \phi_e) = -j^{\text{Li}} \quad \text{in } \mathcal{T} \times \Omega, \qquad (7)$$

$$c_n = c_{n,\max} \times \text{SOC}_{n,\max} \quad \text{in } \mathcal{T}(0) \times \Omega_n, \qquad (8)$$

$$c_p = c_{p,\max} \times \text{SOC}_{p,\min} \quad \text{in } \mathcal{T}(0) \times \Omega_p, \qquad (9)$$

$$D_n^{\text{eff}} \nabla c_n \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{l,o,n-e}, \qquad (10)$$

$$D_p^{\text{eff}} \nabla c_p \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{r,o,p-e}, \qquad (11)$$

$$\phi_n = 0 \quad \text{on } \mathcal{T} \times \Gamma_l, \qquad (12)$$

$$\sigma_n^{\text{eff}} \nabla \phi_n \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{o,n-e}, \qquad (13)$$

$$\sigma_p^{\text{eff}} \nabla \phi_p \cdot \mathbf{n} = -\frac{I}{A_g} \quad \text{on } \mathcal{T} \times \Gamma_r, \qquad (14)$$

$$\sigma_p^{\text{eff}} \nabla \phi_p \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{o,p-e}, \qquad (15)$$

$$\kappa^{\text{eff}} \nabla \phi_e \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{l,r,out}, \qquad (16)$$

where $\epsilon_{a,p}$ is the volume fraction of active material in the composite cathode and $\mathbf{n}$ is the outward normal vector. $c_{n,\max}$ and $c_{p,\max}$ are the upper limits of the Li concentration that can exist in the composite anode and cathode material, respectively.
SOC$_{n,max}$ and SOC$_{p,min}$ are the upper and lower limits of the state of charge in the composite anode and cathode, respectively. A$_{g}$ is the area of the current collector. j$^{Li}$ and η are the reaction current density and overvoltage, respectively, and are defined as follows:

$$j^{\text{Li}} = a_s^{\text{eff}} i_0 \left[ \exp\left(\frac{0.5F}{RT}\eta\right) - \exp\left(-\frac{0.5F}{RT}\eta\right) \right], \quad (17)$$

$$\eta = \begin{cases} \phi_n - \phi_e - U_n(c_n) & \text{in } \Omega_n \\ 0 & \text{in } \Omega_e \\ \phi_p - \phi_e - U_p(c_p) & \text{in } \Omega_p, \end{cases} \quad (18)$$

where i$_{0}$ is the exchange current density. Then, F, R, and T are the Faraday constant, gas constant, and temperature, respectively. In this study, the temperature is assumed to be constant. U$_{n}$ and U$_{p}$ are functions that determine the standard equilibrium electrode potential, which depends on the type of material used in each composite region and the Li concentration.

As shown in Eq. 14, the condition given to the system is assumed to be a constant current. The magnitude of the given current is related to the capacity of the battery, as follows:

$$I = \frac{Q_{\text{cell}}}{3600[\text{s}]} C_{\text{rate}}, \quad (19)$$

where C$_{rate}$ is a dimensionless parameter that determines the magnitude of the battery discharge. By the definition, C$_{rate}$ of 1 represents the magnitude of the current at which the battery is fully discharged in one hour. This study defines the capacity of the battery Q$_{cell}$ based on the amount of active material in the composite cathode, as shown in the following equation:

$$Q_{\text{cell}} = c_{\text{p,max}}(\text{SOC}_{\text{p,max}} - \text{SOC}_{\text{p,min}})\varepsilon_{\text{a,p}}V_{\text{p}}F, \quad (20)$$

where V$_{p}$ is the volume of the composite cathode region. SOC$_{p,max}$ are the upper limits of the state of charge in the composite cathode. Then, the initial conditions for the Li concentrations of the composite anode and cathode are given as Eqs. 8 and 9, respectively. Namely, the initial conditions for this study are set to be fully charged in order to calculate the discharge process.

Next, we explain the termination conditions for the discharge analysis. As mentioned above, the magnitude of the current is given so that the battery should discharge in certain time. However, due to the non-linearity of the governing equations, the time taken for the discharge is not an explicit value, and can change during the optimization process. Therefore, in this study, the discharge calculation is terminated when the terminal voltage, which decreases during the discharge process, falls below a specific threshold value V$_{terminate}$, as shown in the following equation:

$$\frac{\int_{\Gamma_r} \phi_p \, d\Gamma}{\int_{\Gamma_r} d\Gamma} \le V_{\text{terminate}}. \quad (21)$$
Finally, we explain the material interpolation method. That is, the following for-
mulations reveal which variables depend on which material and how they depend on
the volume fraction of the material. In the analytical model in this paper, there are
two types of effective coefficients that require material interpolation, i.e., the transport
coefficient and the surface area of active material per unit volume in the composite
anode, where the design variables are defined, and the composite cathode. For the
transport coefficient, there are five effective coefficients in the governing equations. The
diffusion coefficients of Li in the composite anode and composite cathode depend on
the volume fraction of active material present in each domain. The electronic conduc-
tivity in the composite anode depends on the volume fractions of both active material
and conductive additive material. The electronic conductivity in the composite cath-
ode depends on the volume fraction of the active material. The Li-ion conductivity in
the electrolyte depends on the volume fraction of the electrolyte in each region. Then,
the surface area of active material per unit volume also depends on the volume frac-
tion of active material within each domain and the particle size of the active material.
For simplicity, this study assumes the particle size is constant in the composite region.
These relationships can be expressed in the following equations:

$$
D_n^{\text{eff}} = D_n H_1(\varepsilon_{a,n}(\mathbf{x})), \qquad (22)
$$

$$
D_p^{\text{eff}} = D_p H_1(\varepsilon_{a,p}), \qquad (23)
$$

$$
\sigma_n^{\text{eff}} = \sigma_n H_2(\varepsilon_{a,n}(\mathbf{x}), \varepsilon_c(\mathbf{x})), \quad (24)
$$

$$
\sigma_p^{\text{eff}} = \sigma_p H_1(\varepsilon_{a,p}), \qquad (25)
$$

$$
\kappa^{\text{eff}} = \begin{cases} \kappa H_1(\varepsilon_{e,n}(\mathbf{x})) & \text{in } \Omega_n \\ \kappa H_1(\varepsilon_{e,e}) & \text{in } \Omega_e \\ \kappa H_1(\varepsilon_{e,p}) & \text{in } \Omega_p, \end{cases} \qquad (26)
$$

$$
a_s^{\text{eff}} = \begin{cases} \frac{A(\varepsilon_{a,n}(\mathbf{x}))}{L} & \text{in } \Omega_n \\ \frac{A(\varepsilon_{a,p})}{L} & \text{in } \Omega_p, \end{cases} \qquad (27)
$$

where $H_1$, $H_2$ and $A$ are the dimensionless interpolation functions for the transport coefficients and surface area proposed in this paper, respectively. The details of these functions are described in Sect. 3. $L$ is the length of the unit cell defined in Sect. 3.1. $D_n$ and $D_p$ are the Li diffusion coefficient of active material in the composite anode and cathode, respectively. $\sigma_n$ and $\sigma_p$ are the conductivity of active material in the composite anode and cathode, respectively. Moreover, $\kappa$ is the ionic conductivity of electrolyte material. To clarify our paper description, it should be remarked that $\varepsilon_{a,n}$, $\varepsilon_c$, and $\varepsilon_{e,n}$ depend on the coordinates of the position, while $\varepsilon_{a,p}$, $\varepsilon_{e,p}$, and $\varepsilon_{e,e}$ are constants. This is because the design variable is set to the volume fraction distribution only in the composite anodes in this paper.

2.3 Optimization problem

In this section, we formulate an optimization problem aiming at minimizing the energy
loss during the discharge process. For simplicity of notation, the following variables
are defined:

$$ \tilde{c}_p = \frac{\int_{\Omega_p} c_p \, d\Omega}{\int_{\Omega_p} c_{p,\max} \, d\Omega}, \qquad (28) $$

$$ \tilde{\phi}_s = \frac{\int_{\Gamma_r} \phi_s \, d\Gamma}{\int_{\Gamma_r} d\Gamma}, \qquad (29) $$

where $\tilde{c}_p$ is the average Li concentration in the composite anode, and $\tilde{\phi}_s$ is the average potential of the active material at the rightmost boundary. Using these variables, the optimization problem is formulated as follows:

$$ \inf_{\epsilon_{a,n}, \epsilon_c} J = \frac{\int_{\mathcal{T}} (U_p(\tilde{c}_p) - \tilde{\phi}_s) \, dt}{t_{\max}}, \qquad (30) $$

$$ \text{s.t.} \quad \underline{\epsilon}_{a,n} \le \epsilon_{a,n} \le \overline{\epsilon}_{a,n}, \qquad (31) $$

$$ \underline{\epsilon}_c \le \epsilon_c \le \overline{\epsilon}_c, \qquad (32) $$

$$ \epsilon_{a,n}(\boldsymbol{x}) + \epsilon_c(\boldsymbol{x}) + \epsilon_{e,n}(\boldsymbol{x}) = 1. \qquad (33) $$

The objective function is set as the overvoltage during the discharge process. In other words, the objective function is defined as the difference of the standard electrode potential based on the average Li concentration in the cathode and the terminal voltage. The overvoltage corresponds to the amount of energy required to drive the battery reaction, and a small absolute value of this value means that the energy loss inside the battery is small (Alizadeh et al., 2023). Since this study considers only the discharge process, the overvoltage is negative in the composite cathode. Therefore, the absolute value is removed beforehand by multiplying by -1. It should be noted that the sign of the objective function changes when the charging process is considered. The design variables are $\epsilon_{a,n}$ and $\epsilon_c$. As mentioned in Sect. 2.1, the objective of this study is to obtain the distribution of the volume fractions of active material, electrolyte, and conductive additive in the composite anode. Since the presence of voids does not positively affect the optimized solution, we assume that there are no voids and that Eq. 33 is satisfied. Therefore, the number of the design variables can be reduced. Here, $\overline{\epsilon}_{a,n}$ and $\overline{\epsilon}_c$ are the upper limits of the volume fraction of active material and conductive additive, respectively, and $\epsilon_{a,n}$ and $\epsilon_c$ are the lower limits of those, respectively.

Next, we derive the sensitivity of the objective function formulated above using the Lagrange and continuous adjoint methods. First, the Lagrangian is defined as follows:

$$ \mathcal{L} = J + \int_{\mathcal{T}} \int_{\Omega_n} a_1 \left( \epsilon_{a,n} \frac{\partial (c_n)}{\partial t} - \nabla \cdot (D_n^{\text{eff}} \nabla c_n) + \frac{j^{\text{Li}}}{F} \right) d\Omega dt \\ + \int_{\mathcal{T}} \int_{\Omega_p} a_2 \left( \epsilon_{a,p} \frac{\partial (c_p)}{\partial t} - \nabla \cdot (D_p^{\text{eff}} \nabla c_p) + \frac{j^{\text{Li}}}{F} \right) d\Omega dt \\ + \int_{\mathcal{T}} \int_{\Omega_n} a_3 \left( \nabla \cdot (\sigma_n^{\text{eff}} \nabla \phi_n) - j^{\text{Li}} \right) d\Omega dt $$
$$ 
\begin{aligned}
&+ \int_{\mathcal{T}} \int_{\Omega_{\mathrm{p}}} a_4 (\nabla \cdot (\sigma_{\mathrm{p}}^{\mathrm{eff}} \nabla \phi_{\mathrm{p}}) - j^{\mathrm{Li}}) \mathrm{d}\Omega \mathrm{d}t \\
&+ \int_{\mathcal{T}} \int_{\Omega} a_5 (\nabla \cdot (\kappa_{\mathrm{m}}^{\mathrm{eff}} \nabla \phi_{\mathrm{e,m}}) + j^{\mathrm{Li}}) \mathrm{d}\Omega \mathrm{d}t,
\end{aligned} 
\quad (34) $$

where $a_1, a_2, a_3, a_4$, and $a_5$ are the adjoint variables. This Lagrangian can be differentiated by the design variables and transformed using the chain rule as follows:

$$ 
\begin{aligned}
\frac{d\mathcal{L}}{d\varepsilon_{[\cdot]}} = & \frac{\partial\mathcal{L}}{\partial\varepsilon_{[\cdot]}} + \frac{\partial\mathcal{L}}{\partial c_{\mathrm{n}}} \frac{\partial c_{\mathrm{n}}}{\partial\varepsilon_{[\cdot]}} + \frac{\partial\mathcal{L}}{\partial c_{\mathrm{p}}} \frac{\partial c_{\mathrm{p}}}{\partial\varepsilon_{[\cdot]}} \\
& + \frac{\partial\mathcal{L}}{\partial\phi_{\mathrm{n}}} \frac{\partial\phi_{\mathrm{n}}}{\partial\varepsilon_{[\cdot]}} + \frac{\partial\mathcal{L}}{\partial\phi_{\mathrm{p}}} \frac{\partial\phi_{\mathrm{p}}}{\partial\varepsilon_{[\cdot]}} + \frac{\partial\mathcal{L}}{\partial\phi_{\mathrm{e}}} \frac{\partial\phi_{\mathrm{e}}}{\partial\varepsilon_{[\cdot]}},
\end{aligned} 
\quad (35) $$

where $\varepsilon_{[\cdot]}$ represents $\varepsilon_{\mathrm{a,n}}$ and $\varepsilon_{\mathrm{c}}$. Since it is difficult to derive $\frac{\partial c_{\mathrm{n}}}{\partial\varepsilon_{[\cdot]}}, \frac{\partial c_{\mathrm{p}}}{\partial\varepsilon_{[\cdot]}}, \frac{\partial\phi_{\mathrm{n}}}{\partial\varepsilon_{[\cdot]}}, \frac{\partial\phi_{\mathrm{p}}}{\partial\varepsilon_{[\cdot]}}$, and $\frac{\partial\phi_{\mathrm{e}}}{\partial\varepsilon_{[\cdot]}}$, we aim to find the distribution of adjoint variables that satisfy the conditions $\frac{\partial\mathcal{L}}{\partial c_{\mathrm{n}}} = 0, \frac{\partial\mathcal{L}}{\partial c_{\mathrm{p}}} = 0, \frac{\partial\mathcal{L}}{\partial\phi_{\mathrm{n}}} = 0, \frac{\partial\mathcal{L}}{\partial\phi_{\mathrm{p}}} = 0$, and $\frac{\partial\mathcal{L}}{\partial\phi_{\mathrm{e}}} = 0$ as the following equations show:

$$ - \varepsilon_{\mathrm{a,n}} \frac{\partial a_1}{\partial t} = \nabla \cdot (D_{\mathrm{n}}^{\mathrm{eff}} \nabla a_1) + \tilde{\eta} M \frac{\partial U_{\mathrm{n}}}{\partial c_{\mathrm{n}}} \quad \text{in } \mathcal{T} \times \Omega_{\mathrm{n}}, \quad (36) $$

$$ - \varepsilon_{\mathrm{a,p}} \frac{\partial a_2}{\partial t} = \nabla \cdot (D_{\mathrm{p}}^{\mathrm{eff}} \nabla a_2) + \tilde{\eta} M \frac{\partial U_{\mathrm{p}}}{\partial c_{\mathrm{p}}} - \frac{1}{\partial \tilde{c}_{\mathrm{p}}} L_{\mathrm{h}} c_{\mathrm{p,max}} t_{\mathrm{max}} \quad \text{in } \mathcal{T} \times \Omega_{\mathrm{p}}, \quad (37) $$

$$ \nabla \cdot (\sigma_{\mathrm{n}}^{\mathrm{eff}} \nabla a_3) = -\tilde{\eta} M \quad \text{in } \mathcal{T} \times \Omega_{\mathrm{n}}, \quad (38) $$

$$ \nabla \cdot (\sigma_{\mathrm{p}}^{\mathrm{eff}} \nabla a_4) = -\tilde{\eta} M \quad \text{in } \mathcal{T} \times \Omega_{\mathrm{p}}, \quad (39) $$

$$ \nabla \cdot (\kappa^{\mathrm{eff}} \nabla a_5) = \tilde{\eta} M \quad \text{in } \mathcal{T} \times \Omega, \quad (40) $$

$$ a_1 = 0 \quad \text{in } \mathcal{T}(t_{\mathrm{max}}) \times \Omega_{\mathrm{n}}, \quad (41) $$

$$ a_2 = 0 \quad \text{in } \mathcal{T}(t_{\mathrm{max}}) \times \Omega_{\mathrm{p}}, \quad (42) $$

$$ D_{\mathrm{n}}^{\mathrm{eff}} \nabla a_1 \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{\mathrm{l,o,n-e}}, \quad (43) $$

$$ D_{\mathrm{p}}^{\mathrm{eff}} \nabla a_2 \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{\mathrm{r,o,p-e}}, \quad (44) $$

$$ a_3 = 0 \quad \text{on } \mathcal{T} \times \Gamma_{\mathrm{l}}, \quad (45) $$

$$ \sigma_{\mathrm{n}}^{\mathrm{eff}} \nabla a_3 \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{\mathrm{o,n-e}}, \quad (46) $$

$$ \sigma_{\mathrm{p}}^{\mathrm{eff}} \nabla a_4 \cdot \mathbf{n} = - \frac{1}{L_{\mathrm{h}} t_{\mathrm{max}}} \quad \text{on } \mathcal{T} \times \Gamma_{\mathrm{r}}, \quad (47) $$

$$ \sigma_{\mathrm{p}}^{\mathrm{eff}} \nabla a_4 \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{\mathrm{o,p-e}}, \quad (48) $$

$$ \kappa^{\mathrm{eff}} \nabla a_5 \cdot \mathbf{n} = 0 \quad \text{on } \mathcal{T} \times \Gamma_{\mathrm{l,r,o}}, \quad (49) $$

$$ M = \frac{\partial j^{\mathrm{Li}}}{\partial \eta}, \quad (50) $$
$$ \tilde{\eta} = \begin{cases} \frac{a_1}{F} - a_3 + a_5 & \text{in } \Omega_n \\ 0 & \text{in } \Omega_e \\ \frac{a_2}{F} - a_4 + a_5 & \text{in } \Omega_p. \end{cases} \qquad (51) $$

Finally, the sensitivities can be derived as:

$$ \begin{aligned} \frac{d\mathcal{L}}{d\varepsilon_{a,n}} = \frac{\partial\mathcal{L}}{\partial\varepsilon_{a,n}} &= -\int_{\Omega_n} a_1 c_n|_{t=0} d\Omega + \int_{\mathcal{T}} \int_{\Omega_n} \left( -c_n \frac{\partial a_1}{\partial t} \right. \\ & \qquad \left. + \frac{\partial D_n^{\text{eff}}}{\partial \varepsilon_{a,n}} \nabla a_1 \cdot \nabla c_n - \frac{\partial \sigma_n^{\text{eff}}}{\partial \varepsilon_{a,n}} \nabla a_3 \cdot \nabla \phi_n \right. \\ & \qquad \left. - \frac{\partial \kappa^{\text{eff}}}{\partial \varepsilon_{a,n}} \nabla a_5 \cdot \nabla \phi_e + \left( \frac{a_1}{F} - a_3 + a_5 \right) \frac{\partial j^{\text{Li}}}{\partial \varepsilon_{a,n}} \right) d\Omega dt \\ &:= \int_{\Omega_n} J_1' d\Omega, \end{aligned} \qquad (52) $$

$$ \begin{aligned} \frac{d\mathcal{L}}{d\varepsilon_c} = \frac{\partial\mathcal{L}}{\partial\varepsilon_c} &= \int_{\mathcal{T}} \int_{\Omega_n} \left( -\frac{\partial\sigma_n^{\text{eff}}}{\partial\varepsilon_c} \nabla a_3 \cdot \nabla\phi_n - \frac{\partial\kappa^{\text{eff}}}{\partial\varepsilon_c} \nabla a_5 \cdot \nabla\phi_e \right) d\Omega dt \\ &:= \int_{\Omega_n} J_2' d\Omega, \end{aligned} \qquad (53) $$

The details of the sensitivity derivation are noted in the Appendix of Ishida et al. (2024).

## 3 Material interpolation scheme

This section explains the derivation scheme of functions $H_1$ and $H_2$, related to the transport coefficients that depend on the volume fraction, and function $A$, which is associated with surface area, both derived from microstructure analysis. Since the problem addressed in this study deals with the distribution of the material volume fraction, it is necessary to utilize interpolation functions that accurately reflect actual physical phenomena when interpolating effective properties at intermediate densities. The TO method for ASSBs proposed by Ishida et al. (2024) was based on the homogenization method, which assumes that the unit cells are arranged periodically to calculate macroscopic properties. However, as Kimura (2023) illustrated, the actual configuration of ASSBs has a nonuniform and not periodical material distribution, making it difficult to utilize the homogenization method.

In this study, spherical materials are randomly distributed in a predefined three-dimensional domain to represent the nonuniform material distribution, and the transport coefficient value in three directions and the surface area for the material distribution are evaluated. This procedure is repeated a certain number of times to obtain datasets, and the interpolation functions are obtained by applying the least-squares method to these datasets. Moreover, the presence of Ersatz material is considered to
ensure that the derived interpolation functions do not adversely affect the optimization process. It should be noted that this procedure is conducted before the optimization process performs as Sect. 4 denotes.

In Sect. 3.1, the definition of the computational domain and the representation of the material in the microstructure are explained. Next, Sect. 3.2 describes how the transport coefficients in the three directions and the magnitude of the surface area are evaluated. Then, Sect. 3.3 describes the algorithm for obtaining the interpolation function and the consideration of Ersatz material (Allaire et al., 1997).

## 3.1 Random material distribution schemes

<figure><img src="image_3.png" /><figcaption>Fig. 3: Models of the microstructure analysis</figcaption></figure>

Figure 3 shows the computational domain for the microstructure analysis. In this study, the variables are defined as dimensionless ones to generalize this method. In other words, the size of the cell is set to a dimensionless unit cell with the size of 1 based on the original cell size L [m]. The spatial coordinates are defined by $y_1, y_2$, and $y_3$, as follows:

$$ y = (y_1, y_2, y_3) \in \mathbf{Y} \subseteq \mathbb{R}^3. \qquad (54) $$

When finding $H_1(\varepsilon)$ and $A(\varepsilon)$, which depend on the volume fraction of single type material, a material distribution with a constant diameter is randomly generated in this space, and when finding $H_2(\varepsilon_a, \varepsilon_c)$, which depends on the volume fraction of both active and conductive additive materials, materials with two different diameters and material properties are randomly arranged.

First, the random material generation scheme for single material type is explained. To represent the presence or absence of material, a normalized density function $\rho(y) \in [0, 1]$ is defined, which takes the value 1 in regions where the material is present and 0 in void regions. The following describes how to obtain a random $\rho$ distribution.

Based on an arbitrary random function, we obtain a set of sphere center coordinates $[\tilde{y}_{1,1}, \tilde{y}_{2,1}, \tilde{y}_{3,1}, \dots, \tilde{y}_{1,g}, \tilde{y}_{2,g}, \tilde{y}_{3,g}]$ with values between 0 and 1. Here, g is the number
of spheres to be generated, and the constant diameter of the spheres is $R_s$. Based on this dataset, a density function $\bar{\rho}$ that takes a value of 1 inside the sphere is defined as follows:

$$ \bar{\rho} = \begin{cases} 1 & (\text{if } (y_1 - \tilde{y}_{1,m})^2 + (y_2 - \tilde{y}_{2,m})^2 + (y_3 - \tilde{y}_{3,m})^2 \le \left(\frac{R_s}{2}\right)^2) \\ 0 & (\text{others}). \end{cases} \quad (55) $$

where m is (1, 2, ..., g). Here, $\bar{\rho}$ does not have a smooth isosurface in the discretized space for numerical analysis. Therefore, $\bar{\rho}$ is smoothed using a Helmholtz filter, as shown in the following equations:

$$ -R_s^2 \nabla_Y^2 \hat{\rho} + \hat{\rho} = \bar{\rho} \quad \text{in } \mathbf{Y}, \quad (56) $$

$$ \nabla_Y \hat{\rho} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_{1,2,3,4,5,6}, \quad (57) $$

where $\nabla_Y$ is the gradient in space $\mathbf{Y}$. Next, the smoothed $\hat{\rho}$ is binarized using the following Heaviside function to represent the presence or absence of material:

$$ \rho = \frac{1}{2} \tanh \left\{ \beta_Y \left( \hat{\rho} - \frac{1}{2} \right) \right\} + \frac{1}{2}, \quad (58) $$

where $\beta_Y$ is the steepness of the Heaviside function. Figure 4 illustrates the several material distributions generated using the abovementioned method.

<figure><img src="image_4.png" /><figcaption>Fig. 4: Randomly generated material distributions for single type of material in a unit cell with $R_s$ of 0.2</figcaption></figure>

Next, we explain the random material generation scheme for two types of materials within the unit cell for deriving the function $H_2(\varepsilon_a, \varepsilon_c)$. To represent the presence or absence of active and conductive additive materials, a normalized density function $\rho(\mathbf{y}) \in [0, 1]$ is defined, which takes the value 1 in regions where the conductive additive material is present, 0.5 in regions where the active material is present, and 0 in void regions. Following the abovementioned approach, a normalized density function is generated.
Based on an arbitrary random function, we obtain sets of sphere center coordinates [tilde{y}_{11,1}, tilde{y}_{12,1}, tilde{y}_{13,1}, ..., tilde{y}_{11,g_1}, tilde{y}_{12,g_1}, tilde{y}_{13,g_1}] for the conductive additive material and [tilde{y}_{21,1}, tilde{y}_{22,1}, tilde{y}_{23,1}, ..., tilde{y}_{21,g_2}, tilde{y}_{22,g_2}, tilde{y}_{23,g_2}] for the active material with values between 0 and 1. Here, g_1 and g_2 are the numbers of spheres to be generated for conductive additive and active materials, respectively, and the constant diameters of the spheres for conductive additive material and active material are R_{s1} and R_{s2}, respectively. Based on these datasets, a density function p that takes a value of 1 inside the sphere of conductive additive material and 0.5 inside the sphere of active material is defined as follows:

$$ \bar{\rho} = \begin{cases} 1 & (\text{if } (y_1 - \tilde{y}_{1,m})^2 + (y_2 - \tilde{y}_{2,m})^2 + (y_3 - \tilde{y}_{3,m})^2 > \left(\frac{R_{s2}}{2}\right)^2 \\ \cap (y_1 - \tilde{y}_{1,n})^2 + (y_2 - \tilde{y}_{2,n})^2 + (y_3 - \tilde{y}_{3,n})^2 \le \left(\frac{R_{s1}}{2}\right)^2) \\ 0.5 & (\text{if } (y_1 - \tilde{y}_{1,m})^2 + (y_2 - \tilde{y}_{2,m})^2 + (y_3 - \tilde{y}_{3,m})^2 \le \left(\frac{R_{s2}}{2}\right)^2) \\ 0 & (\text{others}), \end{cases} \quad (59) $$

where n is (1, 2, ..., g_1) and m is (1, 2, ..., g_2), respectively. As mentioned before, p does not have a smooth isosurface in the discretized space for numerical analysis. Therefore, p is smoothed using a Helmholtz filter as defined in Eqs. 56 and 57. After that, p is binarized using the following continuous piece-wise Heaviside projection:

$$ \rho = \frac{1}{2} \left\{ \frac{1}{2} \tanh \left( \beta_Y \left( 2\hat{\rho} - \frac{1}{2} \right) \right) + \frac{1}{2} \right\} + \frac{1}{2} \left\{ \frac{1}{2} \tanh \left( \beta_Y \left( 2\hat{\rho} - \frac{3}{2} \right) \right) + \frac{1}{2} \right\}. \quad (60) $$

Figure 5 illustrates several material distributions with two types of materials generated using the aforementioned method.

<figure><img src="image_5.png" /><figcaption>Fig. 5: Randomly generated material distributions in a unit cell for two material types: R_{s1} = 0.1 (red) and R_{s2} = 0.2 (gray)</figcaption></figure>
## 3.2 Estimation of normalized transport coefficient and surface area

This section describes a method for obtaining the transport coefficients in three directions and the surface area of materials for the given material distribution value $\rho$ obtained in the previous section. In order to obtain the transport coefficients, this paper solves the dimensionless diffusion equations. In other words, the homogeneous Dirichlet boundary condition is given for the direction in which the transport coefficient is to be obtained, and the flux in that direction is calculated. The surfaces in the other directions are applied with the zero Neumann boundary condition. To evaluate the transport coefficient in the $y_1$-direction, we solve the following diffusion equations:

$$ q_1 = -K \nabla_Y C_1, \qquad (61) $$

$$ \nabla_Y \cdot q_1 = 0 \qquad \text{in } Y, \qquad (62) $$

$$ C_1 = 0 \qquad \text{on } \Gamma_1, \qquad (63) $$

$$ C_1 = 1 \qquad \text{on } \Gamma_2, \qquad (64) $$

$$ q_1 \cdot n = 0 \qquad \text{on } \Gamma_{3,4,5,6}, \qquad (65) $$

where $C_1$ is the dimensionless scalar field and $K$ is the dimensionless transport coefficient. As for the material distribution of single type of material, $K$ is defined as 1 where there is material and a smaller value $K_{min}$ where there is no material as follows:

$$ K = \rho + (1 - \rho) K_{\min}. \qquad (66) $$

Whereas, as for the material distribution of both active material and conductive additive material, $K$ is defined as $K_{max}$ where there is conductive additive material, 1 where active material exists, and a smaller value $K_{min}$ where there is no material as follows:

$$ K = (1 - K_{\min}) \left\{ \frac{1}{2} \tanh \left( \beta_Y \left( 2\rho - \frac{1}{2} \right) \right) + \frac{1}{2} \right\} + K_{\min} \\ + (K_{\max} - 1) \left\{ \frac{1}{2} \tanh \left( \beta_Y \left( 2\rho - \frac{3}{2} \right) \right) + \frac{1}{2} \right\}, \qquad (67) $$

where $K_{max}$ is the electronic conductivity of the conductive additive material normalized by that of the active material as follows:

$$ K_{\max} = \frac{\sigma_c}{\sigma_a}, \qquad (68) $$

where $\sigma_c$ is the electronic conductivity of the conductive additive material. The flux of $C_1$, $Q_1$, in the $y_1$ direction is then calculated based on the following equation:

$$ W_1 = Q_1 = \int_{\Gamma_1} q_1 \cdot n \, d\Gamma, \qquad (69) $$
where the dimensionless effective transport coefficient $W_1$ in the $y_1$ direction coincides with $Q_1$ since the size of the unit cell is 1 and the difference between the values of the Dirichlet conditions at boundaries $\Gamma_1$ and $\Gamma_2$ is set to 1.

In the same way, the dimensionless effective transport coefficients, $W_2$ and $W_3$, in the $y_2$ and $y_3$ directions, respectively, are obtained. The following dimensionless diffusion equations for the dimensionless scalar fields, $C_2$ and $C_3$, are solved:

$$ q_2 = -K \nabla_y C_2, \qquad (70) $$

$$ \nabla_y \cdot q_2 = 0 \qquad \text{in } Y, \qquad (71) $$

$$ C_2 = 0 \qquad \text{on } \Gamma_3, \qquad (72) $$

$$ C_2 = 1 \qquad \text{on } \Gamma_4, \qquad (73) $$

$$ q_2 \cdot n = 0 \qquad \text{on } \Gamma_{1,2,5,6}, \qquad (74) $$

$$ q_3 = -K \nabla_y C_3, \qquad (75) $$

$$ \nabla_y \cdot q_3 = 0 \qquad \text{in } Y, \qquad (76) $$

$$ C_3 = 0 \qquad \text{on } \Gamma_5, \qquad (77) $$

$$ C_3 = 1 \qquad \text{on } \Gamma_6, \qquad (78) $$

$$ q_3 \cdot n = 0 \qquad \text{on } \Gamma_{1,2,3,4}. \qquad (79) $$

Then, the flux of $C_2$ and $C_3$ in the $y_2$ and $y_3$ directions, $Q_2$ and $Q_3$, respectively, are evaluated by the following equations, respectively:

$$ W_2 = Q_2 = \int_{\Gamma_3} q_2 \cdot n \, d\Gamma, \qquad (80) $$

$$ W_3 = Q_3 = \int_{\Gamma_5} q_3 \cdot n \, d\Gamma. \qquad (81) $$

Similarly, $Q_2$ and $W_2$, and $Q_3$ and $W_3$ coincide. Based on the transport coefficients obtained from these analyses, isotropic transport coefficient $\tilde{H}$ is derived using the $L_2$ norm as follows:

$$ \tilde{H}(\rho(y; \varepsilon)) = \frac{\sqrt{W_1^2 + W_2^2 + W_3^2}}{\sqrt{3}}, \qquad (82) $$

Next, the surface area of randomly generated material distribution in the unit cell is evaluated using the following equation, which is one of the bases of the Hamilton-Jacobi equation (Osher and Sethian, 1988):

$$ \tilde{A}(\rho(y; \varepsilon)) = \int_{\Omega} |\nabla \rho| \, d\Omega. \qquad (83) $$

### 3.3 Algorithm for obtaining interpolation functions

This section shows the algorithm details for obtaining the interpolation functions, $H_1(\varepsilon)$, $H_2(\varepsilon_a, \varepsilon_c)$ and $A(\varepsilon)$, which depend on the volume fraction of the material. First,
<figure><img src="image_6.png" /><figcaption>Fig. 6: Algorithm for making the material interpolation function $H_1$ and $A$</figcaption></figure>

we explain the algorithm for finding $H_1(\varepsilon)$ and $A(\varepsilon)$ which depend on only one volume fraction. Figure 6 shows the flowchart of the proposed algorithm for obtaining $H_1(\varepsilon)$ and $A(\varepsilon)$. Initially, the integer $i$ is defined as 1 to count the number of loops. Then, a target value for the volume fraction, $\hat{\varepsilon}_i$, which is used when the material distribution is created, is set. Next, the spheres' center coordinates are defined randomly so that the material volume approaches the target value. After that, the actual volume of the generated material distribution $\varepsilon_i$ is obtained. Note that there is no guarantee that $\hat{\varepsilon}_i$, and $\varepsilon_i$ coincide because of the existence of overlaps. Next, filtering and binarization with the Heaviside function are performed on the generated material distribution as shown in Eqs. 56–58, Eqs. 61–66, and Eqs. 70–79 are then solved to obtain the effective transport coefficient $\tilde{H}_{1i}$ and surface area $\tilde{A}_i$ from Eqs. 82 and 83, respectively. If the number of loops $i$ is less than $M$, the target value of the volume fraction $\hat{\varepsilon}_i$ is updated, and the same procedure is conducted.

At the end of $M$ loops, a dataset of volume fractions $\varepsilon_i(i = 1, ..., M)$, transport coefficients $\tilde{H}_{1i}(i = 1, ..., M)$, and surface areas $\tilde{A}_i(i = 1, ..., M)$ for $M$ varieties of material distributions $\rho_i(i = 1, ..., M)$ can be obtained. Here, the functions $\tilde{H}_1$ and $\tilde{A}$ are interpolated by a $U$-dimensional polynomial of the volume fractions, as shown in
the following equations:

$$ \tilde{H}_1 = \sum_{m=0}^{U} h_m \varepsilon^m, \qquad (84) $$

$$ \tilde{A} = \sum_{m=0}^{U} w_m \varepsilon^m, \qquad (85) $$

where $h_m$ and $w_m$ are the coefficients of the respective polynomial equations. Then, the least-squares method is applied to each equation using the acquired dataset to obtain $\tilde{H}_1$ and $\tilde{A}$.

Next, the algorithm for finding $H_2(\varepsilon_a, \varepsilon_c)$ which depends on the volume fractions of both active material and conductive additive material. Initially, the integer $i$ is defined as 1 to count the number of loops. Then, target values for the volume fraction of active material and conductive additive material, $\hat{\varepsilon}_{a,i}$ and $\hat{\varepsilon}_{c,i}$, are set. Next, the spheres' center coordinates are defined randomly so that the material volume approaches the target values. After that, the actual volumes of the generated material distribution, $\varepsilon_{a,i}$ and $\varepsilon_{c,i}$, are obtained. Next, filtering and binarization with the continuous piece-wise Heaviside function are performed on the generated material distribution as shown in Eqs. 56–57, Eq. 60, Eqs. 61–65, Eq. 67, and Eqs. 70–79 are then solved to obtain the effective transport coefficient $\tilde{H}_{2i}$ from Eq. 82. If the number of loops $i$ is less than $M$, the target values of the volume fractions, $\hat{\varepsilon}_{a,i}$ and $\hat{\varepsilon}_{c,i}$, are updated, and the same procedure is conducted.

At the end of $M$ loops, a dataset of volume fractions, $\varepsilon_{a,i}$ and $\varepsilon_{c,i}$ ($i = 1, ..., M$), and transport coefficients $\tilde{H}_{2i}$ ($i = 1, ..., M$) for $M$ varieties of material distributions $\rho_i$ ($i = 1, ..., M$) can be obtained. Here, the function $\tilde{H}_2$ is interpolated by a $U_a$-dimensional polynomial for $\varepsilon_a$ and a $U_c$-dimensional polynomial for $\varepsilon_c$, as shown in the following equation:

$$ \tilde{H}_2(\varepsilon_a, \varepsilon_c) = \sum_{m=0}^{U_a} \sum_{j=0}^{U_c} h_{a,m} h_{c,j} \varepsilon_a^m \varepsilon_c^j, \qquad (86) $$

where $h_{a,m}$ and $h_{c,j}$ are the coefficients of the polynomial equations. At the end, the least squares method is applied to the equation using the acquired dataset to obtain $\tilde{H}_2$.

However, the functions $\tilde{H}_1$ and $\tilde{H}_2$ obtained by this method have a fatal drawback in optimization problems. That is, it does not guarantee to have a monotonic increase property and positive value. The interpolation function obtained from pure least squares may take negative values and may oscillate slightly. These characteristics are different from the physical interpretation and have a harmful impact on the search for an optimization solution that is updated based on the gradient of variables. Therefore, we consider the monotonic increase and positive constraints when the least squares for $\tilde{H}_1$ and $\tilde{H}_2$ are conducted. Moreover, the presence of Ersatz material should be considered to stabilize the convergence of the numerical calculation during
the optimization process as follows:

$$H_1 = \tilde{H}_1(\varepsilon) + H_{1,\min}, \qquad (87)$$

$$H_2 = \tilde{H}_2(\varepsilon_a, \varepsilon_c) + H_{2,\min}, \qquad (88)$$

where $H_{1,\min}$ and $H_{2,\min}$ are the given small values. Similarly, since $\tilde{A}$ has a minimal value around zero, $A_{\min}$ is also added to stabilize the convergence of the numerical calculation as follows:

$$A = \tilde{A}(\varepsilon) + A_{\min}. \qquad (89)$$

Note that any constraint is not considered when $\tilde{A}(\varepsilon)$ is derived using least square.

## 4 Numerical implementations

<figure><img src="image_7.png" /><figcaption>Fig. 7: Algorithm for the optimization</figcaption></figure>

This section describes the numerical implementation details for the optimization problem formulated in Sect. 2.3. Figure 7 shows the optimization algorithm used in this paper. First, the material interpolation function values, $H_1$, $H_2$ and $A$, are computed using the method explained in Sect. 3. Then, initial configurations of the design variables, $\rho_{a,n}$ and $\rho_c$, are given. Next, the governing equations Eqs. 3–18 are solved based on the design variables from $t = 0$ until the termination condition Eq. 21 is satisfied. After that, the objective function $J$ is evaluated. Then, the termination
criterion is checked using the following conditions:

$$ \left\| \frac{J_{r+1} - J_r}{J_r} \right\| \le \omega, \qquad (90) $$

where the index r is the number of the optimization procedure, and ω is the given convergence criterion. If Eq. 90 is satisfied, the convergence is deemed to be achieved and the optimization process is terminated. On the other hand, if the equation is not satisfied, the adjoint equations Eqs. 36–51 are solved. Then, the sensitivities, $J'_1$ and $J'_2$, are calculated based on the state and adjoint variables. Finally, the design variables are updated using the Method of Moving Asymptotes (Svanberg, 1987), and the governing equations are solved again.

# 5 Numerical examples

In this section, several numerical examples are presented to confirm the availability of the proposed method. Section 5.1 presents the parameter settings used in this numerical example. Section 5.2 presents detailed numerical examples and discusses the physical meanings of the obtained material interpolation functions, $H_1$, $H_2m$ and $A$, based on the method described in Sect. 3. Finally, Sect. 5.3 shows the optimal solution obtained by conducting the proposed optimization method, i.e., the optimal distribution of the volume fractions of each three materials existing in the composite anode. In this paper, for the microstructural analysis, FreeFEM++ is used, while the optimization was performed using COMSOL Multiphysics Simulations with Matlab®, a solver with a functional standing for battery analysis.

## 5.1 Numerical settings

Table 2 shows the values of the constants used in the numerical examples. This paper used lithium cobalt oxide for the cathode and graphite carbon for the anode. The equilibrium potential equations, $U_n$ and $U_p$, were based on the material data provided by COMSOL Multiphysics (Multiphysics, 1998). It was assumed that electrochemical reactions are taking place inside the battery based on the following:

$$ \mathrm{LiCoO_2} + 6 \mathrm{C} \rightleftharpoons \mathrm{Li}_{1-x}\mathrm{CoO_2} + \mathrm{Li_xC_6}. \qquad (91) $$

Here, the right-to-left reaction is the discharge reaction, in which Li is deintercalated from the anode region. In contrast, the left-to-right reaction is the charge reaction, in which Li is intercalated into the anode region.

## 5.2 Material interpolation function profile

Fig. 8 shows the discretization mesh used for the microstructural analysis. The tetrahedral elements with size of 1/80 were utilized. The dimensionless numbers $R_s$, $R_{s1}$, and $R_{s2}$ was set to 0.2, 0.1, and 0.2, respectively, where $R_s$ and $R_{s2}$ are the ratio of the actual material diameter of the electrolyte and active materials to the size of the
<table><caption>Table 2: Parameter settings</caption><thead><tr><th>Parameter</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>L<sub>n</sub></td><td>85.0</td><td>µm</td></tr><tr><td>L<sub>e</sub></td><td>45.0</td><td>µm</td></tr><tr><td>L<sub>p</sub></td><td>65.0</td><td>µm</td></tr><tr><td>L<sub>h</sub></td><td>100</td><td>µm</td></tr><tr><td>L</td><td>10</td><td>µm</td></tr><tr><td>A<sub>g</sub></td><td>1</td><td>m</td></tr><tr><td>ε<sub>a,p</sub></td><td>0.44</td><td>–</td></tr><tr><td>ε<sub>e,e</sub></td><td>0.60</td><td>–</td></tr><tr><td>ε<sub>e,p</sub></td><td>0.40</td><td>–</td></tr><tr><td>ε<sub>a,n</sub></td><td>0.55</td><td>–</td></tr><tr><td>ε<sub>a,n</sub></td><td>0</td><td>–</td></tr><tr><td>ε<sub>c</sub></td><td>0.40</td><td>–</td></tr><tr><td>ε<sub>c</sub></td><td>0</td><td>–</td></tr><tr><td>C<sub>rate</sub></td><td>1</td><td>–</td></tr><tr><td>SOC<sub>n,max</sub></td><td>0.98</td><td>–</td></tr><tr><td>SOC<sub>n,min</sub></td><td>0</td><td>–</td></tr><tr><td>SOC<sub>p,max</sub></td><td>1</td><td>–</td></tr><tr><td>SOC<sub>p,min</sub></td><td>0.43</td><td>–</td></tr><tr><td>c<sub>n,max</sub></td><td>31507</td><td>mol/m<sup>3</sup></td></tr><tr><td>c<sub>p,max</sub></td><td>56250</td><td>mol/m<sup>3</sup></td></tr><tr><td>F</td><td>96485.3</td><td>C/mol</td></tr><tr><td>R</td><td>8.3143</td><td>J/(mol·K)</td></tr><tr><td>T</td><td>298</td><td>K</td></tr><tr><td>V<sub>terminate</sub></td><td>2.7</td><td>V</td></tr><tr><td>ω</td><td>1.0 × 10<sup>-3</sup></td><td>–</td></tr><tr><td>i<sub>0,n</sub></td><td>36</td><td>A/m<sup>2</sup></td></tr><tr><td>i<sub>0,p</sub></td><td>26</td><td>A/m<sup>2</sup></td></tr><tr><td>D<sub>n</sub></td><td>2.09 × 10<sup>-15</sup></td><td>m<sup>2</sup>/s</td></tr><tr><td>D<sub>p</sub></td><td>5.00 × 10<sup>-13</sup></td><td>m<sup>2</sup>/s</td></tr><tr><td>σ<sub>n</sub></td><td>1.00 × 10<sup>-1</sup></td><td>S/m</td></tr><tr><td>σ<sub>p</sub></td><td>1.13 × 10<sup>-2</sup></td><td>mS/cm</td></tr><tr><td>σ<sub>c</sub></td><td>100</td><td>S/m</td></tr><tr><td>κ</td><td>2.00 × 10<sup>-2</sup></td><td>S/m</td></tr><tr><td>M</td><td>(80, 400)</td><td>–</td></tr><tr><td>β<sub>Y</sub></td><td>30</td><td>–</td></tr><tr><td>H<sub>1,min</sub></td><td>0.001</td><td>–</td></tr><tr><td>H<sub>2,min</sub></td><td>0.001</td><td>–</td></tr><tr><td>A<sub>min</sub></td><td>1</td><td>–</td></tr><tr><td>K<sub>min</sub></td><td>0.001</td><td>–</td></tr><tr><td>K<sub>max</sub></td><td>1000</td><td>–</td></tr></tbody></table>

unit cell L to be considered, and $R_{s1}$ is the ratio of the actual material diameter of the conductive additive material to the size of the unit cell. Here, as an example in this study, we obtained the characteristics of the physical properties observed in a unit cell with a size L of 10 [µm], assuming that the diameters of the electrolyte and active materials were 2 [µm] and that of the conductive additive was 1 [µm]. The number of microstructure analysis loops M was set to 80 for $H_1$ and A, and 400 for $H_2$. The ratio of the electronic conductivity of the conductive additive material to that of the active material, $K_{max}$, was set to 1000.
<figure><img src="image_8.png" /><figcaption>Fig. 8: Mesh for the unit cell</figcaption></figure>

Figure 9 presents the obtained data points and the constructed interpolation functions, $\tilde{H}_1$, $\tilde{H}_2$, and $\tilde{A}$, as a function of the volume fraction. As Fig. 9 shows, monotonically increasing functions $\tilde{H}_1$ and $\tilde{H}_2$, and an upward convex function $\tilde{A}$ interpolate the acquired datasets well. It should be noted that an additional data-point of $\tilde{H}_1 = 1$ at $\varepsilon = 1$ was considered when creating the interpolation function $\tilde{H}_1$, since the transport coefficient at $\varepsilon = 1$ is the same as that of the bulk material.

Next, we consider the effect of the orders of the interpolation function shown in Eqs. 84–86 on the function configurations. The order used in the least square is determined by trial and error approach, and there is no existing method to determine it uniquely. Therefore, in this numerical example, Akaike information criterion (AIC) was utilized as a criterion.

<table><caption>**Table 3:** The coefficients $h$ for $\tilde{H}_1$ obtained by least squares for different values $U$ and the values of AIC</caption><thead><tr><th>U</th><th>h<sub>0</sub></th><th>h<sub>1</sub></th><th>h<sub>2</sub></th><th>h<sub>3</sub></th><th>h<sub>4</sub></th><th>h<sub>5</sub></th><th>h<sub>6</sub></th><th>h<sub>7</sub></th><th>AIC</th></tr></thead><tbody><tr><td>4</td><td>1.9 × 10<sup>-9</sup></td><td>5.7 × 10<sup>-2</sup></td><td>-0.74</td><td>3.4</td><td>-1.6</td><td>-</td><td>-</td><td>-</td><td>-711</td></tr><tr><td>5</td><td>1.5 × 10<sup>-10</sup></td><td>0.12</td><td>-1.3</td><td>5.2</td><td>-4.0</td><td>0.98</td><td>-</td><td>-</td><td>-713</td></tr><tr><td>6</td><td>4.2 × 10<sup>-9</sup></td><td>9.2 × 10<sup>-8</sup></td><td>0.27</td><td>-3.4</td><td>17</td><td>-21</td><td>8.3</td><td>-</td><td>-728</td></tr><tr><td>7</td><td>3.7 × 10<sup>-10</sup></td><td>2.6 × 10<sup>-2</sup></td><td>-0.24</td><td>0.11</td><td>5.2</td><td>-2.1</td><td>-6.5</td><td>4.6</td><td>-726</td></tr></tbody></table>

Tables 3 and 4 present the coefficients and AIC values for $\tilde{H}_1$ and $\tilde{A}$ obtained through least squares fitting for various orders $U$. For $\tilde{H}_1$, the sixth-order function yielded the smallest AIC value, while for $\tilde{A}$, the smallest AIC value was observed at $U = 3$. Accordingly, the orders of $\tilde{H}_1$ and $\tilde{A}$ were set to 6 and 3, respectively.

Table 5 presents the AIC values for $\tilde{H}_2$ obtained through least squares fitting with various orders of $U_a$ and $U_c$. The combination of $U_a = 2$ and $U_c = 6$ yielded the smallest AIC value. Accordingly, this study adopts this pair for the least squares fitting
<figure><img src="image_9.png" /><figcaption>(a) For the dimensionless transport coefficient which depends on single type of material, $U = 6$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) For the dimensionless surface area, $U = 3$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) For the dimensionless transport coefficient which depends on two types of materials, $U_a = 2$ and $U_c = 6$</figcaption></figure>

<caption>Fig. 9: Comparison of data-points from microstructure analysis and the approximated material interpolation functions, $\tilde{H}_1$, $\tilde{H}_2$ and $\tilde{A}$</caption>

of $\tilde{H}_2$, resulting in the following equation:

$$ \tilde{H}_2(\varepsilon_a, \varepsilon_c) = 0.15 - 23\varepsilon_c^1 + 1100\varepsilon_c^2 - 14000\varepsilon_c^3 + 70000\varepsilon_c^4 - 130000\varepsilon_c^5 + 87000\varepsilon_c^6 \\ + \varepsilon_a(-0.051 - 2.1\varepsilon_c^1 + 420\varepsilon_c^2 - 8000\varepsilon_c^3 + 54000\varepsilon_c^4 - 140000\varepsilon_c^5 + 130000\varepsilon_c^6) \\ + \varepsilon_a^2(0.18 + 64\varepsilon_c^1 - 2700\varepsilon_c^2 + 35000\varepsilon_c^3 - 180000\varepsilon_c^4 + 390000\varepsilon_c^5 - 310000\varepsilon_c^6). \tag{92} $$
<table><caption>**Table 4:** The coefficients *w* for Ā obtained by least squares for different values *U* and the values of AIC</caption><thead><tr><th><i>U</i></th><th><i>w</i><sub>0</sub></th><th><i>w</i><sub>1</sub></th><th><i>w</i><sub>2</sub></th><th><i>w</i><sub>3</sub></th><th><i>w</i><sub>4</sub></th><th><i>w</i><sub>5</sub></th><th>AIC</th></tr></thead><tbody><tr><td>2</td><td>-0.57</td><td>40</td><td>-35</td><td>-</td><td>-</td><td>-</td><td>-199</td></tr><tr><td>3</td><td>0.015</td><td>32</td><td>-13</td><td>-16</td><td>-</td><td>-</td><td>-260</td></tr><tr><td>4</td><td>-0.018</td><td>32</td><td>-17</td><td>-9.2</td><td>-39</td><td>-</td><td>-258</td></tr><tr><td>5</td><td>0.022</td><td>31</td><td>-4.7</td><td>-46</td><td>43</td><td>-21</td><td>-257</td></tr></tbody></table>

<table><caption>**Table 5:** AIC values for varying orders of least squares, $U_a$ and $U_c$</caption><thead><tr><th><i>U</i><sub>a</sub></th><th><i>U</i><sub>c</sub></th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><td>2</td><td>1151</td><td>1140</td><td>1132</td><td>1134</td></tr><tr><td>3</td><td>1161</td><td>1151</td><td>1145</td><td>1150</td></tr></tbody></table>

Figure 10 illustrates the plots of the proposed function $H_1$ and the SIMP method with penalty parameters $p = 1.5$, 2, and 3 as a function of the volume fraction. The plots indicate that the behavior of the proposed function $H_1$ differs from that of the SIMP method. For small volume fractions, $H_1$ closely resembles the distribution with a penalty parameter of $p = 3$, whereas for larger volume fractions, it approaches the configuration with $p = 2$. In the context of this study, where intermediate material properties are incorporated into the optimal solution, these differences significantly influence the optimized configuration, as discussed in Section 5.3.2. Consequently, the material interpolation function derived through the physical interpretation presented here plays a critical role in optimization problems for battery designs.

## 5.3 Optimization results

In this section, the optimized configurations of the material volume fraction obtained by the proposed optimization method are shown to confirm the effectiveness and characteristics of the proposed method. Section 5.3.1 exhibits the performance of the optimized configuration by comparing the performance with a uniform distribution. Section 5.3.2 shows the importance of proposing the material interpolation functions. Section 5.3.3 demonstrates the generality of the proposed method by presenting the variations in the optimal solution when the conductivities of the conductive additive material and active material are altered.

### 5.3.1 Performance evaluation of optimized material distribution: Comparison with uniform distribution

Figure 11 shows the initial configurations of the three materials, and Fig. 12 shows the optimized configurations obtained by the proposed method. The average volume fractions of the electrolyte, electrolytic additive, and active material for the obtained
<figure><img src="image_12.png" /><figcaption>**Fig. 10:** Plots of the proposed material interpolation function $H_1$ for the dimensionless transport coefficient and SIMP method with $p = 1.5, 2$, and $3$ as a function of volume fraction</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11:** The initial material distributions for the optimization</figcaption></figure>

results after optimization were 0.595, 0.085, and 0.320, respectively. The obtained results exhibit that the electrolyte is distributed mainly on the side close to the electrolyte phase, while the active and conductive additive materials are distributed near the current collector. This tendency is an attempt to reduce the resistance of electron and Li ion conduction during the discharge process and is consistent with the results described in Clausnitzer et al. (2023). The active material is distributed in a wavy pattern. This distribution attempts to supply Li evenly in the composite anode by arranging the active material, which is responsible for holding Li, in a wave-like form. Such a material distribution design was realized through the high degree of freedom offered by the TO method.

Next, the performance of the obtained configuration is demonstrated by comparing it to a commonly used homogeneously mixed battery, where this configuration is named 'uniform material distribution' in this paper. To ensure consistent battery conditions, including the theoretical capacity, the volume of each material in the uniform material distribution is set to be equal to the obtained optimal solution. The uniform
<figure><img src="image_14.png" /><figcaption>Fig. 12: The obtained material distributions by the proposed method</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13: The uniform material configurations for the comparison</figcaption></figure>

material distributions are shown in Fig. 13. The objective functions of the optimized structure and the comparison example were 0.661 and 0.738, respectively, indicating that the proposed method reduces overvoltage by 10%. Furthermore, the SOC of the anode at the termination time for the optimized structure and the uniform material distribution were 0.0448 and 0.138, respectively. This result demonstrates that the proposed TO method enhances the capacity by 11% compared to the uniform distribution, while maintaining the same theoretical capacity.

Figure 14 shows a detailed history of each objective function term during the discharge process. The horizontal axis indicates the SOC of the anode, and the perpendicular axis indicates the voltage. The diamond symbols in the graph indicate the history of $U_p(\tilde{c}_p)$ and $\tilde{\phi}_s$ for the optimized configuration, and the circle plots the history of $U_p(\tilde{c}_p)$ and $\tilde{\phi}_s$ for the uniform material distribution. As shown in Eq. 30, the objective function is the difference between $U_p(\tilde{c}_p)$ and $\tilde{\phi}_s$ integrated over time during the discharge process. As the plot shows, the optimized configuration could reduce the difference between the two evaluation variables compared to the uniform distribution. This result implies that the proposed optimization method provides a design guideline for a suitable material distribution.
<figure><img src="image_16.png" /><figcaption>Fig. 14: Histories of $U_p(\tilde{c}_p)$ and discharge curve $\tilde{\phi}_s$ during the discharge process for the uniform and optimized material distributions as functions of SOC$_n$</figcaption></figure>

### 5.3.2 Comparison and analysis of optimal solutions obtained using different material interpolation functions

This paper formulated the interpolation function using microstructure analysis instead of using the SIMP method as the material interpolation function. To demonstrate the importance of the proposed interpolation functions, here we show the comparison between the optimized configuration obtained using the proposed interpolation functions and one obtained using the SIMP method. The material interpolation functions, $H_1$ and $H_2$, can be expressed using the SIMP method as follows:

$$H_1(\varepsilon) = (1 - H_{1,\min})\varepsilon^p + H_{1,\min}, \qquad (93)$$

$$H_2(\varepsilon_a, \varepsilon_p) = (1 - H_{2,\min})\varepsilon_a^p + K_{\max}(1 - H_{2,\min})\varepsilon_c^p + H_{2,\min}. \qquad (94)$$

Figure 15 illustrates the optimized configurations obtained using the proposed functions and the SIMP method where the penalty parameter $p$ is set to either 3 or 1.5. As can be seen in Fig. 15, the optimization configuration significantly depends on the interpolation functions used in the optimization calculation. These changes are due to the inclusion of intermediate volume fraction materials in the optimal structure. In other words, in optimization, where design variables are updated based on the gradient of variables, interpolation functions have a significant influence on the direction of the solution search.

Next, the optimization method was confirmed to investigate the validity of our implementation using cross-checking for the obtained objective function value. For each optimized configuration obtained by the proposed interpolated functions, and SIMP method $p$ of 3 and 1.5, the objective function is evaluated by computing the discharge process with the respective interpolation function. The horizontal axis in Table 6 denotes the interpolation function values used for optimization, and the vertical axis shows the interpolation function values used when the recalculation was performed.
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** The obtained material distributions by the optimization method using the proposed function $H_1$ and $H_2$, and the penalty parameter $p$ of 3 and 1.5 with SIMP method</figcaption></figure>

Next, the optimization method was validated by cross-checking the obtained objective function values to confirm the accuracy of our implementation. For each optimized configuration derived from the proposed interpolated functions and the SIMP method with $p = 3$ and $p = 1.5$, the objective function was re-evaluated by simulating the discharge process using the respective interpolation functions. In Table 6, the columns represent the interpolation functions used for optimization, while the rows indicate those applied for recalculation. When the objective function values of each line are compared, the objective function value evaluated using the interpolation functions applied to the optimization is smaller than the objective function values re-evaluated using other methods. This suggests that the implemented optimization algorithm is
capable of deriving the optimal solution for each formulation. These results demonstrate that the proposed implementation scheme is valid and has been executed correctly.

Thus, as shown in Fig. 15, the optimized configuration highly depends on the interpolation function profile, and it is important to use interpolation functions, $H_1$, $H_2$ and $A$, matching the physical interpretation when applying the proposed TO method to the design of actual batteries.

<table><caption>**Table 6:** The re-calculated objective function values using the proposed material interpolation function $H_1$ and $H_2$, and $p$ of 3 and 1.5 with the SIMP method for the obtained material distributions</caption><thead><tr><th>calculation</th><th>optimization</th><th>proposed functions</th><th>$p=3$</th><th>$p=1.5$</th></tr></thead><tbody><tr><td>proposed functions</td><td></td><td>0.661</td><td>0.683</td><td>0.697</td></tr><tr><td>$p=3$</td><td></td><td>0.695</td><td>0.694</td><td>0.728</td></tr><tr><td>$p=1.5$</td><td></td><td>0.624</td><td>0.645</td><td>0.613</td></tr></tbody></table>

### 5.3.3 Relation between electronic conductivity values of materials and optimized structures

As many researchers devote much effort to materials development, battery performance strongly depends on the physical properties of the materials (Chabot et al., 2013). That is, the physical properties of the materials used in batteries are constantly evolving. For the proposed optimization method to be used in battery development, it is important to show the possibility of deriving the optimized structure with various material properties. Therefore, in this section, the optimized configurations are investigated, setting the various electronic conductivity values of the active material and conductive additive material as (0.001 [S/m], 1 [S/m]), (0.01 [S/m], 10 [S/m]), (0.1 [S/m], 100 [S/m]), and (1 [S/m], 1000 [S/m]). Figure 16 shows the obtained optimal material distributions for each pair of electronic conductivities. Comparing Figs. 16b, 16e, 16h and 16k, it can be seen that the volume of the conductive additive in the optimized material distribution decreases as the electronic conductivity increases. This tendency is obtained because the resistance due to electron transfer decreases as the electronic conductivity increases, and the volume of the electrolyte and active material can be increased to compensate. Thus, the proposed method can derive the optimal structure for different material properties.

## 6 Conclusions

This paper proposed a TO method that provides the optimized material distribution in the composite anode of ASSBs incorporating microstructure analysis. We achieved the following remarks:
<figure><img src="image_18.png" /><figcaption>Fig. 16: The obtained material distributions for the various electronic conductivity values of the active material and conductive additive material: $(\sigma_n, \sigma_c)=(0.001 \text{ [S/m]}, 1 \text{ [S/m]})$, $(0.01 \text{ [S/m]}, 10 \text{ [S/m]})$, $(0.1 \text{ [S/m]}, 100 \text{ [S/m]})$, and $(1 \text{ [S/m]}, 1000 \text{ [S/m]})$</figcaption></figure>
* We formulated an optimization problem to minimize the overvoltage, corresponding to the energy loss in the battery charge-discharge reaction. The design variables were set to the normalized densities representing the distributions of the volume fractions of three materials, i.e., electrolyte, conductivity additive, and active materials in the composite anode, assuming the absence of voids. The design sensitivities used for updating the design variables were derived using the Lagrange and continuous adjoint methods.

* Based on microstructural analysis, material interpolation functions were formulated to describe the relationships between the volume fractions of the materials and their transport coefficients, as well as the relationships between volume fractions and the surface area of the material distributions. To develop these functions, 80 material distributions with a single material type and 400 distributions with two material types were randomly generated. The interpolation functions were then derived using the least-squares method under constraints of monotonicity and positivity. Additionally, the Ersatz material concept was introduced to ensure numerical stability during the optimization process. The validity for determining the order used in the least squares was confirmed using AIC values.

* The optimization algorithm was developed based on the formulation for the optimization problem. The optimized configurations of material distributions in the composite anode represent an unprecedented result that could only be obtained using the TO method. A comparison of the derived optimized distribution with a uniform distribution under the same conditions indicated that the proposed method could reduce the overvoltage by 10% and enhance the capacity by 11%.

* The optimized results derived using the SIMP method were compared with those obtained by the proposed method to demonstrate the necessity of appropriate material interpolation functions. The obtained distributions were utterly varied, even though the same physical properties were used in the optimization. Namely, it was found that the material interpolation function, which represents the intermediate material's physical properties, significantly affects the optimization when the intermediate material is assumed to be included in the solution.

* The proposed method's effectiveness in battery development, where new materials are developed by various ways and the physical properties of materials evolve daily, was confirmed by changing the electronic conductivity values of the conductive additive material and active material, and deriving the optimized material distributions. As the electronic conductivity increased, the conductive additive's volume in the optimized distributions decreased. This result was obtained because the resistance to electron transfer decreases as conductivity increases, and the volume is replaced by materials that play another role, such as active and electrolyte material. Thus, the proposed method could derive the optimized configurations using a variety of different material properties.

**Acknowledgments.** This work was supported by JSPS KAKENHI Grant Number JP22H01384.
## Declarations

* Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.

* Replication of results
The presented results were obtained using COMSOL Multiphysics Simulations with Matlab® and we can provide them on reasonable request.

## References

Alizadeh, M., Charoen-amornkitt, P., Suzuki, T., et al. (2023) A numerical simulation of evolution processes and entropy generation for optimal architecture of an electrochemical reaction-diffusion system: Comparison of two optimization strategies. Journal of The Electrochemical Society 170(11):114520. https://doi.org/10.1149/1945-7111/ad0a7c, URL https://dx.doi.org/10.1149/1945-7111/ad0a7c

Allaire, G., Bonnetier, E., Francfort, G., et al. (1997) Shape optimization by the homogenization method. Numerische Mathematik 76:27–68

Bendsøe, M. P. (1989) Optimal shape design as a material distribution problem. Structural Optimization 1:193–202. https://doi.org/10.1007/BF01650949

Bendsøe, M. P., Sigmund, O. (1999) Material interpolation schemes in topology optimization. Archive of Applied Mechanics 69:635–654. https://doi.org/10.1007/s004190050248

Bendsøe, M. P., Kikuchi, N. (1988) Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering 71(2):197–224. https://doi.org/10.1016/0045-7825(88)90086-2

Borrvall, T., Petersson, J. (2003) Topology optimization of fluids in stokes flow. International Journal for Numerical Methods in Fluids 41(1):77–107

Cano, Z. P., Banham, D., Ye, S., et al. (2018) Batteries and fuel cells for emerging electric vehicle markets. Nature Energy 3(4):279–289. https://doi.org/10.1038/s41560-018-0108-1

Chabot, V., Farhad, S., Chen, Z., et al. (2013) Effect of electrode physical and chemical properties on lithium-ion battery performance. International Journal of Energy Research 37(14):1723–1736

Chung, D.-W., Ebner, M., Ely, D. R., et al. (2013) Validity of the bruggeman relation for porous electrodes. Modelling and Simulation in Materials Science and Engineering 21(7):074009. https://doi.org/10.1088/0965-0393/21/7/074009, URL
Clausnitzer, M., Mücke, R., Al-Jaljouli, F., et al. (2023) Optimizing the composite cathode microstructure in all-solid-state batteries by structure-resolved simulations. Batteries & Supercaps 6(11):e202300167.
https://doi.org/https://doi.org/10.1002/batt.202300167, URL https://chemistry-europe.onlinelibrary.wiley.com/doi/abs/10.1002/batt.202300167, https://chemistry-europe.onlinelibrary.wiley.com/doi/pdf/10.1002/batt.202300167

Deng, J., Ren, X., Lin, H., et al. (2024) Functionally gradient materials for sustainable and high-energy rechargeable lithium batteries: Design principles, progress, and perspectives. Journal of Energy Chemistry 99:426-449. https://doi.org/https://doi.org/10.1016/j.jechem.2024.07.062, URL https://www.sciencedirect.com/science/article/pii/S2095495624005436

Deysher, G., Ridley, P., Ham, S.-Y., et al. (2022) Transport and mechanical aspects of all-solid-state lithium batteries. Materials Today Physics

Fuller, T. F., Harb, J. N. (2018) Electrochemical Engineering. John Wiley & Sons

Giraldo, S., Nakagawa, K., Vásquez, F. A., et al. (2021) Preparation of composite electrodes for all-solid-state batteries based on sulfide electrolytes: An electrochemical point of view. Batteries 7(4). https://doi.org/10.3390/batteries7040077, URL https://www.mdpi.com/2313-0105/7/4/77

Guo, C., Liu, H., Guo, Q., et al. (2022) Investigations on a novel cold plate achieved by topology optimization for lithium-ion batteries. Energy 261:125097. https://doi.org/https://doi.org/10.1016/j.energy.2022.125097, URL https://www.sciencedirect.com/science/article/pii/S0360544222019922

Ishida, N., Furuta, K., Kishimoto, M., et al. (2024) Topology optimization for all-solid-state-batteries using homogenization method. Structural and Multidisciplinary Optimization 67(9):162. https://doi.org/https://doi.org/10.1007/s00158-024-03864-1

Iwai, H., Kuroyanagi, A., Saito, M., et al. (2011) Power generation enhancement of solid oxide fuel cell by cathode/electrolyte interface modification in mesoscale assisted by level set-based optimization calculation. Journal of Power Sources 196(7):3485-3495. https://doi.org/https://doi.org/10.1016/j.jpowsour.2010.12.024

Janek, J., Zeier, W. G. (2016) A solid future for battery development. Nature Energy 1(9):1-4. https://doi.org/https://doi.org/10.1038/nenergy.2016.141

Kanno, R. (2023) Developments in new materials for electrochemistry and energy storage devices. Electrochemistry 91(10):102001-102001. https://doi.org/10.5796/electrochemistry.23-00064
Kimura, Y. (2023) Development of measurement/analysis techniques to bridge micro-
scopic and macroscopic electrochemical phenomena in energy storage devices.
Electrochemistry 91(10):101005–101005. https://doi.org/10.5796/electrochemistry.
23-00074

Lin, T. Y., Baker, S. E., Duoss, E. B., et al. (2022) Topology optimization of 3d flow fields for flow batteries. Journal of The Electrochemical Society 169(5):050540. https://doi.org/10.1149/1945-7111/ac716d, URL https://dx.doi.org/10.1149/1945-7111/ac716d

Liu, P., Sherman, E., Jacobsen, A. (2009) Design and fabrication of multifunctional structural batteries. Journal of Power Sources 189(1):646–650. https://doi.org/ https://doi.org/10.1016/j.jpowsour.2008.09.082, URL https://www.sciencedirect.com/science/article/pii/S0378775308018405, selected Papers presented at the 14th INTERNATIONAL MEETING ON LITHIUM BATTERIES (IMLB-2008)

Lu, L., Yamamoto, T., Otomori, M., et al. (2013) Topology optimization of an acoustic metamaterial with negative bulk modulus using local resonance. Finite Elements in Analysis and Design 72:1–12. https://doi.org/https://doi.org/10.1016/j.finel.2013.04.005, URL https://www.sciencedirect.com/science/article/pii/S0168874X13000607

Lv, C., Zhou, X., Zhong, L., et al. (2022) Machine learning: an advanced platform for materials development and state prediction in lithium-ion batteries. Advanced Materials 34(25):2101474. https://doi.org/https://doi.org/10.1002/adma.202101474

Matsumori, T., Kondoh, T., Kawamoto, A., et al. (2013) Topology optimization for fluid-thermal interaction problems under constant input power. Structural and Multidisciplinary Optimization 47:571–581. https://doi.org/https://doi.org/10.1007/s00158-013-0887-8

Minnmann, P., Quillman, L., Burkhardt, S., et al. (2021) Editors’ choice—quantifying the impact of charge transport bottlenecks in composite cathodes of all-solid-state batteries. Journal of The Electrochemical Society 168(4):040537

Mo, X., Zhi, H., Xiao, Y., et al. (2021) Topology optimization of cooling plates for battery thermal management. International Journal of Heat and Mass Transfer 178:121612. https://doi.org/https://doi.org/10.1016/j.ijheatmasstransfer.2021.121612, URL https://www.sciencedirect.com/science/article/pii/S0017931021007158

Multiphysics, C. (1998) Introduction to comsol multiphysics®. COMSOL Multi-physics, Burlington, MA, accessed Feb 9:2018

Nanda, S., Gupta, A., Manthiram, A. (2021) Anode-free full cells: a pathway to high-energy density lithium-metal batteries. Advanced Energy Materials 11(2):2000804. https://doi.org/https://doi.org/10.1002/aenm.202000804
Ong, S. P. (2019) Accelerating materials science with high-throughput computations and machine learning. Computational Materials Science 161:143-150. https://doi.org/10.1016/j.commatsci.2019.01.013, URL https://www.sciencedirect.com/science/article/pii/S0927025619300138

Osher, S., Sethian, J. A. (1988) Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations. Journal of Computational Physics 79(1):12-49. https://doi.org/10.1016/0021-9991(88)90002-2, URL https://www.sciencedirect.com/science/article/pii/0021999188900022

Pejman, R., Najafi, A. R. (2023) Multiphysics topology optimization of a multifunctional structural battery composite. Structural and Multidisciplinary Optimization 66(3):46. https://doi.org/10.1007/s00158-023-03490-3

Reisacher, E., Kaya, P., Knoblauch, V. (2023) Percolation behavior of a sulfide electrolyte-carbon additive matrix for composite cathodes in all-solid-state batteries. Batteries 9(12). https://doi.org/10.3390/batteries9120595, URL https://www.mdpi.com/2313-0105/9/12/595

Ren, Y., Danner, T., Moy, A., et al. (2023) Oxide-based solid-state batteries: a perspective on composite cathode architecture. Advanced Energy Materials 13(1):2201939. https://doi.org/10.1002/aenm.202201939

Roy, T., Troya, M. A., Worsley, M. A., et al. (2022) Topology optimization for the design of porous electrodes. Structural and Multidisciplinary Optimization 65(6):171

Shahzad, K., Mardare, A. I., Hassel, A. W. (2024) Accelerating materials discovery: combinatorial synthesis, high-throughput characterization, and computational advances. Science and Technology of Advanced Materials: Methods 4(1):2292486. https://doi.org/10.1080/27660400.2023.2292486

Shi, T., Tu, Q., Tian, Y., et al. (2020) High active material loading in all-solid-state battery electrode via particle size optimization. Advanced Energy Materials 10(1):1902881. https://doi.org/10.1002/aenm.201902881

Snyder, J., Baechle, D., Wetzel, E., et al. (2008) Multifunctional structural composite batteries for us army applications. In: Proceedings of the 26th Army Science Conference, pp 1-8

Snyder, J. F., Carter, R. H., Wetzel, E. D. (2007) Electrochemical and mechanical behavior in mechanically robust solid polymer electrolytes for use in multifunctional structural batteries. Chemistry of Materials 19(15):3793-3801. https://doi.org/10.1021/cm070213o

Svanberg, K. (1987) The method of moving asymptotes—a new method for structural optimization. International Journal for Numerical Methods in Engineering
24(2):359-373. https://doi.org/10.1002/nme.1620240207

Tjaden, B., Cooper, S. J., Brett, D. J., et al. (2016) On the origin and application of the bruggeman correlation for analysing transport phenomena in electrochemical systems. Current Opinion in Chemical Engineering 12:44-51. https://doi.org/10.1016/j.coche.2016.02.006, URL https://www.sciencedirect.com/science/article/pii/S2211339816300119, nanotechnology / Separation Engineering

Wadbro, E., Berggren, M. (2006) Topology optimization of an acoustic horn. Computer Methods in Applied Mechanics and Engineering 196(1-3):420-436

Wang, R., Cui, W., Chu, F., et al. (2020) Lithium metal anodes: Present and future. Journal of Energy Chemistry 48:145-159. https://doi.org/10.1016/j.jechem.2019.12.024, URL https://www.sciencedirect.com/science/article/pii/S2095495620300012

Wong, E., Baechle, D., Xu, K., et al. (2007) Design and processing of structural composite batteries. Proceedings of SAMPE pp 3-7

Yaji, K., Yamasaki, S., Tsushima, S., et al. (2018) Topology optimization for the design of flow fields in a redox flow battery. Structural and Multidisciplinary Optimization 57:535-546. https://doi.org/10.1007/s00158-017-1763-8

Yong-Gun, L., Satoshi, F., Changhoon, J., et al. (2020) High-energy long-cycling all-solid-state lithium metal batteries enabled by silver-carbon composite anodes. Nature Energy 5(4):299-308

Zhang, Y., Zuo, T.-T., Popovic, J., et al. (2020) Towards better li metal anodes: Challenges and strategies. Materials Today 33:56-74. https://doi.org/10.1016/j.mattod.2019.09.018, URL https://www.sciencedirect.com/science/article/pii/S1369702119307837