# Mesoscale simulations of aqueous suspension-based electrophoretic coating process

Mélina Noblesse, Lisa Pin, Sandrine Duluard, Florence Ansart, Céline Merlet

► To cite this version:

Mélina Noblesse, Lisa Pin, Sandrine Duluard, Florence Ansart, Céline Merlet. Mesoscale simulations of aqueous suspension-based electrophoretic coating process. Open Ceramics, 2025, 21, pp.100722. 10.1016/j.oceram.2024.100722. hal-04870522

HAL Id: hal-04870522

https://cnrs.hal.science/hal-04870522v1

Submitted on 8 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mesoscale simulations of aqueous suspension-based electrophoretic coating process

Mélina Noblesse<sup>a,b</sup>, Lisa Pin<sup>b</sup>, Sandrine Duluard<sup>a</sup>, Florence Ansart<sup>a</sup>, Céline Merlet<sup>a,*</sup>

## ARTICLE INFO

**Keywords:**

Electrophoretic deposition

Mesoscopic simulations

Brownian molecular dynamics

## ABSTRACT

Electrophoretic deposition appears as a method of choice to generate coatings thanks to its ease of implementation and ability to produce deposits of relatively large thicknesses in a single step process. Research in the field has mainly focused on organic suspensions but is now moving towards suspensions in water, non-toxic and more environmentally acceptable. Here, we use Brownian molecular dynamics simulations to investigate the influence of electrophoretic deposition parameters on the properties of suspensions and deposits of mullite particles in water. We show that concentrations of particles and stabilizing ions have a large effect on inter-particle distances and electrical conductivities in the suspensions. To investigate the influence of such differences on the deposits, electrophoretic depositions with two electric fields were simulated. We demonstrate that while inter-particle distances in the deposits are much more homogenous across suspension parameters, different Debye lengths can lead to disparities in ordering (e.g. square or hexagonal).

## 1. Introduction

In the aeronautics sector, environmental issues represent a major challenge, closely linked to greenhouse gas emissions, atmospheric pollution and the consumption of natural resources. One of the main strategies for improving turbomachinery is to reduce the weight of its components. Particular attention is paid to turbine blades, usually made of nickel-based alloys. Despite their attractive properties, nickel-based alloys show their limits at high temperatures, leading to the exploration of new materials such as ceramic matrix composites (CMC) and more specifically silicon carbide (SiC) matrix and SiC fiber called SiC/SiC. While they have interesting thermomechanical properties, these materials are sensitive to water and oxygen at high temperatures [1], which can affect their durability and performance. To overcome this problem, environmental barrier coatings are developed [2]. These coatings are specifically designed to protect turbine blade materials from damage resulting from extreme operating conditions. In this context, mullite (3Al₂O₃·2SiO₂), a material resistant to corrosion by molten salts, is used to coat CMCs.

Electrophoretic deposition (EPD) is a well-established technique in the realm of wet deposition methods, offering advantages such as shaping at low temperatures, precise control of the deposited material

stoichiometry, and enhanced applicability to coat complex shape parts compared to dry processes. These advantages are valuable for various applications such as: electronic devices [3], health care [4], energy [5], and many more [6–8]. In this deposition process, an electric field is applied to a suspension of charged particles which migrate towards a selected substrate and adhere to it. By applying a coating, the substrate is shielded and its inherent properties are retained. These coatings can offer protection from chemical alterations such as corrosion, erosion, and oxidation [9–11]. EPD is also attractive due to the various coating morphologies and deposition kinetics which can be obtained depending on the electrophoretic parameters. There are different categories of parameters that can be examined: suspension parameters (zeta potential, electrophoretic mobility, particle size, presence of additives, etc.), substrate parameters (substrate shape, roughness, conductivity etc.), process parameters (electric signal shape, voltage, duration, frequency, etc.) and post-treatments (drying and sintering).

In addition to the choice of deposition material and electrodes, a critical component of the process, both ecologically and logistically, is the solvent used. Two types of solvent, organic and aqueous, can be studied. Electrophoretic deposition in organic suspensions is more commonly employed in experiments and industry due to the broad electrochemical stability window and the absence of electrode heating

<footnote><sup>a</sup> CIRIMAT, Université Toulouse 3 Paul Sabatier, Toulouse INP, CNRS, Université de Toulouse, 118 Route de Narbonne, 31062 Toulouse cedex 9, France</footnote>

<footnote><sup>b</sup> Safran Ceramics, 105 avenue Marcel Dassault, Mérignac, France</footnote>

<footnote>* Corresponding author.
E-mail address: celine.merlet@univ-tlse3.fr (C. Merlet).</footnote>
in these solvents [12]. Nevertheless, organic solvents often come with drawbacks related to their volatility, and sometimes toxicity. The utilization of water as a solvent has the potential to alleviate or eliminate environmental and toxicological concerns. However, when water is used, two phenomena can arise at the electrodes as a result of water electrolysis: i) local pH fluctuations that can destabilize the suspension and ii) the formation of bubbles at the electrodes that can introduce porosity and/or degrade the coating. It is worth noting that in some aqueous suspensions, it is possible to realize EPD without detecting any water electrolysis or bubble formation. A wide variety of electrophoretic parameters can be used to tune the topography of the coatings, e.g. to obtain rough/smooth, dense/porous, homogeneous/heterogeneous coatings. However, the link between process parameters and deposit characteristics such as homogeneity, microstructure, thickness, is poorly understood. Numerical tools can provide useful information to understand the deposition mechanisms and predict the properties of the coating as a function of the selected set of parameters. The effect of such phenomena on the mesoscopic structure of the deposit can only be investigated through mesoscopic simulations.

Currently, several modeling techniques, focusing on different scales (atomic, particle, cell, module and system) are being used to investigate phenomena related to electrophoresis. Vermet et al., used a Finite Element Method to study macroscopically the coverage of the electrophoretic paint coating in the automotive sector [13–15]. Their aim was to develop a method for studying deposit coverage in order to predict possible corrosion sites. More recently, Salazar de Troya et al. also used a Finite Element Method to study flow-based electrophoretic deposition and understand specific morphologies observed in such conditions [16]. However, such macroscopic models do not consider certain key experimental parameters such as the number of particles, the zeta potential, or the granulometry of the particles. Giera et al. investigated the organic deposition of monodispersed polystyrene using Brownian molecular dynamics [17]. The study showcased the evolution of deposit thickness over deposition time, reaching 2 µm within a 2 ms timeframe. Furthermore, the research delves into the examination of material organization influenced by the Debye length and Peclet number. More recently, Giera et al. proposed the inclusion of local physical phenomena in their particle-based model, such as the immobilization of particles in the deposit or the charge transfer between the particles and the electrode substrate [18].

In this work, as a first step towards the simulation of electrophoretic deposition of industrially relevant particles from aqueous suspensions, we adapt the model developed by B. Giera et al. [17] which focused on polystyrene particles in organic suspensions, to the case of mullite particles, from environmental barrier coatings in water. The parameters are adjusted on experimental values when available. This allows us to investigate the effect of the concentration of particles, the concentration of stabilizing ions and the applied voltage on the microstructure of the deposit. We determine pair distribution functions and particle electrical conductivities for the various suspensions studied and show that the temperature does not affect the inter-particle distance but has a large influence on the dynamics. Such a detailed simulation study of the suspensions has not yet been reported in the literature. We then characterize deposits through pair distribution functions and coordination numbers and demonstrate that the concentration of particles and the applied electric field can have some influence on the compactness of the deposit. Different organizations of particles are observed, namely hexagonal and square, with more or less defects.

## 2. Methods

Brownian molecular dynamics simulations have been realized on two distinct systems: (a) suspensions of charged particles with 3D periodic boundary conditions and (b) the same suspensions undergoing electrophoretic depositions with 2D periodic boundary conditions (in x and y directions). All simulations have been conducted with the

LAMMPS software (Large-scale Atomic/Molecular Massively Parallel Simulator [19]). The simulations are conducted following several steps illustrated in Fig. 1: i) generation of initial positions for the particles in a random organization, ii) relaxation of the particles positions in the case of deposits, iii) production of the trajectories of the particles through Brownian molecular dynamics. For suspensions, the analyses are done on part of the trajectory ensuring that the system is equilibrated. For deposits, only the deposits at the end of the simulations are characterized.

To investigate the effect of the concentration of particles on the dynamical and structural properties of the suspensions as well as their effect on the final deposit, three concentrations were considered. Suspensions with weight percentages equal to 3, 5 and 8, corresponding to suspensions which can be prepared experimentally, were simulated. The volume of the simulation box, constant during the simulations and the same for all systems, is equal to 7 × 7 × 13.5 µm³. The corresponding numbers of particles and volume fractions are given in Table S1.

### 2.1. Interaction potentials

Following the work of Giera et al. [17], the overall energy of the system is described as the sum of an internal potential, $U_{internal}$, and an external potential, $U_{external}$:

$$U(r^{col}) = \sum_{i=1}^{N_{col}} \sum_{j>i} U_{internal}(r_{ij}) + \sum_{i=1}^{N_{col}} U_{external}(z_i) \quad (1)$$

where $N_{col}$ is the total number of colloidal particles in the system, $r_{ij}$ is the distance between particles i and j, and $z_i$ is the distance between particle i and the substrate.

In the model chosen here, the solvent is considered in an implicit way and $U_{internal}$ is composed of two terms: $U_{DLVO}$ for the interactions between colloid particles and $U_{solvent}$ for the interactions between the particles and the solvent. Eq. (2) identifies each term of the internal potential, where $\eta$ is the viscosity and $T$ is the temperature.

$$U_{internal}(r) = U_{solvent}(r, \eta, T) + U_{DLVO}(r) \quad (2)$$

The presence of an implicit solvent is introduced through the consideration of hydrodynamic and Brownian forces, the latter entailing a random and disorganized movement. The hydrodynamic forces are introduced through the isotropic terms of the Fast Lubrication Dynamics method [20]. Details regarding these specific forces can be found in the LAMMPS documentation [17]. It is worth noting that the use of an explicit solvent representation is prohibited here as the number of water molecules needed would be on the order of 10¹³, orders of magnitude larger than the largest recent molecular simulations [21–23].

Following the Derjaguin, Landau, Verwey, and Overbeek (DLVO) theory [24,25], which accounts for the overall interaction potential between two particles, $U_{DLVO}$ is composed of a van der Waals potential $U_{VDW}$ and a repulsive potential $U_{Yuk}$ stemming from electrostatic forces.

$$U_{DLVO}(r) = U_{Yuk}(r, a, \lambda_D, A_{Yuk}) + U_{VDW}(r, a, \sigma_{LJ}, A_{col}) \quad (3)$$

The Yukawa potential $U_{Yuk}$, describing the electrostatic double layer repulsion, depends on the Debye length $\lambda_D$, as defined by Eq. (5), the particle radius $a$, and the Yukawa force prefactor defined in Eq. (6).

$$U_{Yuk} = A_{Yuk} \lambda_D \exp\left(\frac{r - 2a}{\lambda_D}\right) \quad (4)$$

$$\lambda_D = \sqrt{\frac{\epsilon\epsilon_0 k_B T}{(Ze)^2 \rho_{\pm}}} \quad (5)$$

$$A_{Yuk} = 32\pi \frac{a\epsilon\epsilon_0}{\lambda_D} \left(\frac{k_B T}{Ze}\right)^2 \tanh^2\left(\frac{\zeta Ze}{k_B T}\right) \quad (6)$$

In the above equations, $\epsilon_0$ is the vacuum permittivity, $\epsilon$ is the relative
(a) Suspensions

<figure><img src="image_2.png" /></figure>

PARTICLES GENERATION:
random positions

<figure><img src="image_3.png" /></figure>

PARTICLES MOTION:
Brownian dynamics

(b) Electrophoretic deposition

<figure><img src="image_4.png" /></figure>

PARTICLES GENERATION:
random positions in a smaller box

<figure><img src="image_5.png" /></figure>

SUSPENSION STABILIZATION:
relaxation of the positions

<figure><img src="image_6.png" /></figure>

PARTICLES MOTION:
Brownian dynamics under E

Fig. 1. Steps followed for the simulation of (a) suspensions and (b) electrophoretic deposition. For the EPD simulations, the positions of the particles are first generated in a slightly smaller box to allow for the inclusion of the substrate.

permittivity of the solvent, e the elementary charge, k_B the Boltzmann constant, ρ_± corresponds to the bulk concentration of small ions within the suspension (e.g. OH⁻) and Z their valency (here considered equal to 1), ζ is the zeta potential of the colloid particles, q_col is their charge and σ_LJ stands for the Lennard-Jones diameter. The other term of the DLVO potential, U_VDW, follows a commonly used expression [26] and contains both attractive and repulsive terms.

$$U_{\text{VdW}}(\mathbf{r}) = U_{\text{repulsive}}(\mathbf{r}, a, \sigma_{\text{LJ}}, A_{\text{col}}) + U_{\text{attractive}}(\mathbf{r}, a, A_{\text{col}}) \quad (7)$$

External potentials, only present when there is a substrate, correspond to the interactions between particles and the substrate, and the application of an electric field.

$$U_{\text{external}}(z^{N_{\text{col}}}) = \sum_{i=1}^{N_{\text{col}}} U_{\text{wall}}(z_i, \sigma_{\text{LJ}}, A_{\text{wall}}) + U_{\text{Field}}(z_i) \quad (10)$$

The substrate is modelled through a van der Waals potential where the Hamaker constant and the Lennard-Jones diameter are taken equal to the values for the colloid particles. The electrophoretic deposition is modeled by the application of a uniform and constant field, E_field, acting

$$U_{\text{repulsive}}(\mathbf{r}) = A_{\text{col}} \sigma_{\text{LJ}}^6 \left\{ \frac{[16a^6(525r^{10} + 13696a^8r^2 - 10528a^6r^4 + 4536a^4r^6 - 840a^2r^8 - 7680a^{10})]}{[4725r^8(r^2 - 4a^2)^7]} \right\} \quad (8)$$

$$U_{\text{attractive}}(\mathbf{r}) = \frac{A_{\text{col}}}{6} \left( \ln\left(1 - \frac{4a^2}{r^2}\right) - \frac{8a^4 - 4a^2r^2}{r^4 - 4a^2r^2} \right) \quad (9)$$

on colloids with an effective charge q_col determined through Eq. (12) [27]. The electric field is oriented along the z axis, in the direction perpendicular to the substrate.

$$U_{\text{Field}}(z) = E_{\text{Field}}q_{\text{col}}z \quad (11)$$

$$q_{\text{col}} = \frac{a\pi\epsilon_0\epsilon_z}{24\lambda_D^6} \left( a^5\lambda_D - a^4\lambda_D^2 - 10a^3\lambda_D^3 + 6a^2\lambda_D^4 + 96\lambda_D^6 + a^4(12\lambda_D^2 - a^2)e^{a/\lambda_D} \int_1^\infty \frac{e^{-ag/\lambda_D}}{\hat{g}} d\hat{g} \right) \quad (12)$$

Where A_col is the Hamaker constant which is a fundamental physical constant associated with the strength of van der Waals interactions.
<caption>**Table 1**</caption>

Parameters defined directly following known experimental quantities. The different zeta potentials correspond to the different concentrations of particles.

<table><thead><tr><th>Variable</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Concentration</td><td>[wt%]</td><td>3; 5; 8</td><td>wt%</td></tr><tr><td>Zeta potential</td><td>ζ</td><td>-60.8; -65.1; -72.0</td><td>mV</td></tr><tr><td>Radius</td><td>a</td><td>85</td><td>nm</td></tr><tr><td>Density</td><td>d</td><td>3.1</td><td>g cm<sup>-3</sup></td></tr><tr><td>Permittivity</td><td>ε</td><td>80</td><td>n.u</td></tr><tr><td>Viscosity</td><td>η</td><td>1.0, 0.8, 0.7</td><td>mPa s</td></tr><tr><td>Temperature</td><td>T</td><td>290, 300, 310</td><td>K</td></tr></tbody></table>

## 2.2. Parametrization of the interaction potentials

Some of the parameters considered in the model are known experimentally and set up accordingly. Table 1 summarizes the quantities used in this case. No specific experimental details are given in this work as the aim here is to assess trends rather than reaching an absolute comparison with experiments. The Hamaker constant, $A_{col}$, is not known precisely for mullite particles but a value of $10^{-20}$ J, similar to values measured for other particles [28], can be considered reasonable. The Lennard Jones radius $\sigma_{L,J}$ is defined as 0.3a following Giera et al. [17].

While the pH of the aqueous suspension is known and close to 9, the concentration of small ions cannot be taken equal to $10^{-5}$ mol L$^{-1}$ as other ions are added to stabilize the suspension. The exact concentration of small ions is not known and several concentrations are considered. As the concentration of small ions has a large influence on electrostatic screening, these concentrations lead to quite different Debye lengths and resulting forms of the interaction potentials. The various Debye lengths obtained for different concentrations of small ions and temperatures are summarized in Table 2. The Yukawa prefactors and total charges of the particles determined for the various concentrations in particles and small ions, and the different temperatures are all given in Table S2.

The inter-particle interaction potential, $U_{DLVO}$, resulting from the model and parametrization chosen is shown in Fig. 2. The van der Waals and Yukawa contributions are given in Figure S1. The van der Waals potential decays to zero for relatively small distances. A cut-off of 6a, beyond which $U_{vdW}$ is not computed and considered to be zero, is chosen. Different cut-off values are chosen for the Yukawa potential depending on the concentration of small ions: $20\lambda_D$ for $5.10^{-3}$ mol L$^{-1}$, $40\lambda_D$ for $5.10^{-4}$ mol L$^{-1}$, and $50\lambda_D$ for $5.10^{-5}$ mol L$^{-1}$. While the temperature appears in the equations describing the Debye length and the Yukawa prefactor, its effect on $U_{DLVO}$ is very limited here as can be seen in Figure S2.

The inter-particle interaction potential has a very different shape for the largest concentration of small ions, corresponding to the smallest Debye length. This suggests a different behavior in terms of organization of particles in the suspension, and potentially in the deposit.

## 2.3. Characterization of suspensions and deposits

The organization of the particles in suspensions or in the deposits is analyzed using pair distribution functions which quantify the probability of encountering a pair of particles i and j separated by a distance $r_{ij}$, in comparison to the probability anticipated in a fully random

<caption>**Table 2**</caption>

Debye lengths calculated for different concentrations of small ions and temperatures considered.

<table><thead><tr><th>[mol L<sup>-1</sup>]</th><th>λ<sub>D</sub> (nm) at 290K</th><th>λ<sub>D</sub> (nm) at 300K</th><th>λ<sub>D</sub> (nm) at 310K</th><th>Notation</th></tr></thead><tbody><tr><td>5 10<sup>-3</sup></td><td>6.1</td><td>6.2</td><td>6.3</td><td>6 nm</td></tr><tr><td>5 10<sup>-4</sup></td><td>19.2</td><td>19.5</td><td>19.8</td><td>20 nm</td></tr><tr><td>5 10<sup>-5</sup></td><td>60.6</td><td>61.7</td><td>62.7</td><td>60 nm</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 2.** Inter-particle interaction potential, $U_{DLVO}$, as a function of inter-particle distance $r$ for the various concentrations of small ions considered.</figcaption></figure>

distribution with an equivalent density. Pair distribution functions are calculated following:

$$g(r) = \frac{1}{N_{col}^2} \sum_{i=1}^{N_{col}} \sum_{j=1}^{N_{col}} \langle \delta|r_{ij}| - r \rangle \quad (13)$$

where $r_{ij}$ is the distance between particles i and j and r is the distance considered. According to this definition, peaks in the pair distribution functions indicate probable inter-particle distances.

The electrical conductivity, due to particles motion, was determined for all suspensions according to:

$$\sigma = \frac{e^2}{k_B T V} \lim_{t \to \infty} \frac{1}{6t} \left\langle \left| \sum_{i=1}^{N_{col}} q_{col} \delta r_i(t)^2 \right| \right\rangle \quad (14)$$

where V is the volume of the simulation box, $\delta r_i(t)^2$ is the mean-square displacement of particle i, t is the time and $\langle \rangle$ denotes an average over all colloidal particles [29]. This conductivity, which takes into account correlated motions of charged particles, will be designated as the “real” conductivity. The Nernst-Einstein conductivity, which assumes independent motions of particles, was also calculated following:

<figure><img src="image_8.png" /><figcaption>**Fig. 3.** Pair distribution functions at 300 K for the 3wt%, 5wt% and 8wt% suspensions simulated. r is the inter-particle distance.</figcaption></figure>
$$ \sigma_{NE} = \frac{e^2}{k_B T V} N_{col} q_{col}^2 D_p \quad (15) $$

where $D_p$ is the individual diffusion coefficient of the particles determined from mean-square displacements:

$$ D_p = \lim_{t \to \infty} \frac{1}{6t} \langle |\delta r_i(t)|^2 \rangle \quad (16) $$

# 3. Results and discussions

## 3.1. Suspensions

To investigate how the concentration of particles and small stabilizing ions influence the structure of the suspension, and in particular inter-particle distances, pair distribution functions have been obtained for 3wt%, 5wt% and 8wt%, and for Debye lengths of 6 nm, 20 nm, and 60 nm corresponding to various concentrations of small ions (see Table 2). Such a characterization of the distribution of particles in suspension has, to the best of our knowledge, not been reported so far. Fig. 3 shows pair distribution functions calculated for a temperature of 300K.

The pair distribution functions for $\lambda_D$ equal to 60 nm show a markedly different behavior with a main peak at distances varying with the concentration of particles and subsequent fluctuations (see Figure S3 for larger inter-particles distances). In this case, the DLVO potential becomes strongly repulsive at distances close to 800 nm (much larger than the particle size) and the particles organize themselves in a quite ordered fashion, as far from each other as possible, according to their concentration. The most probable distance thus decreases dramatically with

increasing concentration from 635 nm at 3wt% to 485 nm at 8wt%.

The different behaviors identified in the pair distribution functions can be visualized in Fig. 4 showing regions of the 3wt% and 8wt% suspensions. For short Debye lengths, particles form clusters and appear to agglomerate while for a Debye length of 60 nm the suspensions are more homogeneous. Indeed, groups of particles at close distances and relative large voids are visible in the extracted configurations of the simulations with short Debye lengths.

Fig. 5a details the evolution of the most probable inter-particle distance, $r_{peak}$, as a function of the Debye length for all concentrations of particles and temperatures considered. The influence of the Debye length, identified in the pair distributions functions at 300 K, is confirmed for all temperatures. At short Debye lengths, $r_{peak}$ values are similar regardless of mass percentage. In contrast, at high Debye lengths, an effect of mass percentage is observed. The effect of the temperature on $r_{peak}$ is very limited for all conditions.

Fig. 5b shows the electrical conductivity due to the particles motion for all suspensions at the different temperatures studied. The mean-square displacements, from which the conductivities are determined, are shown in Figure S4. It is worth noting that the values obtained are much smaller than what is measured in typical aqueous suspensions (~50–70 µS cm⁻¹). In fact, measurements correspond to the overall electrical conductivity due to both particles and small ions, which move fast, leading to much higher experimental values.

The conductivity increases with the concentration of particles. This is expected as the higher the concentration, the higher the number of charges carriers and thus the conductivity. Unsurprisingly, the conductivity also increases with temperature, in accordance with an

<figure><img src="image_9.png" /><figcaption>Fig. 4. Snapshots of a region of the simulation box for the all suspensions studied: particles with $z < 500$ nm, full section of the box in x and y - 7000 × 7000 nm².</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 5. (a) Most probable inter-particle distance, r<sub>peak</sub>, as a function of the temperature. (b) Real conductivity, σ, as a function of the temperature for 3wt%, 5wt% and 8wt%.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Comparison of real (see Eq. (14)) and Nernst-Einstein (see Eq. (15)) conductivities, due to particles mobility, for 3wt% (left panel, black) and 8wt% (right panel, blue) suspensions with various Debye lengths and at various concentrations.</figcaption></figure>

increase in particles mobility accentuated by the decrease in viscosity. Interestingly, the conductivity decreases with an increasing Debye length.

To check that the increase of conductivity with concentration of particles and temperature is simply due to the number of charge carriers and the increased mobility, the Nernst-Einstein conductivity, which is a function of the number of particles and diffusion coefficients (see Eq. (15)) was determined. While relatively common in the study of electrolytes, such an analysis is rarely conducted on suspensions of particles. Indeed, in many materials such as glasses, and solid and liquid electrolytes, the Haven ratio, which is the ratio between the real and Nernst-Einstein conductivity, has shown large deviations from unity usually indicating correlations between motions of charged species [30–32]. The comparison between both conductivities, for 3wt% and 8wt%, is plotted in Fig. 6 (values for 5wt% are shown in Figure S5).

For small concentrations and Debye lengths, the agreement between real and Nernst-Einstein conductivities is very good showing that the motion of the particles is uncorrelated. This indicates that the light agglomeration in these cases is reversible. For larger concentrations or Debye lengths, the real conductivity is larger than the Nernst-Einstein conductivity showing that the particles move in a much more correlated fashion. The effect of the Debye length can be clearly seen in **Figure S6**. In particular, for 8wt% and a Debye length of 60 nm, the motion of the particles is minimal but the conductivity is still significant. This suggests that particles do small motions but all in the same direction at a given time. This is consistent with a system with a strong repulsion where all particles interact strongly and each particle moving pushes nearby particles.

The characterization of the various suspensions studied show significant differences in terms of structural and dynamical properties. In particular, the most probable inter-particle distance can go from 250 nm, for a Debye length of 6 nm, to 635 nm, for a Debye length of 60 nm

and a concentration of 3wt%. Some results also indicate a light reversible agglomeration for some of the suspensions. In addition, the electrical conductivities span a large range of values from 0.5 to almost 6 µS cm⁻¹. Such variations can influence the electrophoretic deposition and this is what will be explored next.

## 3.2. Coatings

Simulations of electrophoretic depositions were realized at 300 K for all suspensions and two values of applied electric field (3700 and 7400 V cm⁻¹). These values are not very comparable with those used experimentally, but experimental values do not allow for reasonable calculation times. Electric fields 100 times larger have therefore been applied to obtain acceptable calculation times, in agreement with previous studies [17]. To get a better idea of the influence of the applied electric field, some simulations have been conducted with applied electric fields of 370 and 740 V cm⁻¹. The results suggest that while the compaction is slightly smaller for lower electric fields, the influence of the suspension parameters is the same as for the larger electric fields. Fig. 7 shows side and top views of the deposits obtained for all Debye lengths for the 8wt% suspension. The first layer, i.e. the particles directly in contact with the substrate, and the second layer, i.e. the particles in the layer slightly further away from the substrate, are displayed. Visualizations of the deposits for smaller concentrations of particles are shown in **Figure S7**. As expected, the particles are much more densely packed in the deposit compared to the suspensions (see Fig. 4). The particles seem to be mostly organized following a hexagonal structure but regions with a square packing can be observed, especially for the intermediate Debye length of 20 nm. Such ordering of coatings deposited through electrophoresis have been observed both from simulations [17] and experiments [33, 34] with various suspensions.

Domains with different orientations of the hexagonal or square
<figure><img src="image_12.png" /><figcaption>Fig. 7. Side and top views of the deposit obtained with a 7400 V cm⁻¹ electric field for each Debye length for the 8wt% suspension.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 8. Number of particles in the first two layers of the deposit as a function of the Debye length for all suspensions simulated and an applied electric field of (a) 3700 V cm⁻¹ and (b) 7400 V cm⁻¹.</figcaption></figure>

structures are visible in all cases for the 8wt% suspension. At the boundaries, defects are observed where particles are undercoordinated. For the 3wt% suspension, the deposit seems less dense, especially for the second layer of particles. In the less dense layers, the structure looks more homogeneous and with less defects. It is worth noting that the number of the particles in the system is quite limited, leading to a full depletion at the end of deposition. In the case of the 3wt% suspension, the structure of the deposit is probably largely influenced by the number

of particles. Indeed, in this case, there are not enough particles to generate three layers of deposits.

To compare the compactness of the deposits obtained for different concentrations of particles with 3700 or 7400 V cm⁻¹ electric fields applied, the number of particles in the two first layers of deposit was determined. Results are shown in Fig. 8. For all suspensions, the number of particles is lower in the second layer than in the first layer. However, the difference between layers decreases with the increase in the
concentration of particles. This is probably due to a stronger pressure exerted by additional layers being deposited when there are more particles in the system. It is worth reminding here that the electric field is constant in the simulations (see Eq. (11)). A different potential drop could exist with deposits of varied densities, this is out of the scope of this work and will be investigated in the future.

On Fig. 8, one can also see that the number of particles decreases with the increase in Debye length. This is consistent with an increase in repulsion between colloidal particles when the concentration in small ions decreases. For large concentrations of particles and strong electric fields, the influence of the Debye length and the proximity to the substrate decrease. Interestingly, the highest density of particles, for 8wt%, correspond to a high coverage of the substrate, around 60%.

To confirm the influence of the “outside” layers, i.e. far from the substrate, on the “inner” layers, i.e. in contact or close to the substrate, the number of particles in different regions of the deposit are determined in two more cases. For the 3wt% suspension case, for which the depletion in particle is strong, we conducted additional simulations where the initial suspension is made twice larger in the z direction (corresponding to around 5000 particles). These simulations, designated as “with z extension”, are done for all Debye lengths considered in this work and results are shown in Figure S8. The addition of particles leads to the formation of a deposit with five layers instead of three. The first and second layers are indeed denser when more particles are present. Nevertheless, the qualitative trend observed with the change in Debye length is the same for the original and extended systems.

For the 8wt% suspension case, it is possible to look at the number of particles in five layers going from direct contact with the substrate to the end of the deposit. Results of this analysis are given in Figure S9. It is clear that in most cases, the compactness decreases when the distance from the substrate increases. For a Debye length of 6 nm, the variation is less dramatic than for larger Debye lengths as it is favorable to have particles in close proximity.

Going beyond the characterization of the density of particles in the deposit, it is possible to characterize its microstructure more precisely. Fig. 9 shows the pair distribution functions for the particles in the first layer of the deposits obtained with an applied electric field of 3700 V cm⁻¹ (results for 7400 V cm⁻¹ are shown in Figure S10). Well defined peaks are observed for all concentrations of particles and Debye lengths. The most probable inter-particle distance is close to 200 nm, much smaller than in most of the suspensions. Most deposits show additional smaller peaks at larger distances indicating an ordering at intermediate range, in agreement with the front views of the deposits in Figs. 7 and S7.

For the smallest Debye length of 6 nm, neighboring particles adopt a similar distance of 205–215 nm for all deposits. For larger Debye lengths, the most probable inter-particle distance depends on the concentration of particles in the suspensions. Larger concentrations leading to closer packing. For an applied electric field of 7400 V cm⁻¹, most inter-particles distances are the same or slightly smaller than at 3700 V cm⁻¹.

Overall, the results show that the electric fields applied are sufficient to overcome the repulsion between particles and lead to a dense deposit. In addition, in terms of inter-particle distances, the application of an electric field reduces the differences between suspensions with distinct Debye lengths or concentration of particles. It is worth reminding here that the applied electric fields are quite high. In the future, simulations with smaller electric fields, still compatible with the observation of a deposition on the time scales of calculations should be explored.

As described previously from Figs. 7 and S7 showing front views of first and second layers of the deposits, the particles can adopt a hexagonal or square structure with domains having different orientation leading to different quantities of boundaries and defects. To assess the prevalence of structures depending on the conditions of deposition, distributions of coordination numbers were calculated for all deposits. These are shown in Figs. 10 and S11. Coordination numbers of equal to 6 indicate a hexagonal structure while coordination numbers equal to 4 suggest a square organization or zones of low density of particles.

<figure><img src="image_14.png" /><figcaption>**Fig. 9.** Pair distribution functions for the first layer of the deposits simulated with an applied electric field of 3700 V cm⁻¹: for (a) the 3wt% suspension, (b) the 5wt% suspension and (c) the 8wt% suspension. (d) Most probable inter-particle distance in the first and second layers for all deposits obtained with an applied electric field of 3700 V cm⁻¹.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 10.** Distributions of coordination numbers in the first (a and c) and second (b and d) layers of the deposits obtained for 3wt% (a and b) and 8wt% (c and d) suspensions.</figcaption></figure>

Coordination numbers equal to 5 suggest defects and boundaries.

Fig. 10 confirms that in all deposits, the coordination numbers mainly observed are 4, 5 and 6 as expected. In rare cases, a few particles with 3 or 7 neighboring particles are present. In most cases, the most probable coordination state is 6, as can also be seen on the front views of the deposits. The second layer usually shows less particles with a coordination 6, and more with a coordination 4 or 5, than the first layer. This is consistent with a lower coverage and with more defects in the second layer than in the first layer.

While no straightforward trend is observed, it seems that the fraction of particles with a coordination 4 increases with the concentration of particles. The combination of a small (6 nm) or intermediate (20 nm) Debye length leads to the largest fractions of coordination 4 in dense layers. The presence of a square structure can be identified in the pair distribution functions in such cases. Indeed, a small peak at around 300 nm, corresponding to the next nearest neighbors in the square structure, can be identified in Figs. 9.c, S10b and S10c.

Overall, in terms of particle packing, it seems that the differences of structural and dynamical properties identified in the various suspensions have little effect. This suggests that, as long as the concentration of particles and electric fields are high enough, experiments with different concentrations of particles and small ions, and different applied electric fields could lead to similar microstructures.

## 4. Conclusions

Following the recent interest for electrophoretic deposition in aqueous suspensions, a less toxic and more environmentally acceptable alternative to organic suspensions, we have used Brownian molecular dynamics simulations to investigate the influence of some process parameters on the structural and dynamical properties of suspensions as well as on the microstructure of the deposits. To this aim, we have adapted a previously published model to the case of mullite particles in water.

We have conducted simulations of suspensions with different concentrations of particles and different concentrations of stabilizing ions at various temperatures and demonstrated that the suspensions obtained have a wide range of inter-particle distances and electrical conductivities. In particular, the most probable inter-particle distance can go from 250 nm, for the smallest Debye length of 6 nm, to 635 nm, for the largest Debye length of 60 nm and a concentration of 3wt%. In some cases, a light reversible agglomeration is observed, consistent with an energy well in the DLVO interaction potential. The good agreement between the real and Nernst-Einstein electrical conductivities for small concentrations and Debye lengths suggest an uncorrelated motion of the particles. In contrast, for larger concentrations or Debye lengths, the real conductivity is larger than the Nernst-Einstein conductivity which suggests a correlated motion between particles, consistent with the strong repulsion existing in these cases. Such a detailed characterization of the suspensions of mullite particles has not been previously reported in the literature.

We then conducted simulations of deposition on a smooth substrate under applied electric fields of 3700 and 7400 V cm⁻¹. These electric fields appear sufficient to overcome the particle-particle repulsion and to obtain dense deposits with surface coverages reaching 60% close to the substrate. While significant differences were observed between structural and dynamical properties of suspensions, no straightforward trends can be extracted for the microstructure of the deposits.

The current study focused on a few suspension and process parameters. In the future, it would be interesting to investigate the effect of more complex properties and phenomena such as the polydispersity of the particle sizes, the roughness of the substrate or the sometimes observed electrolysis at the electrodes.

## Supporting information

Supplementary Information available: total number of particles in the different simulation boxes, additional model parameters for the
various suspensions studied, plots of interaction potentials not shown in the main text, pair distribution functions in suspensions and deposits, mean square displacements, electrical conductivities, visualizations of suspensions and deposits, number of particles and distributions of coordination numbers.

## Author contributions

CM, SD, and FA supervised and guided the project. All authors designed the research. MN prepared the input files for the simulations and conducted all calculations and analyses. CM and MN produced the necessary software, when missing, for conducting the analyses. CM and MN drafted the manuscript. All authors participated in the interpretation of the results and in the manuscript revision.

## Funding sources

The authors thank the national research and technology agency (ANRT) as well as the SAFRAN Ceramics company for supporting part of this research work.

## CRediT authorship contribution statement

**Mélina Noblesse:** Writing – review & editing, Writing – original draft, Software, Methodology, Investigation, Formal analysis. **Lisa Pin:** Writing – review & editing, Resources, Project administration, Funding acquisition, Conceptualization. **Sandrine Duluard:** Writing – review & editing, Supervision, Investigation, Funding acquisition, Conceptualization. **Florence Ansart:** Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization. **Céline Merlet:** Writing – review & editing, Writing – original draft, Supervision, Software, Methodology, Investigation, Formal analysis, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The authors acknowledge Thibaud Derieux, Remi Martin for useful discussions. This work was granted access to the HPC resources of the CALMIP supercomputer center under the allocation P21014.

## Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.o ceram.2024.100722.

## References

[1] E.J. Opila, H. R, Paralinear oxidation of CVD SiC in water vapor, J. Am. Ceram. Soc. 80 (1997) 197.

[2] H. Eaton, G.D. Linsey, K.L. More, J.B. Kimmel, J.R. Price, N. Miriyala, EBC protection of SiC/SiC composites in the gas turbine combustion environment, Turbo Expo: Power Land, Sea Air 4 (2000). V004T02A018.

[3] S. Yokoyama, I. Suzuki, K. Motomiya, T. Hideyuki, K. Tohji, Aqueous electrophoretic deposition of citric-acid-stabilized copper nanoparticles, Colloid Surf. 545 (2017) 93.

[4] T. Yoshioka, A. Chavez-Valdez, J.A. Roether, D.W. Schubert, A.R. Boccaccini, AC electrophoretic deposition of organic-inorganic composite coatings, J. Colloid Interf. Sci. 392 (2013) 167.

[5] R.N. Basu, C.A. Randall, M.J. Mayo, Fabrication of dense zirconia electrolyte films for tubular solid oxide fuel cells by electrophoretic deposition, J. Am. Ceram. Soc. 84 (2001) 33.

[6] L. Besra, A review on fundamentals and applications of electrophoretic deposition (EPD), Prog. Mater. Sci. 52 (2007) 1.

[7] P. Amrollahi, J.S. Krasinski, R. Vaidyanathan, L. Tayebi, D. Vashaee, Electrophoretic deposition (EPD): fundamentals and applications from nano- to microscale structures. Handbook of Nanoelectrochemistry, Springer, 2015, p. 561.

[8] A.R. Boccaccini, I. Zhitomirsky, Application of electrophoretic and electrolytic deposition techniques in ceramics processing, Curr. Opin. Solid State Mater. Sci. 6 (2002) 251.

[9] M. Prioux, S. Duluard, F. Ansart, G. Pujol, P. Gomez, L. Pin, Advances in the control of electrophoretic process parameters to tune the ytterbium disilicate coatings microstructures, J. Am. Ceram. Soc. 103 (2020) 6724.

[10] J. Wagner, C. Josse, L. Gani, S. Knittel, P.-L. Taberna, F. Ansart, Electrophoretic deposition of aluminum particles from pure propan-2-ol suspensions, Res. Mater. 13 (2022) 100259.

[11] C.A. Patricio Magalhaes, F. Ansart, P.-L. Taberna, J.-P. Bonino, Electrophoretic deposition of hybrid film on aluminum 2024 using sol-gel boehmite nanoparticles, Surf. Coat. Tech. 289 (2016) 165.

[12] M. Ammam, Electrophoretic deposition under modulated electric fields: a review, RSC. Adv. 2 (2012) 7633.

[13] K. Verma, H. Cao, P. Mandapalli, R. Wille, Modeling and simulation of electrophoretic deposition coatings, J. Comput. Sci. 41 (2020) 101075.

[14] K. Verma, L. Ayuso, R. Wille, Parallel simulation of electrophoretic deposition for industrial automotive applications, in: International Conference on High Performance Computing & Simulation, 2018, p. 468.

[15] K. Verma, J. Oder, R. Wille, Simulating industrial electrophoretic deposition on distributed memory architectures, in: 27th Euromicro International Conference on Parallel Distributed and Network-Based Processing, 2019, p. 414.

[16] M.A. Salazar de Troya, J.R. Morales, B. Giera, A.J. Pascall, M.A. Worsley, R. Landingham, W.L. Du Frane, J.D. Kuntz, Modeling flow-based electrophoretic deposition for functionally graded materials, Mater. Des. 209 (2021) 110000.

[17] B. Giera, L.A. Zepeda-Ruiz, A.J. Pascall, T.H. Weisgraber, Mesoscale particle-based model of electrophoretic deposition, Langmuir. 33 (2017) 652-661.

[18] J.J. Karnes, A.-J. Pascall, C. Rehbock, V. Ramesh, M.A. Worsley, S. Barcikowski, E. Lee, B. Giera, Particle-based simulations of electrophoretic deposition with adaptive physics models, Comput. Phys. Commun. 297 (2024) 109062.

[19] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1.

[20] A. Kumar, J.J.L. Higdon, Origins of the anomalous stress behavior in charged colloidal suspensions under shear, Phys. Rev. E. 82 (2010) 051401.

[21] X. Zhou, Z. Wei, H. Lu, J. He, Y. Gao, X. Hu, C. Wang, Y. Dong, H. Liu, Large-scale molecular dynamics simulation based on heterogeneous many-core architecture, J. Chem. Inf. Model. 64 (2024) 851.

[22] J. Jung, C. Kobayashi, K. Kasahara, C. Tan, A. Kuroda, K. Minami, S. Ishiduki, T. Nishiki, H. Inoue, Y. Ishikawa, M. Feig, Y. Sugita, New parallel computing algorithm of molecular dynamics for extremely huge scale biological systems, Comput. Chem. 42 (2020) 231.

[23] J. Jung, C. Tan, Y. Sugita, GENESIS CGDYN: large-scale coarse-grained MD simulations with dynamic load blancing for heterogeneous biomolecular systems, Nature Commun. 15 (2024) 3370.

[24] B. Derjaguin, L. Landau, The theory of stability of highly charged lyophobic sols and coalescence of highly charged particles in electrolyte solutions, Acta Physicochim. URSS 14 (1941) 633.

[25] E. Verwey, G. Overbeek, Theory of the Stability of Lyophobic Colloids, Elsevier, Amsterdam, 1948.

[26] R. Everaers, M.R. Ejtehadi, Interaction potentials for soft and hard ellipsoids, Phys. Rev. E 67 (2003) 041710.

[27] B. Giera, L.A. Zepeda-Ruiz, A.J. Pascall, J.D. Kuntz, C.M. Spadaccini, T. H. Weisgraber, Mesoscale particle-based model of electrophoresis, J. Electrochem. Soc. 162 (2015) D3030.

[28] S.-L. Tsaur, R.M. Fitch, Preparation and properties of polystyrene model colloids: II. Effect of surface charge density on coagulation behavior, J. Colloid Interf. Sci. 115 (1987) 463.

[29] D. Frenkel, B. Smit, Understanding Molecular Simulation, from Algorithms to Applications, Academic Press, 2002.

[30] G.E. Murch, The haven ratio in fast ionic conductors, Solid. State Ion. 7 (1982) 177.

[31] J.O. Isard, The Haven ratio in glasses, J. Non. Cryst. Solids. 246 (1999) 16.

[32] B.R.N.M. Vargas-Barbosa, Dynamic ion correlations in solid and liquid electrolytes: how do they affect charge and mass transport? ChemElectroChem. 7 (2019) 367.

[33] M. Böhmer, In situ observation of 2-dimensional clustering during electrophoretic deposition, Langmuir. 12 (1996) 5747.

[34] M. Trau, D. Saville, I. Aksay, Field-induced layering of colloidal crystals, Science (1979) 272 (1996) 706.