Image Anal Stereol 2018;37:145-158 doi:10.5566/ias.1894
Original Research Paper
MORPHOLOGICAL MODELING OF COLD SPRAY COATINGS
VINCENT BORTOLUSSI1, BRUNO FIGLIUZZIB,2, FRANÇOIS WILLOT2, MATTHIEU
FAESSEL2 AND MICHEL JEANDIN1
1Centre for Materials Engineering, Mines ParisTech - PSL Research University, Evry, F-91003, France; 2Centre
for Mathematical Morphology, Mines ParisTech - PSL Research University, Fontainebleau, 77300, France
e-mail: vincent.bortolussi@mines-paristech.fr, bruno.figliuzzi@mines-paristech.fr,
francois.willot@mines-paristech.fr, matthieu.faessel@mines-paristech.fr, michel.jeandin@mines-paristech.fr
(Received January 22, 2018; revised March 9, 2018; accepted April 13, 2018)
ABSTRACT
In this article, we study the microstructure of cold sprayed films of copper particles deposited onto a carbon
fiber reinforced polymer. The microstructure of the coating is made of a packing of seemingly round-shaped
particles of varying sizes embedded in a polymer matrix. The copper particles are separated by thin interstices.
The coating is designed to cover the body of recent commercial aircrafts. Its role is to protect the aircraft from
lightning impact by ensuring that the surface is conductive enough to evacuate electrical charges. A high
resistivity contrast is observed between the copper particles and the polymer matrix. Therefore, the global
resistivity of the material is highly dependent on the microstructure geometry.
Following an approach commonly used in materials science, to investigate its influence on the electrical
properties of the global material at the macroscopic scale, we design a 3D multiscale stochastic model that
enables us to simulate the microstructure. The model is based upon a generalization of the classical Johnson-
Mehl tessellation, which accounts for the interstices that appear between copper particles. The method is very
general and could potentially be applied to model any microstructure exhibiting similar interstices between
aggregates of particles.
Keywords: cold spray, mathematical morphology, microstructure simulation, stereology, stochastic geometry.
INTRODUCTION
Modern materials manufacturing has evolved
toward a better control and optimization at the
microscopic scale. Hence, microstructures commonly
exhibit complex morphologies mixing various
materials including carbon fibers arrangements into
polymer matrix. The physical and mechanical features
of heterogeneous materials at the macroscopic scale
are largely dictated by phenomena appearing at the
microscopic level. Therefore, the development of
numerical tools deriving the physical properties of a
material at the macroscopic scale from its geometry at
the microscopic scale has been a very active topic of
research over the past few years (Torquato, 2002).
Current computational capabilities and imaging
techniques can fill the gap between microstructure
observation and computation of effective physical and
mechanical properties through simulation. However,
if recent tomography processes can entirely recreate
complex microstructures, some materials remains
inadapted to imaging process involving high energy
beams. In this case, an alternative is to rely
on stochastic models to describe and simulate
microstructures. Computational tools and models
brought by mathematical morphology (Serra, 1982;
Soille, 2003) allow to generate and simulate complex
microstructures on which one can perform physical
simulations (Zeng et al., 2008). This method can
cope very well with multiscale and multiphased
random sets (Jeulin, 2012), as demonstrated by
recent studies investigating the effects of rigid fillers
into soft matrix such as black carbon particles
embedded into a polymer matrix (Jean et al., 2011;
Figliuzzi et al., 2016) or shells of argan nut in
polypropylene (El Moumen et al., 2014). Recently,
similar methods were used to describe polycristalyne
microstructures (Gasnier et al., 2015).
In this study, our aim is to simulate the
electrical properties of a biphased coating obtained
by thermal spraying. The coating is designed to
cover the body of recent commercial aircrafts. Its
role is to protect the aircraft from lightning impact
by ensuring that the surface is conductive enough
to evacuate electrical charges. Numerous physical
and mechanical properties of cold sprayed coatings
have already been investigated in the literature,
including oxidation (Cochelin et al., 1999), Young’s
modulus (Amsellem et al., 2008), porosity (Beauvais
et al., 2008) or thermal conductivity (Bobzin et al.,
2012). Jeandin et al. (1999) presented a model of
coating build up in plasma spray. More recently,
Delloro et al. (2017) developed a morphological
model for the cold spray process that accounts for
145
BORTULOSSI V ET AL: Morphological Modeling of cold spray coatings
the morphology of the deposit and for mechanical,
physical and dynamic phenomena. However, to our
knowledge, this study constitutes one of the first
attempt to model the electrical properties of cold spray
deposits from their microstructure.
In this paper we present a 3D stochastic model
describing the coating microstructure. The model is
derived from experimental images of slices of the
material taken with a microscope. The microstructure
is constituted of copper particles embedded in a PEEK
polymer matrix. The particles are separated by thin
interstices of PEEK polymer. In terms of resistivity,
a high constrast prevails between the polymer and the
copper particles. Hence, a key issue toward an effective
modeling of the microstructure is to accurately
reproduce the interstice thicknesses and repartition. To
that end, we introduce a new stochastic model derived
from the classical Johnson-Mehl tessellation (Johnson
and Mehl, 1939; Møller, 1989; 1992; 1994) that
enables us to properly reproduce these geometrical
characteristics. The method is very general and could
potentially be applied to model any microstructure
exhibiting similar interstices between aggregates of
particles.
The article is organized as follows. In section
“Materials”, we describe the materials on which
we conduct our study and we briefly recall the
basics of cold spray deposition. In section “Image
segmentation and characterization”, we discuss the
segmentation techniques and the statistical features
used to process the experimental microscopic images.
The 3D stochastic model is described in section
“Stochastic model for the microstructure”, and we
comment on the results of the study in section “Results
and discussion”. Conclusions are drawn in the last
section.
MATERIALS
MATERIALS
The body of recent commercial aircrafts is made
of carbon fiber-reinforced polymers. These materials
ensure high mechanical properties while being lighter
than regular aluminium alloys. However, the polymer
matrix is generally highly electrically insulating,
which causes security concerns with lightning. Higher-
grade aerospace composites are made with a matrix
of PEEK (Poly-Ether-Ether-Ketone), a thermoplastic
polymer offering good mechanical and thermal
properties. In addition the PEEK is a very good
insulator with a resistivity of 1.0×1014 Ωm.
To evacuate electrical charges in case of lightning
impact, a mesh of copper is ordinary applied onto
the composite body. Copper is a readily available
electrical conductor, that is easily machinable and
corrosion resistant. The electrical resistivity of copper
is 1.68×10−8 Ωm.
To avoid complex manufacturing and assembly
of copper meshings, we developed a new coating
method which relies on copper powder thermally
sprayed onto composite parts. More precisely, to
achieve an adherent and electrically conductive layer
we use a powders mixture. The mixture contains 80%
volumetric of spherical copper powder (10-35 µm)
and 20% of irregular PEEK particles (35-65 µm),
increasing adherence with the composite. The mixture
is sprayed using the cold-gas dynamic spraying or
“cold spray” process.
PROCESS
Cold spray is a thermal spraying process using
high speed spraying of powders to create a dense
coating. It relies on a gas flow to drive the powder
toward the substrate. The gas is usually nitrogen
under pressure (1 to 4 MPa), which is heated (200
to 800◦C) and accelerated through a convergent-
divergent (De Laval type) nozzle. The powder is
injected in the accelerated gas flow and can reach
supersonic speed. Due to relatively low spraying
temperature and processing time, powder particles
remain solid. When powder particles impact the
substrate, both of them undergo plastic strain. The
adhesion of the coating is usually guaranteed by
mechanical and chemical bonding. This process can
lead to very dense coating due to dynamic impact,
as well as low oxide content due to low temperature.
In this study, we use manufactured coating made
of a copper and PEEK mixture to investigate the
microstructure of cold spray coatings.
Samples preparation
Cold spray coatings were cross-sectioned and
polished previous to observation. Coating samples
were cut in two directions, along the spraying
path and orthogonal to the spraying path. Cutting
and polishing debond copper particles due to poor
mechanical anchorage in the matrix, leaving dark
holes at the surface. Manual polishing severely
influences debonding phenomenon. Samples were then
metallized with a 2 nm layer of gold-palladium using
a Cressington sputter coater. This is a key operation
as the layer modifies the colour of the PEEK matrix,
by greatly enhancing colour gradient between phases.
We observed cross-sections at ×20 magnification in
bright field using a Leica optical microscope with
146
Image Anal Stereol 2018;37:145-158
a resolution of 0.2428 µm per pixel. The selected
observation scale was chosen to obtain a representative
fraction of copper while highlighting PEEK interstices.
Microstructure
The microstructure of the coating can be observed
in Fig. 1 with copper particles in yellow, dark
footprints and a grey matrix of PEEK. The resulting
cold spray coating contains copper particles embedded
in a PEEK matrix. The matrix is obtained from
irregular PEEK particles highly deformed upon high
speed impact. It is considered dense with no pores
appearing at this scale. Copper particles are deforming
solely when impacting each other resulting in a limited
plastic strain. These particles are forming a network of
copper clusters.
Fig. 1: Optical microscope cross-section of the
coating microstructure with debonded particles in
black (2560×1920 pixels equals to 620×476 µm2).
Fig. 2: Magnified optical microscope cross-section of
the microstructure of the coating.
Common cold spraying of metal particles onto
metal substrate generally involve chemical bonding
and interdiffusion at the interfaces between particles
and substrate. In our case, optical observation at
larger scale highlights thin (<10 µm) PEEK layers
lying between deformed copper particles, as shown
in Fig. 2. These interstices prevent direct contact
between copper particles. Filled with electrically
insulating PEEK, they allegedly increase resistivity,
thus lowering the coating conductivity.
IMAGE SEGMENTATION AND
CHARACTERIZATION
In this section, we describe the segmentation
process developed for the images of the microstructure
of the coating.
SEGMENTATION PROCESS
The segmentation process is performed with a
Python script based upon the free image library
SMIL (Faessel and Bilodeau, 2014). It aims at
identifying the two phases of the material, namely the
PEEK matrix and the copper aggregates. Previous to
any segmentation operation, we transform the RGB-
images from the microscope into 8-bits images, by
extracting the second channel, as shown in Fig. 3a.
Then, we threshold the image with an automatic Otsu
process (Otsu, 1975) involving two thresholds. This
leads to three-phased coloured images, as shown in
Fig. 3b, where copper particles appear in yellow,
footprints of the debonded particles in black and the
PEEK matrix in red. Due to low greylevel gradient,
this first thresholding operation is not able to separate
all close particles.
In a second step, we merge the footprints into the
matrix by thresholding the particles (Fig. 3c). Next,
small isolated particles and observation artefacts with
an area smaller than 100 pixels are removed. Then we
select inner markers for the particles with by keeping
the h-maxima of a distance function applied on the
particles. This operation can be seen on Fig. 3d and
the result on Fig. 3e.
Once the markers are selected, particles are
separated with a watershed process (Beucher and
Lantuéjoul, 1979; Vincent and Soille, 1991), using
the core of the particles and a grey scale colour
gradient map to define the boundaries of each
cell. Separation achieved with the watershed process
is assessed on Fig. 3f. Finally we label each
particles obtained after the watershed. The watershed
process is able to achieve good separation of close
particles, even if the interstices do not always
147
BORTULOSSI V ET AL: Morphological Modeling of cold spray coatings
(a) (b)
(c) (d)
(e) (f)
Fig. 3: Segmentation process of particles. (a) Grey scale image. (b) Thresholded image. (c) Particles selection.
(d) Distance function on particles. (e) Cores of particles. (f) Watershed procedure on particles.
148
Image Anal Stereol 2018;37:145-158
(a) (b)
(c) (d)
Fig. 4: Segmentation process of debonded particles. (a) Thresholded image. (b) Debonded particles selection. (c)
Cores of debonded particles. (d) Watershed procedure on debonded particles.
(a) (b)
Fig. 5: (a) Final binary image. (b) Labelled particles.
149
BORTULOSSI V ET AL: Morphological Modeling of cold spray coatings
appear on thresholded images. The main steps depicted
on Fig. 3 use smaller images than the real samples.
The same exact process is applied on debonded
particles. This time real particles are merged into
the matrix, and footprints are filled by changing
pixels values of thresholded images (Fig. 4a and
Fig. 4b). As one footprint could have contained several
particles, this filling process is not able to reconstruct
hypothetical interstices separating them. Then we
apply the segmentation process used on the particles to
label the debonded particles, as shown in Fig. 4c and
Fig. 4d.
Finally we superimpose the two labelled images
(Fig. 5b). By colouring the matrix in black and
each particle in white we obtain a binary image like
Fig. 5a, used for morphological measurements. The
segmentation process was applied on 13 randomly-
selected images of 2560×1960 pixels representing an
area of 620 µm×476 µm.
MORPHOLOGICAL DATA
In this section, we conduct morphological
measurements on the binary images to characterize
the spatial distribution of copper particles. We focus
more specifically on three morphological descriptions:
the covariance, the granulometry and the thicknesses
distribution of the interstices. These morphological
measurements will be further exploited to compute the
parameters of the microstructure model.
Covariance
The covariance is defined as the probability:
C(r) = P{x ∈H ,x+ r ∈H } , (1)
where H is the union of copper particles in the
coating, x some point, and r some vector in R2.
For a stationary random set H , the covariance
does not depends on x anymore but only on r. It
depends on the orientation and on the modulus of r and
represents the probability for two points separated by
distance r to belong to the same phase of the material.
For r = 0, the covariance is equal to the surface fraction
of H (in 2D). At large distance |r|  1 the two events
x ∈H and x+r ∈H become uncorrelated and C(r =
∞)≈C(0)2.
(a)
(b)
(c)
(d)
Fig. 6: Separation of close particles. (a) and (b)
Examples of close particles. (c) Separated particles
without interstices. (d) Separated particles with
interstices.
In this study, the covariance is measured for
vectors of increasing moduli along the horizontal and
vertical directions. Note that if the microstructure
is isotropic, the covariance does not depend on the
orientation of r. In our case, isotropy and stationarity
have been assessed using all sample images. Fig. 7
shows the covariances measured on the segmented
image represented in Fig. 1.
150
Image Anal Stereol 2018;37:145-158
Fig. 7: Covariance curve as measured on the
microstructure displayed in Fig. 5a.
The covariances obtained on the sample
microstructure images show small dispersion, with
a mean surface fraction of 55% and a standard
deviation of 4.9%. We compared the mean covariances
measured vertically, horizontally and diagonally on the
segmented images.
Granulometry
The cumulative granulometry is calculated by
opening the copper phase with probability
G(s) =
P{x ∈H }−P{x ∈H (S;s)}
P{x ∈H }
, (2)
where x is a point in the image, H is the copper phase
and H (S;s) is the morphological opening of H by
a structural element S dilated by size s. In this case
the structuring element is a 8-connexity square. The
cumulative granulometry is measured on each image,
giving a mean granulometry. However due to limited
interstices segmentation, the measured granulometry
does not completely reflect the real granulometry of
the microstructure. In addition, it is important to note
that the latter is different from the initial powder’s
granulometry because of the spraying process: some
of the particles bounce at impact, thus missing in the
microstructure. In addition, particules deform upon
impact.
Fig. 8 shows the granulometry as measured on
the segmented image shown in Fig. 1. The mean
cumulative granulometry shows that 10% of the
microstructure is removed after an opening of size
1.7 µm (' 7 pixels), 50% of the microstructure is
Fig. 8: Granulometry curve as measured on Fig. 5a.
removed after an opening of size 5.7 µm (' 24 pixels)
and 90% of the microstructure is removed after an
opening of size 10.6 µm (' 44 pixels). When measured
on the thirteen samples, it is less dispersed compared
to the covariance.
Interstice thickness
Interstices appearing between copper particles are
a key feature of the microstructure that strongly impact
the conductivity of the material at the macroscopic
scale. Therefore, they must be carefully described.
The labelled image issued from the segmentation
process enables us to determine the number of non-
connected particles n0. After a dilation of 1 pixel,
subtraction of the number of remaining isolated
particles n1 to the initial quantity gives the number of
particles separated by a distance of 2 pixel: ninterstice =
n0− n1. Actually this process only measure distance
of an even number of pixels. This method is repeated
for 20 dilations on all sample images, thus yielding the
size distribution of the interstices thickness.
We apply the interstices measuring process to
obtain the thickness distribution displayed in Fig. 9.
The average distribution is fitted with an
exponential law. The probability density function of
this distribution is:
p(x,k) = ke−xk . (3)
The exponential law with k =1.19 µm−1 displayed in
Fig. 10 shows good agreement with the experimental
measurements. Interstices thickness were also
measured on magnified images, showing similar
distribution.
151
BORTULOSSI V ET AL: Morphological Modeling of cold spray coatings
(a)
(b)
Fig. 9: Interstices measurement process on a quarter
of an image. (a) Initial labelled particles. (b) Labelled
particles after 2 dilation.
Fig. 10: Interstices mean thickness distribution.
Copper phase perimeter and area
The perimeter of the copper phase can be roughly
estimated by eroding the copper particles using a
8-connexity square of size one and considering the
residue from the initial copper particles. Enclosure of
the particles remains with a thickness of one pixel. An
area measurement of the enclosure yields the perimeter
of the particles per unit area. The perimeter is very
sensitive to the segmentation of the interstices. The
area is equal to C(0), C being the covariance.
For experimental samples that have been cut
in directions parallel to the covered substrate, the
measurements yield a mean perimeter per unit area
of 0.19 µm−1 for the copper phase with a standard
deviation of 0.01 µm−1. The mean area of the copper
phase is 0.55 with a standard deviation of 0.056. For
experimental samples that have been cut in the vertical
direction, the measurements yield a mean area for the
copper phase of 0.59 with a standard deviation of
0.047.
STOCHASTIC MODEL FOR THE
MICROSTRUCTURE
In this section, we introduce a probabilistic model
for the 3D microstructures of the coating. We assume
that the microstructure is described by copper spheres
embedded in a PEEK matrix and separated by
thin PEEK interstices. Covariance, copper fraction,
granulometry, and interstice thickness are measured on
2D slices of the coating, thus 3D model’s parameters
must be inferred from these 2D informations. The
final two-scale model relies on a two-step simulation
process, namely 1/ Boolean spheres implantation of
intensity θ , whose radii follow a Gamma distribution
law with parameters λ and a; and 2/ Interstices
implantation based on a modified Johnson-Mehl
tessellation.
The boolean spheres represent copper particles,
forming clusters due to interpenetration. Simulated
volumes are sliced in 1000×1000 pixels images.
BOOLEAN MODEL OF SPHERES
We assume that the distribution of radii of the
spheres is described by a Gamma law. The probability
density function of the Gamma law is given by
p(r,λ ,a) =
ra−1
Γ(a)λ a
exp
(
− r
λ
)
, (4)
where Γ denotes the Gamma function. The average
radius of the spheres is aλ . Its variance is aλ 2.
We can easily show that the average surface of a
grain is
152
Image Anal Stereol 2018;37:145-158
Sv =
∫ +∞
0
4πra+1
Γ(a)λ a
exp
(
− r
λ
)
dr = 4πλ 2a(a+1) .
(5)
Similarly, the average volume of a grain is
Vv =
∫ +∞
0
4πra+2
3Γ(a)λ a
exp
(
− r
λ
)
dr
=
4π
3
πλ 3a(a+1)(a+2) . (6)
Let us denote by θ the intensity of the 3D Boolean
model of spheres, and by θ2 the intensity of the
corresponding Boolean model of disk as observed
on 2D slices extracted from the original model. To
determine the 3D parameters of the model from the
2D measurements, we use the stereological formulae
θVv = θ2Ā (7)
and
θSv =
4
π
θ2L̄ , (8)
where Ā is the mean area of the sliced spheres and L̄
their perimeter. We need to link the 2D measurements
to the parameters of the boolean model. To that end,
we rely on Miles’ formulae (Miles, 1972; Chiu et al.,
2013; Schneider and Weil, 2008)
AA = 1− e(−θ2Ā) , (9)
and
LA = θ2L̄(1−AA) , (10)
where Aa is the mean surface fraction of copper on the
segmented images and LA is the mean perimeter of the
copper phase on the segmented images divided by the
total surface.
Using Miles’ formulae in conjunction with
stereological formulae (Chiu et al., 2013; Schneider
and Weil, 2008), we find, for the Boolean model
AA = 1− exp(−θVv) , (11)
and
LA =
π
4
θSV exp(−θVv) . (12)
Overall, we have three unknowns in these
equations, namely the intensity θ of the Boolean
model and the parameters a and λ of the Gamma
distribution. Hence, we can express all parameters as
functions of a. Using Eqs. 11 and 12, we find
θ =− 3
4πλ 3a(a+1)(a+2)
ln(1−AA) , (13)
and
λ =− 3π
4(a+2)LA
(1−AA) ln(1−AA) . (14)
To determine the parameters of the stochastic
model, we rely on a maximum likehood approach
to find the parameters that minimize the least-square
distance between the covariance of the simulated
microstructure and the covariance that is measured
on the available experimental dataset. However LA
is highly influenced by interconnection between
particles. As many particles remains in contact due to
poor interstices segmentation, computing θ only from
λ and a provides a more robust algorithm.
INTERSTICES IMPLANTATION
To simulate the second scale of the microstructure,
that corresponds to the interstices between the particles
of the same aggregate, we rely on random Johnson-
Mehl tessellations restricted to each aggregate, or
connected component of the first scale of the
microstructure.
Let Ω denote a given volume in R3. A Voronoı̈
tessellation is a tessellation built from a Poisson point
process P in the space R3. Every point x of R3 is
associated to the class Ci containing all points of R3
closer from the point xi of P than from any other point
of P . Hence, the classes Ci, i = 1, ..,N are defined by
Ci =
{
y ∈ R3,∀ j 6= i,‖xi− y‖ ≤ ‖x j− y‖
}
. (15)
With probability one, Voronoı̈ tessellations are
normal and face-to-face. Voronoı̈ tessellations are
characterized by one single parameter, namely the
intensity of the underlying point process.
The Johnson-Mehl tessellation is a sequential
version of the Voronoı̈ model (Johnson and Mehl,
1939; Møller, 1989; 1992; 1994), where the Poisson
points are implanted sequentially according to a
parameter t which can conveniently be interpreted
as an implantation time. All classes grow then
isotropically with the same rate ν , and the growth
of crystal boundaries is stopped when they meet. All
Poisson points falling in an existing cell are removed.
From a mathematical perspective, a Johnson-Mehl
tessellation is constructed from a sequential Poisson
point process where the points xi, i = 1, . . . ,N are
implanted sequentially at a time ti, i = 1, . . . ,N. The
classes Ci, i = 1, . . . ,N corresponding to the points
xi, i = 1, . . . ,N are defined by
Ci =
{
y ∈ R3,∀ j 6= i, ti +
‖xi− y‖
v
≤ t j +
‖x j− y‖
v
}
.
(16)
153
BORTULOSSI V ET AL: Morphological Modeling of cold spray coatings
Note that when all times are set to zero, we recover the
classical Poisson-Voronoı̈ tessellation model.
In this study, to account for the interstices between
grains of the same aggregate, we rely on a modified
version of the Johnson-Mehl tessellation For each
grain n, we simulate a random number ζn according
to an exponential law with some mean k. The classes
are now defined, for a given aggregate A , by
Ci =
{
y ∈A ,∀ j 6= i, ti +
‖xi− y‖
v
+ζi ≤ t j +
‖x j− y‖
v
}
,
(17)
where ‖‖ is the geodesic distance on the set B defined
by the first scale of the model. With this definition, we
note that some points of the aggregates don’t belong
to any class of the Johnson-Mehl tessellation. We
consider that these points form the interstices between
the grains of the microstructure.
The main challenge faced using Johnson-Mehl
tesselation was to quickly calculate the position
of the interstices. To do so we rely on the fast
marching method in 2D or in 3D (Sethian, 1996). This
method solves boundary value problems of the eikonal
equation in the spatial domain Ω:
|∇u(x)|= 1
f (x)
, for x ∈Ω ,
u(x) = 0 , for x ∈ ∂Ω ,
(18)
in which u(x) is the arrival time of the front and f (x)
its propagation time. This method allows to quickly
compute propagation time of closed curves from nuclei
in space.
A question remains, which is how to select the
initial nuclei of the tessellation and the nucleation
times. While selecting the nuclei, our aim is to
preserve the geometrical shape of the grains of
the microstructure. The goal is to set the nuclei
in the center of connected components to simulate
a granulary microstructure. Therefore, we rely on
the h-maxima of the distance function to generate
the nuclei. The h-maxima of the distance function
form connected components. For each component, we
select its barycentre to be the location of a nucleus.
The threshold for the h-maxima is selected after an
optimization procedure that aims at minimizing the
distance between the granulometries. In addition, we
only keep a proportion p of the h-maxima, p being an
additional parameter of the model. For each nucleus
n, we denote by dn the value of the distance function
at the location of the nucleus. The nucleation time
associated to nucleus n is defined to be
tn = max
m
dm−dn . (19)
MAXIMUM LIKEHOOD FOR
PARAMETERS ESTIMATION
The parameters estimation follows a two-step
process. For estimating the parameters a and λ of
the first scale of the model, we compute a slice of a
3D microstructure resulting from the Boolean model.
We use 1000× 1000 pixels slices, taken from a larger
volume to account for border effects. The spheres
are implanted within this given volume, requiring it
to be large enough to reach a representative copper
fraction. We measure then the covariance associated
to the slices. The covariance is compared with the
experimental one. We aim at minimizing the least
square distance between both covariances. To do so,
we rely on a Nelder-Mead procedure optimizing on
the parameters of the radii distribution: a and λ .
Each guess of these parameters gives a value for θ
according to relation (Eq. 13). With good initial values,
it takes less than 40 iterations to fit the experimental
covariance.
For the second scale of the model, we rely
on a similar procedure, this time involving the
granulometry measurements, to estimate the threshold
value h of the h-maxima distance that is used to select
the nuclei locations, the proportion p of h-maxima
that needs to be removed and the parameter of the
exponential law describing interstice thickness.
NUMERICAL IMPLEMENTATION OF
THE MODEL
In this section, we briefly recall the main
steps of the developed algorithm and we describe
the numerical implementation of the two-scale
model. This implementation aims at simulating 3D
microstructures while reducing calculation time. The
algorithm proceeds as follows:
1. The first step of the algorithm is to simulate
a Boolean model of spheres with intensity θ .
The radii of the implanted spheres follow the
Gamma law (Eq. 4). The boolean model of
spheres is simulated using the 3D vectorial VTK-
based library vtkSim (Faessel, 2016; Faessel
and Jeulin, 2011). To accelerate volume building
from simulation, building was performed using
subvolumes. This method allows to simulate
500×500×500 voxels microstructures in 400 s
(proc time) on a PC (Ubuntu 6.14 i7 2.20 GHz
6 GB RAM). The complexity of the algorithm is
O(N2) with N the number of voxels in the volume.
2. The second step of the algorithm is the simulation
of the interstices. The algorithm follows the
following steps.
154
Image Anal Stereol 2018;37:145-158
(a) Calculation of the distance function for the
microstructure obtained from the Boolean
model of spheres
(b) Determination of the h-maxima of the distance
function for a specified h. When the h-maxima
form a connected component, the algorithm
extracts its barycenter.
(c) Random selection of nuclei from the
barycenters. The algorithm selects a given
barycenter with a specified probability p.
Each nucleus is then given a nucleation time
according to relation (Eq. 19).
(d) The propagation time of each nucleus is
computed using the fast marching method
through the Python extension module
scikit-fmm (Furtney, 2015). The algorithm
complexity of the fast marching method is
O(N lnN). The module takes as input the
coordinates of the nuclei and a velocity map
to compute the distance function. The velocity
is equal to one at the locations included in the
Boolean spheres and zero otherwise.
(e) The simulated microstructure is obtained from
the propagation times according to relation
(Eq. 17).
RESULTS AND DISCUSSION
INTERSTICES SEGMENTATION
Due to large physical difference between
PEEK and copper, the solutions for imaging the
microstructure are very restricted. In fact optical
microscopic analysis remains the best way to obtain
images of satisfactory quality. For instance, scanning
Electron Microscope would induce white artefacts
and blur around the copper particles. In addition,
in our problem, the relatively low resolution of the
images makes the selection and the segmentation
of the individual particles a rather difficult task.
In particular, it is quite difficult to achieve a good
segmentation of the PEEK interstices because of the
low image resolution and poor grey gradient at the
interfaces, so that some of them remain partially
replaced by copper. One can see on Fig. 6 that the
watershed process achieves particles separation but
is not able to include each PEEK layer between
them. Once segmented, separated particles in the
center of Fig. 6a and Fig. 6b do feature an interstice
on 6d and do not on Fig. 6c, because particles are
closer. This phenomenon could be avoided with
larger magnification. However it would signify larger
images to keep the surface representative, hence more
memory and longer processes. One must therefore
keep in mind that these phenomenon induces biased
assessment of morphological parameters. First, the
computation of the perimeter of the copper phase,
which remains highly sensitive to interstices presence,
must be used carefully. This is mostly the reason
which explains why we did not rely on Eq. 14
for the parameter identification phase. Second, the
measurements slightly underestimate the contribution
of the thinnest interstices to the thickness distribution
displayed in Fig. 10.
MODEL
To estimate the adequation between the model
and the experimental images, we considered several
morphological criteria, namely the covariance, the
granulometry and the thickness distribution of the
interstices. The joint use of the covariance and the
granulometry is relatively common in morphological
models (Jean et al., 2011). The granulometry provides
a characterization for the size distribution of the
objects that are present in the microstructure, while
the covariance provides a second order statistics which
characterize the scales present in the microstructure
as well as their superimposition. The distribution of
interstices is more specific to our model. It is a key
parameter of the model: since the conductivity contrast
is high between the copper particles and the polymer
matrix, large interstices lead to a significant decrease
of the electric conductivity of the microstructure at the
macroscopic scale.
We show in Fig. 11 a slice of the simulated
microstructure obtained with the parameters estimated
with the maximum likehood approach. The parameters
are reported in 1, with Aa the measured surface fraction
of copper on segmented images, θ the intensity of
the poisson point process, λ and a the parameters of
the Gamma distribution of the spheres radii, h the h-
maxima of the distance function used for selecting the
nuclei of the tesselation and p the proportion of h-
maxima that are randomly selected. The covariance
and the granulometry of the simulated and of the
experimental microstructures are shown in Fig. 12. A
comparison between the thickness distribution of the
interstices is displayed in Fig. 13.
Covariances and granulometries taken from
simulated and segmented images shows good
agreement, assessing the validity of the model. We
can observe a small discrepancy, around 0.01, between
the experimental and the simulated covariances for
distances ranging from 10 µm to 40 µm. We attribute
this discrepancy to the existence of a few large
155
BORTULOSSI V ET AL: Morphological Modeling of cold spray coatings
Table 1: Model parameters.
Aa θ λ a k h p
63% 4.8×10−3 µm−3 2.3 µm 0.8 1.2 µm 7.26 µm 0.5
Fig. 11: Simulated microstructure after parameters optimization (left) and experimental microstructure (right),
with resolution 0.484 µm.
Fig. 12: Simulated and experimental microstructures
comparisons between covariances (a) and
granulometries (b).
Fig. 13: Comparison between the interstice thickness
distribution, as estimated on the simulated and on the
experimental images.
zones in the polymer that do not contain copper
particles. These large zones can be observed in
Fig. 11). To accurately simulate these empty areas,
a potential solution is to introduce an additional
Boolean model aimed at defining so-called exclusion
zones in the microstructure. However, this requires
defining additional parameters for describing the
model. These parameters are difficult to estimate
from the available data: since these empty zones
are rarely observed, the experimental data does not
provide us with a statistically relevant description of
these microstructure areas. In addition, due to their
156
Image Anal Stereol 2018;37:145-158
scarcity, these zones have little influence on the electric
conductivity of the microstructure at the macroscopic
scale. Hence, since the covariance curves remain
relatively close, we discarded the existence of these
empty zones in our model.
CONCLUSION
In this article, we developed an original method
to model the microstructure of a cold spray coating.
The two-scale microstructure is modeled by Boolean
spheres separated from each other by interstices
implanted using a Johnson-Mehl tesselation. The
covariance and the granulometry measured on slices
of simulated microstructure fit well those from
segmented images. The parameters of the model are
optimized by fitting distribution law on experimental
data and with a Nelder-Mead optimization algorithm.
The model allows us to simulate representative
volumes of the microstructure. We assessed the
representativeness of the model with the covariance,
granulometry and also the distribution of interstice
size.
The robust and simple segmentation process allow
us to use a large amount of images without any
complex preparation. The model is flexible enough to
easily simulate various shape of particles or volume
fraction. Three-dimensional simulation require an fast
computation time to generate a small number of
representative volumes of coating. These simulated
microstructures will eventually allow us to investigate
the physical properties of the coating without having
to spray a large amount of various coating. This will
be done by simulations using finite element method or
Fast Fourier Transform.
ACKNOWLEDGEMENTS
This work was partly carried out within the
“C.O.MET” FUI program. The authors would like
to thank all the members of the program consortium
(under coordination of Dassault Aviation/Argenteuil-
France) and BPI France, Pôle Astech, Pôle de la
céramique, Pôle des microtechniques and DGAC for
financial support. BF and FW were supported by grant
No FA9550-15-10461 DE from the Air Force Office of
Scientific Research (AFOSR).
REFERENCES
Amsellem O, Madi K, Borit F, Jeulin D, Guipont V, Jeandin
M, Boller E, Pauchet F (2008). Two-dimensional (2D)
and three-dimensional (3D) analyses of plasma-sprayed
alumina microstructures for finite-element simulation of
Young’s modulus. J Mater Sci 43(12):4091–8.
Beauvais S, Guipont V, Jeandin M, Jeulin D, Robisson
A, Saenger R (2008). Study of the porosity in
plasma-sprayed alumina through an innovative three-
dimensional simulation of the coating buildup. Metall
Mater Trans A 39(11):2711–24.
Beucher S, Lantuéjoul C (1979). Use of watersheds in
contour detection. Proc Int Worksh Image Process Real-
Time Edge Motion Detect. Rennes, France. Sept 17–21.
Bobzin K, Kopp N, Warda T, Öte M (2012). Determination
of the effective properties of thermal spray coatings
using 2D and 3D models. J Therm Spray Techn
21(6):1269–77.
Chiu SN, Stoyan D, Kendall WS, Mecke J (2013).
Stochastic geometry and its applications. Chichester:
John Wiley & Sons.
Cochelin E, Borit F, Frot G, Jeandin M, Decker L, Jeulin
D, Al Taweel B, Michaud V, Noel P (1999). Oxidation
and particle deposition modeling in plasma spraying of
Ti-6Al-4V/SiC fiber composites. J Therm Spray Techn
8(1):117–24.
Couka E, Willot F, Jeulin D, Ben Achour M, Chesnaud A,
Thorel A (2015). Modeling of the multiscale dispersion
of nanoparticles in a hematite coating. J Nanosci
Nanotechn 15(5):3515–21.
Delloro F, Jeandin M, Jeulin D, Proudhon H, Faessel M.
Bianchi L, Meillot E, Helfen L (2017). A morphological
approach to the modeling of the cold spray process. J
Therm Spray Techn 26:1838–50.
Faessel M, Jeulin D (2011). 3D multiscale vectorial
simulations of random models. Proc Int Congr Stereol
(ICS13). Beijing, China. 19–22.
Faessel M, Bilodeau M (2014). SMIL simple morphological
image library. LRDE Tech Rep.
Faessel M (2016). vtkSim software. http://cmm.ensmp.fr/
∼faessel/vtkSim/demo/
Figliuzzi B, Jeulin D, Faessel M, Willot F, Koishi M,
Kowatari N (2016). Modelling the microstructure and
the viscoelastic behaviour of carbon black filled rubber
materials from 3D simulations. Techn Mech 36:32–56.
Furtney J (2015). Scikit-fmm software. https://github.com/
scikit-fmm
Gasnier JB, Willot F, Trumel H, Figliuzzi B, Jeulin
D, Biessy M (2015). A Fourier-based numerical
homogenization tool for an explosive material. Matér
Tech 103(3):308.
Hermann H (1991). Stochastic models of heterogeneous
materials. Zurich: Trans Tech Publ.
157
BORTULOSSI V ET AL: Morphological Modeling of cold spray coatings
Jean A, Jeulin D, Forest S, Cantournet S, N’guyen F (2011).
A multiscale microstructure model of carbon black
distribution in rubber. J Microsc 241(3):243–60.
Jeandin M, Borit F, Guipont V, Decker L, Jeulin D, Suzuki
M, Sodeoka S (1999). Lattice gas modelling in thermal
spraying. Surf Eng, 15(3):191–4.
Jeulin D (2012). Morphology and effective properties
of multi-scale random sets: A review. Compt Rend
Mécanique 340(4):219–29.
Johnson WA, Mehl RF (1939). Reaction kinetics in
processes of nucleation and growth. T Metall Soc AIME
135(8):396–415.
Miles RE (1972). The random division of space. Adv Appl
Probab 4:243–66.
El Moumen A, Kanit T, Imad A, El Minor H (2013). Effect
of overlapping inclusions on effective elastic properties
of composites. Mech Res Commun 53:24–30.
El Moumen A, Imad A, Kanit T, Hilali E, El Minor
H (2014). A multiscale approach and microstructure
design of the elastic composite behavior reinforced with
natural particles. Compos Part B Eng 66:247–54.
Møller J (1989). Random tessellations in Rd . Adv Appl
Probab 21(1):37–73.
Møller J (1992). Random Johnson-Mehl tessellations. Adv
Appl Probab 24:814–44.
Møller J (1994). Lectures on random Voronoi tessellations.
Lect Not Stat 87.
Otsu N (1975). A threshold selection method from gray-
level histograms. Automatica 11:23–7.
Schneider R, Weil W (2008). Stochastic and integral
geometry. Berlin: Springer.
Serra J (1982). Image analysis and mathematical
morphology. London: Academic Press.
Sethian JA, Popovici AM (1993). 3-D traveltime
computation using the fast marching method.
Geophysics 64(2):516–23.
Sethian JA. A fast marching level set method for
monotonically advancing fronts (1996). Proc Natl Acad
Sci 93(4):1591–5.
Soille P (2003). Morphological image analysis: principles
and applications. New York: Springer.
Torquato S (2002). Random heterogeneous materials:
microstructure and macroscopic properties. New York:
Springer.
Tsitsiklis JN (1995). Efficient algorithms for globally
optimal trajectories. IEEE T Autom Contr 40(9):1528–
38.
Vincent L, Soille P (1991). Watersheds in digital spaces:
an efficient algorithm based on immersion simulations.
IEEE T Pattern Anal 6:583–98.
Zeng QH, Yu AB, Lu GQ (2008). Multiscale modeling and
simulation of polymer nanocomposites. Prog Polym Sci
33(2):191–269.
158