Numerical and analytical studies of
kinetics, equilibrium, and stability of
the chemical reaction fronts in
deformable solids
vorgelegt von
Aleksandr Morozov, M.Eng.
von der FakultΓ€t V β Verkehrs- und Maschinensysteme
der Technischen UniversitΓ€t Berlin
zur Erlangung des akademischen Grades
Doktor der Ingenieurwissenschaften
β Dr.-Ing. β
genehmigte Dissertation
Promotionsausschuss:
Vorsitzender: Prof. Dr.-Ing. U. von Wagner
Gutachter: Prof. Dr. rer. nat. W. H. MΓΌller
Prof. Dr. Sc. A. B. Freidin
Prof. Dr. hab. V. A. Eremeev
Prof. Dr.-Ing. habil. Dr. h. c. mult. H. Altenbach
Tag der wissenschaftlichen Aussprache: 22. Dezember 2020
Berlin 2021
Numerical and analytical studies of the chemical reaction front kinetics in solids III
Abstract
In the present work a chemical reaction between a solid and a diffusing con-
stituent is considered. Many experimental observations show a coupling be-
tween mechanical stresses and the kinetics at a chemical reaction interface that
propagates in a deformable solid. In general, two major processes control the
chemical reaction: diffusion of a reactant through the reacted material and
the consumption of the diffusing reactant by a chemical reaction. Mechanical
stresses may affect the reaction front propagation via the direct influence on
the diffusion process, the reaction rate, or both.
The reaction localized at a sharp interface is considered. The kinetics of
the interface is modeled by using the chemical affinity tensor concept, which,
in a thermodynamically consistent way, couples mechanical stresses and the
chemical reaction rate. The tensorial nature of the affinity tensor follows from
the fact that the reaction occurs at the oriented surface element.
The chemical interface can be retarded or even blocked by mechanical stresses.
When modeling the moving interface, especially when approaching the blocking
position, the question of stability arises naturally. A kinetic stability approach
for the analytical stability analysis of interfaces at a chemical equilibrium is
formulated and applied in several examples.
A finite element procedure for numerical simulations of the reaction front
propagation is developed based on remeshing algorithms. In order to check
as to whether the numerical procedure can be used for the stability analysis
of the interface, several problems are examined numerically, and the results
are verified with analytical solutions and predictions of stability. Also, the
developed procedure is cross-validated with another numerical method for
simulating the moving interface, namely CutFEM.
The last part of the manuscript reports experimental findings of IMC growth
in microchips during a high-temperature storage test. The growth kinetics is
modeled by applying a continuum model based on the chemical affinity tensor.
By evaluating and combining a real experiment with theory, values for the
diffusion coefficient and the chemical reaction constants are obtained. These
values demonstrate the consistency of the developed theoretical models.
Numerical and analytical studies of the chemical reaction front kinetics in solids V
Zusammenfassung
In dieser Arbeit wird die chemische Reaktion zwischen einem FestkΓΆrper und
einem diffundierenden Anteil untersucht. Viele experimentelle Beobachtungen
zeigen eine Kopplung der mechanischen Spannungen mit der Kinetik an der
chemischen ReaktionsflΓ€che auf, die sich in einem verformbaren FestkΓΆrper
ausbreiten kann. HauptsΓ€chlich steuern zwei Prozesse die chemische Reaktion:
die Diffusion eines Reaktionspartners innerhalb des reagierenden Materials und
die Umsetzung des diffundierenden Reaktionspartners aufgrund der chemis-
chen Reaktion. Die mechanischen Spannungen kΓΆnnen somit eine direkte
Auswirkung auf die Ausbreitung der Reaktionsfront haben, indem sie entweder
den Diffusionsprozess beeinflussen, oder indem sie die Reaktionsgeschwindigkeit
manipulieren, wobei auch beides gleichzeitig mΓΆglich ist.
Der Reaktionsvorgang an einer ausgezeichneten Grenzschicht wird unter-
sucht. Die Kinetik der Grenzschicht wird mit Hilfe des Konzeptes des chemis-
chen AffinitΓ€tstensors modelliert, der, aufbauend auf den thermodynamischen
Grundgesetzen, die mechanischen Spannungen mit der chemischen Reaktions-
geschwindigkeit in Verbindung setzt. Der AffinitΓ€tstensors orientiert sich am
entsprechenden OberflΓ€chenelement, an dem der Reaktionsvorgang stattfindet.
Die Ausbreitung der chemischen Grenzschicht kann durch die mechanische
Beanspruchungen behindert oder sogar komplett blockiert werden. Bei der
Modellierung der sich bewegenden Grenzschicht, insbesondere beim Erreichen
der Blockade, stellt sich die Frage nach der StabilitΓ€t. Ein kinetischer Ansatz
fΓΌr die analytische StabilitΓ€tsuntersuchung von Grenzschichten im chemischen
Gleichgewicht wird formuliert und in mehreren Beispielen angewendet.
Ein Finite-Elemente-Verfahren zur numerischen Simulation der Ausbreitung
der Reaktionsfront wird mit Hilfe von Algorithmen zur Wiedervernetzung
entwickelt. Um zu prΓΌfen, ob das numerische Verfahren fΓΌr die StabilitΓ€ts-
analyse an den Grenzschichten hinreichend genau ist, werden verschiedene
Problemstellungen auch numerisch untersucht. Diese Ergebnisse werden mit
den analytischen LΓΆsungen und den StabilitΓ€tsvorhersagen bestΓ€tigt. ZusΓ€t-
zlich wird das entwickelte Verfahren zur Simulation der Ausbreitung der sich
bewegenden Grenzschicht auch mit einer anderen numerischen Methode, der
sogenannten CutFEM, validiert.
Im letzten Teil der Arbeit werden experimentelle Neuheiten des IMC-
Wachstums in Mikrochips wΓ€hrend eines Hochtemperatur-Storage-Tests
vorgestellt. Die Kinetik des Wachstums wird mit Hilfe eines Kontinuumsmod-
VI
ells, das auf dem chemischen AffinitΓ€tstensor basiert, nachgestellt. Durch
die Auswertung und Kombination von realen Experimenten mit der Theorie,
werden der Diffusionskoeffizient und die Konstanten der chemischen Reak-
tion genau bestimmt. Diese Werte zeigen die Konsistenz der entwickelten
theoretischen Modelle.
Numerical and analytical studies of the chemical reaction front kinetics in solids VII
Preface
The research reported in this doctoral dissertation has been carried out in the
Technical University of Berlin, Institute of Mechanics, Chair of Continuum
Mechanics and Constitutive Theory, during 2017β 2020. The work was sup-
ported by DFG/RFFI grants (No. MU 1752/47-1 and 17-51-12055) on topic
βMechanochemistry of advanced anode designs in Li-ion batteriesβ.
I would like to thank both my supervisors, Professor Wolfgang H. MΓΌller
and Professor Alexander Freidin, for allowing me to perform research in this
fascinating interdisciplinary field of science. I am grateful for their guidance
and the kind support which they gave me during the last years. They created a
working environment in which I felt comfortable and could fully concentrate on
my research. Further, I would like to express my gratitude to Professor Victor
Eremeev and Professor Holm Altenbach for accepting my work for review and
providing valuable feedback, and to Professor Utz von Wagner for being a
chairman at the dissertation procedure.
The main numerical results of this work were obtained during my extremely
productive research visits to the Aalto University and to the University of
Warwick. For these opportunities and warm hosting, I want to thank Professor
Jarkko Niiranen and Professor Εukasz Figiel. Many significant improvements for
my numerical procedure were investigated during our brainstorming discussions
with Dr. Mikhail Poluektov, for which I am very grateful. I am also thankful to
Dr. Alexander Semencha, whose help in organizing and running the experiments
cannot be overestimated.
Special gratitude should be addressed to my friends and colleagues Dr.
Viacheslav Balobanov, Dr. Sergei Khakalo for their inspiration, curiosity, and
fruitful discussions about this research, the science in general, and far beyond
these topics. Also, I would like to thank the group of Professor MΓΌller at
the TU-Berlin, especially M.Sc. Wilhelm Rickert, M.Sc. Gregor Ganzosch,
Arion Juritza, and Miriam Ziert, who helped me a lot to organize my work
and private life in Germany.
And last but not least, I want to thank all my family for the belief and
support and my dear Dr. Galina Lavrenteva, who was waiting for me so long.
Numerical and analytical studies of the chemical reaction front kinetics in solids IX
Contents
Abstract III
Zusammenfassung V
Preface VII
1 Introduction 1
2 Chemical reaction front kinetics and stability 9
2.1 Chemical affinity tensor and general problem statement . . . 9
2.2 On stability of a chemical interface . . . . . . . . . . . . . . . 16
2.2.1 Variation of the normal . . . . . . . . . . . . . . . . . 17
2.2.2 Variation of a jump across the interface . . . . . . . . 18
2.2.3
Linearized kinetic equation for perturbed phase interfaces
20
2.2.4
Linearized kinetic equation for perturbed chemical in-
terfaces .......................... 22
2.3 Stability of a planar chemical interface . . . . . . . . . . . . . 23
2.4 Stability of a cylindrical chemical interface . . . . . . . . . . . 30
2.5 Remark on phase transition zones . . . . . . . . . . . . . . . . 37
2.6 Conclusions............................ 39
3 Numerical simulation and study of the chemical reaction front
kinetics 43
3.1 Description of the numerical procedure . . . . . . . . . . . . . 44
3.2
Kinetics of the cylindrical chemical reaction front approaching
theblockingstate......................... 48
3.3 A note about stresses caused by the interface instability . . . 51
3.4 Comparison of the CutFEM and remeshing results . . . . . . 55
3.4.1 CutFEM overview . . . . . . . . . . . . . . . . . . . . 55
3.4.2
Stable and unstable configurations for planar chemical
interface.......................... 56
3.4.3
Stable configuration for the chemical interface under
shear............................ 59
3.4.4 Closed interface in the square domain . . . . . . . . . 60
3.5
Kinetics of a chemical reaction front in a body with a pore or
aninclusion............................ 63
Contents
X
3.6 Conclusions............................ 66
4 Experimental and theoretical studies of Cu-Sn intermetallic phase
growth 69
4.1 Overview on intermetallic compound growth . . . . . . . . . . 69
4.2 Experiment overview . . . . . . . . . . . . . . . . . . . . . . . 72
4.2.1 Specimen preparation . . . . . . . . . . . . . . . . . . 72
4.2.2 Experimental procedure . . . . . . . . . . . . . . . . . 73
4.3 Experimental results . . . . . . . . . . . . . . . . . . . . . . . 74
4.3.1 Chemical composition of the compounds . . . . . . . . 75
4.3.2 Thickness data for the IMC layers . . . . . . . . . . . 77
4.4 Theoretical model . . . . . . . . . . . . . . . . . . . . . . . . 82
4.4.1 Analytical solution of a model problem . . . . . . . . 82
4.4.2 Fitting the model parameters . . . . . . . . . . . . . . 86
4.5 Conclusions............................ 90
5 Conclusions 91
List of Figures XI
List of Tables XV
Bibliography XVII
A Analytical solutions for planar and cylindrical reaction fronts kiξΏ
netics i
A.1 Planar reaction front kinetics . . . . . . . . . . . . . . . . . . i
A.2 Cylindrical reaction front kinetics . . . . . . . . . . . . . . . . iii
Contents
Numerical and analytical studies of the chemical reaction front kinetics in solids 1
1 Introduction
During the last decades, much attention in continuum mechanics has been paid
to the study of materials, which change their structure due to phase or chemical
transformations under thermomechanical actions (see, e.g., Grinfeld [1991];
Gurtin [2000]; Abeyaratne and Knowles [2006]; MΓΌller, Vilchevskaya, and
Freidin [2015]; Freidin and Vilchevskaya [2020] and the references therein). The
peculiarity of such studies is their interdisciplinarity, where coupled problems
of mechanics, physics, and chemistry arise. As an example of mechanical and
physical phenomena coupling, one can consider austenite-martensite phase
transformation or orientation (ordering) transformation in polymers. In the
present work, coupled problems of mechanics and chemistry are considered.
There are two terms to describe such a coupling: mechanochemistry and
chemomechanics (see, e.g., Freidin and Vilchevskaya [2020]). Both of them
used to emphasize that mechanical processes influence the chemical ones
and vice versa. In this work, stress-affected chemical reactions in solids are
considered. Without further elaboration on etymology, the coupling of chemical
and mechanical phenomena is defined by the word βmechanochemistryβ in the
following text.
Some examples of mechanochemical problems come from MicroElectroMe-
chanical Systems (MEMS) and microelectronics. Oxidation processes intercon-
nected with crack growth in polycrystalline silicon microscale parts determine
the lifetime of MEMS (e.g., Muhlstein, Stach, and Ritchie [2001]; Muhlstein,
Brown, and Ritchie [2002]; Muhlstein and Ritchie [2003]). In Buttner and
Zacharias [2006] it was shown that on nanoscale mechanical stresses can retard
the oxidation of silicon nanowires.
In microelectronics, the growth of the Intermetallic Compound (IMC) phases
in Pb-free solders is of great interest. The main technological process for
creating an electrical contact between components in a (micro-) electric circuit
is soldering. After soldering, an IMC layer appears and establishes a mechanical
contact between eutectic tin-silver solder bumps and Cu interconnects, see,
e.g., Callister Jr. and Rethwisch [2010]. IMC formation is a result of diffusion
and chemical reaction processes, which involve a change in shape and volume
between the products and reactants. Strictly speaking, IMC formation is based
on a multicomponent diffusion in solids, including vacancies as a migrating
species leading to Kirkendall voiding (Dybkov [2010]; Ross, Vuorinen, and
Paulasto-KrΓΆckel [2016]). In addition to mechanical stress it can be enhanced
2
by electric currents (Chao et al. [2009]) Fig. 1.1. Consequently, the rate of
IMC growth has a strong implication on solder joint reliability.
(a) (b)
Fig. 1.1:
Cross section of a solder bump prior (a) and after (b) current stressing,
reprinted from Chao et al. [2007]
Other examples of such mechanochemical studies are Lithium-Ion Batteries
(LIBs) with novel anode materials, such as silicon. The charge capacity of a
silicon-based anode is about ten times higher than the capacity of graphite
anodes, which are most commonly used in the industry (Kasavajjula, Wang,
and Appleby [2007]). However, the lithiation and de-lithiation with large
amounts of Li result in dramatic volumetric changes of 300% in the Si anodes,
which can lead to fracture of the anode, Fig. 1.2. In order to cope with the large
(a) (b)
Fig. 1.2:
SEM morphology of 250 nm a-Si film on Cu after 1 (a), and 30 cycles (b),
reprinted from Kasavajjula, Wang, and Appleby [2007]
volume change and, therefore, to obtain better capacity retention and cycle life
for Si anodes, various designs of different anode structures have been proposed,
such as thin films (Bourderau, Brousse, and Schleich [1999]), nanoparticles (Liu
et al. [2012a]), nanowires or hollow nanowires (Jia and Li [2015]; Zhao et al.
[2012b]), and even remarkable morphologies, such as honeycomb structures
(Baggetto, Danilov, and Notten [2011]). Unless accommodated by appropriate
compensating deformation, the large volumetric change generates mechanical
Chapter 1. Introduction
Numerical and analytical studies of the chemical reaction front kinetics in solids 3
stresses, which, in turn, may affect the diffusion (Yang [2010], Chang, Moon,
and Cho [2015]). Recent experimental observation shows (Liu et al. [2012a];
Liu et al. [2012b]; McDowell et al. [2013]; Liu et al. [2013]) that there is a sharp
interface between lithiated and unlithiated phases, that the lithiation process
may be controlled by an interfacial chemical reaction (Zenga et al. [2016]), and
that the chemical reaction rate is affected by mechanical stresses as well (Yang
[2010]; Cui, Gao, and Qu [2012a]; Levitas and Hamed [2013]). As a result, one
faces a complex coupled problem of mechanochemistry with a moving chemical
reaction front.
A distinctive feature of the aforementioned examples is that the thickness
of the interface is negligible when compared to the characteristic dimensions
of considered solid bodies. For instance, in the case of silicon lithiation, the
interface thickness is roughly 1nm when the silicon nanoparticle or nanowire
has a diameter of 150 -250 nm, see, e.g., McDowell et al. [2013].
In general, two major processes control the propagation of a sharp chemical
reaction front: (i) diffusion of a reactant in the body undergoing a chemical
reaction, (ii) consumption of the diffusing reactant by a chemical reaction at
the reaction front. Thus, stresses may affect the reaction front propagation
via the influence on diffusion, or the reaction rate, or both. To couple the
chemical reaction rate with the mechanical stresses, some models include
additional stress-dependent cross-effect terms in the diffusion flux which appears
based on generalized expressions of stress-dependent scalar chemical potentials:
Knyazeva [2003]; Loeffel and Anand [2011]; Cui, Gao, and Qu [2012b]; Bower
and Guduru [2012]; Levitas and Hamed [2013]; Brassart and Suo [2013].
However, the velocity of the reaction front may be controlled by the reaction
rate rather than by the diffusion (see, e.g., Zhao et al. [2012a]; Jia and
Li [2015] and references therein). In this case, the influence of stresses on
the reaction rate becomes important. In classical physical chemistry, the
reaction rate is determined by a scalar chemical affinity, which is equal to a
combination of scalar chemical potentials of the reaction constituents (Prigogine
and Defay [1954]). Scalar chemical potentials were derived for the case of phase
transformations in gases and liquids where stresses were reduced to a scalar
pressure (Gibbs [1948]).
In the last decades of the XXth century, it was recognized, that chemical
potential is a tensor in the case of solid phases. The tensorial nature followed
from the fact that the equilibrium was considered not just in a point of the
interface but at oriented surface element, see, e.g., Bowen [1967]; Truesdell
[1969]; Grinfeld [1980]; James [1981]. Later, a tensorial nature of the chemical
potential was discussed in Rusanov [2005]; Rusanov [2006]. These results
could be considered as a prelude to the tensorial chemical affinity, and in
Freidin [2013]; Freidin, Vilchevskaya, and Korolev [2014] it was shown, that
in the case of a propagating reaction front the driving force is equal to the
4
normal component of the chemical affinity tensor which, in turn, equals to the
combination of the chemical potential tensors of solid constituents and the
chemical potential of the diffusing constituent. By this approach, mechanical
stresses affect the reaction rate via the chemical affinity tensor. A kinetic
equation for the propagating chemical reaction front was formulated in the
form of the dependence of the normal component of the reaction front velocity
on the normal component of the chemical affinity tensor.
This work can be considered as a step of further utilization of the chemical
affinity tensor concept in a form developed in Freidin [2009]; Freidin [2013];
Freidin, Vilchevskaya, and Korolev [2014]; Freidin [2015]; Freidin et al. [2016],
for studying stress-controlled chemical reaction front propagation. The kinetics
of the front propagation is described within the frames of the mechanics of
configurational forces. It is known that in the case of quasi-static stress-
induced phase transformations the configurational force driving the interface
is given by the jump of the normal component of the Eshelby stress tensor
(see, e.g., monographs Silhavy [1997]; Abeyaratne and Knowles [2006]; Maugin
[2010]; Gurtin [2000] and reference therein). As it was derived and approbated
during the last years, in the case of stress-assisted chemical reactions, the
configurational force is determined by the normal component of a chemical
affinity tensor. In both cases, based on irreversible thermodynamics reasons, a
kinetic equation can be formulated in the form of the dependence of the interface
velocity on the corresponding configurational force. Then stresses and strains
affect the interface propagation as they are presented in the configurational
force.
In the case of phase transformation, a zero configurational force corresponds
to phase equilibrium, but this is only a necessary condition of phase equilibrium,
and the interface may be unstable in such a state (Grinfeld [1982]; Grinfeld
[1990]; Grinfeld [1991]). In the case of a chemical reaction, zero configura-
tional force corresponds to chemical equilibrium at which the reaction front
propagation is blocked. Therefore, by analogy with phase transformations, the
question about the stability of the reaction front in the vicinity of the blocking
state arises naturally. Chemical reaction retardation and even blocking by
mechanical stresses were experimentally observed in, e.g., Marcus [1982]; Kao
et al. [1988]; Mihalyi, Jaccodine, and Delph [1999]; Heidemeyer et al. [2000];
Buttner and Zacharias [2006]. Special attention was paid to the influence of
mechanical stresses on the stability of the growing interface in Ortiz, Repetto,
and Si [1999]; Phan et al. [2001]; Barvosa-Carter and Aziz [2004]; Zeeshan and
Venkatasubramanian [2017]; Ahmad and Viswanathan [2017]. These problems
occur in the solid electrolyte interfaces modeling for Li-ion batteries, Natsiavas
et al. [2016]; HΓΌter et al. [2017].
The importance of the stability analysis of chemical reaction fronts also
follows from the following circumstances. The difference between phase and
Chapter 1. Introduction
Numerical and analytical studies of the chemical reaction front kinetics in solids 5
chemical transformations is that in the latter case one deals with an open
system with diffusion. One more equation β a diffusion equation β is included in
the analysis in comparison with phase transformations. But a more important
difference may be that a stress-induced phase transformation starts only if
it provides minimization of the energy, and the origin of the transformation
may be inside of a body or at the outer surface of a body depending on energy
preferences and stability. For example, in the case of isotropic phases, for a
spherical particle undergoing a phase transformation under increasing external
strain, an equilibrium and stable new phase domain can grow only from the
center of the particle if a shear module of a new phase is larger than a shear
module of an initial phase. Other interfaces are unstable even if they satisfy
phase equilibrium conditions. If the shear module of the new phase is less than
the shear module of the initial phase, then a stable new phase grows from the
outer surface of the particle (Eremeev, Freidin, and Sharipova [2003]; Eremeev,
Freidin, and Sharipova [2007]). The origin of phase transformations can also
be a new phase nucleus, which starts at imperfections or heterogeneities (see,
e.g., Vilchevskaya, Filippov, and Freidin [2013]). Then the position of the
interfaces is affected by the position of the imperfections. Moreover, two-phase
microstructures may be more preferential then two-phase configurations with
one smooth interface (Ball and James [1987]; Chenchiah and Bhattacharya
[2008]; Antimonov, Cherkaev, and Freidin [2016]; Freidin and Sharipova [2018]).
At given boundary conditions, it is a material βthat decidesβ where to initiate
a new phase and how it will grow.
The situation is different in the case of chemical reaction fronts. Here the
direct reaction can go only if a diffusing reactant is delivered to the reaction
front through the transformed material. In the case of the aforementioned
spherical particle, the supply of the diffusing constituent is possible only
through the surface of the particle. Therefore, the direct chemical reaction
can start only at the outer surface and the front can propagate only inwards,
irrespectively of the shear moduli of the solid reactants. In the case of the
reverse reaction, the reactant is taken away from the front, and the front
propagates outwards.
Therefore, starting place of the reaction and direction of the front propagation
are controlled by the supply or removal of the diffusing reactant. It means that
the reaction front may be forced to go in the direction of an unstable blocking
state and this instability may affect the propagation of the interface. Thus,
in the case of chemical reactions unstable reaction fronts are of interest not
only from the mathematical point of view but also because they are physically
meaningful.
Stability of the interface between material phases for phase transformation
was previously studied in, e.g., Grinfeld [1982]; Gurtin [1983]; Grinfeld [1990];
Grinfeld [1991]; Fried [1993]; Osmolovsky [2000]. In the current work, the
6
stability of the reaction front is studied by using the procedure of so-called
linear kinetic stability analysis and says that the interface is unstable if its
small perturbations grow due to the kinetic equations. The utilized procedure
was developed and explored earlier for the case of equilibrium phase interfaces
in Eremeev, Freidin, and Sharipova [2003]; Fu and Freidin [2004]; Freidin et al.
[2006]; Eremeev, Freidin, and Sharipova [2007].
Note that kinetic stability analysis performs not only an energy-based check-
ing of the stability and states the fact of the stability loss, but also gives
hints on types (or modes) of instabilities formation and the tendencies of
further kinetics of the perturbations. In Phan et al. [2001]; Barvosa-Carter
and Aziz [2004] the authors numerically considered the growth stability of
crystalline silicon from the amorphous phase. Based on a kinetic model from
Barvosa-Carter et al. [1998], the authors considered the influence of mechanical
stresses and growth kinetic anisotropy on the interface roughening. Consid-
ering the kinetic model of film growth resulting from deposition and mass
transport proposed in Ortiz, Repetto, and Si [1999], the authors in Natsiavas
et al. [2016] analyzed the influence of the elastic pre-stress on the stability of
planar growth asymptotic analysis of a nearly flat interface for the problem of
anode-electrolyte interaction. Based on these results, the authors in HΓΌter et al.
[2017] considered the stability of the electrode-coating-electrolyte interface
depending on the thickness of the thin film interlayer and the magnitude of
the elastic pre-stresses.
In Zeeshan and Venkatasubramanian [2017] the stability criteria for elec-
trodeposition at solid-solid interfaces is derived using linear stability analysis
assuming that the solids are linearly elastic isotropic materials based on a
kinetic model proposed by Monroe and Newman [2004]. It described later
cathodic roughening and dendritic growth by taking into account mechanical
stresses and their influence on the current exchange densities and potentials at
the roughened interface.
Independently of the kinetic stability analysis, the concept of phase transition
zones (PTZ) was introduced in Freidin and Chiskis [1994a]; Freidin and Chiskis
[1994b]. The PTZ is a zone in a strain space formed by all strains which can
exist at the equilibrium interfaces in a given material. The procedure was
developed for the case of nonlinear-elastic materials in Freidin and Chiskis
[1994a]; Freidin and Chiskis [1994b]; Freidin et al. [2006], and then for the
case of linear elastic phases in Freidin [1999]; Morozov and Freidin [1998].
Correlation of the results obtained by the kinetic stability analysis and the
PTZ approach is discussed in Fu and Freidin [2004]; Eremeev, Freidin, and
Sharipova [2003]; Freidin et al. [2006]; Eremeev, Freidin, and Sharipova [2007];
Vilchevskaya, Filippov, and Freidin [2013]. In these works, it was shown that
instability of phase interface is not observed if the strains at the interface in a
body belong to the boundary of the PTZ. Then in the papers Grabovsky and
Chapter 1. Introduction
Numerical and analytical studies of the chemical reaction front kinetics in solids 7
Truskinovsky [2011]; Grabovsky and Truskinovsky [2013], it was proved that
belonging of strains to the external PTZ-boundaries is a necessary stability
condition.
In the current work, the kinetic stability analysis and the PTZ procedure
are compared for the interfaces at the equilibrium position. However, the PTZ
approach cannot be explicitly extrapolated to the case of the chemical interface
at the reaction blocking state. The linear stability analysis procedure for phase
transformations is extended to the case of a chemical reaction. Note that
analytical investigation of the interface stability based on the perturbed kinetic
equation is rather complicated, even for simple geometries. Therefore, in this
work, a numerical procedure for simulating the reaction front propagation is
developed. This procedure is verified with analytical solutions, and special
attention is paid to the interfaces approaching the equilibrium position. It is
done to check whether the developed numerical methods can reveal physically
stable and unstable configurations.
Many numerical simulations of the transformation fronts propagation have
been made earlier in problems of different nature. For the case of phase
transformations see, e.g., Finite Element Analysis (FEA) applications in Mueller
and Gross [1998]; Mueller and Gross [1999]; Gross, Mueller, and Kolling [2002];
Mueller, Gross, and Lupascu [2006]. Stress-field analysis in Li-ion batteries
during propagation of the lithiation front propagation was performed using FE
in Jia and Li [2015]. A description of the FE simulation of the propagation
of a chemical reaction front with the reaction front kinetics controlled by the
chemical affinity tensor was implemented in Freidin et al. [2016].
In this work, an FEM procedure for numerical simulations of the reaction
front propagation is developed based on remeshing algorithms. Utilizing
the proposed procedures, the propagation of the transformation fronts is
simulated, and the consistency between the numerical and analytical results
is checked. The method itself is compared to other numerical approaches
to simulate the interface propagation, namely the CutFEM approach. A
more detailed discussion about this is given in Chapter 3. The numerical
simulations demonstrate how interfaces propagate if the equilibrium positions,
found analytically, are stable or unstable. In the case of stable configuration, the
interface converges smoothly to the equilibrium position. If the configuration
is unstable, the interface also propagates smoothly toward the equilibrium
position, but in the vicinity of the equilibrium, when the interface velocity
is almost zero, the instabilities become visible and start to grow. This effect
may be explained by the competition between global and local kinetics of
the interface propagation. It should be emphasized that the existence and
propagation of the chemical reaction fronts far from equilibrium is a natural
process and stress-induced front retardation and these instabilities, which might
8
grow considerably, may be the source of further fracture, or another type of
failure.
This manuscript is organized as follows. In Chapter 2, an overview of the
chemical affinity tensor concept is given. A kinetic stability approach for the
analytical stability analysis of interfaces at a chemical equilibrium is formulated.
It is done by extending the linear stability analysis procedure developed for
phase transformations. In this chapter two simplest problems are studied for
linear elastic solid phases of reactants: stability of the plane problem with
planar interface in the infinite layer and axially-symmetric plane problem
with a cylindrical interface. The influence of the material parameters on the
stability of the interface is studied, and, as a result, βstableβ and βunstableβ
sets of parameters are proposed for the following numerical simulations. An
FE procedure for numerical simulations of the reaction front propagation based
on remeshining algorithms is introduced in Chapter 3. Different approaches
to simulating numerically the moving interface are discussed and compared.
To check whether the numerical procedure for simulating the reaction front
propagation can adequately reveal stable and unstable behavior of the interface,
several problems are examined numerically and compared with the analytical
solutions and stability predictions from the previous chapter. Also, the remesh-
ing procedure is cross-validated with the other numerical method for simulating
the moving interface, namely the CutFEM. Chapter 4 contains theoretical and
experimental studies of the IMC growth at the interface between copper pads
and tin-based solder alloys. The first part of Chapter 4 reports experimental
findings of IMC growth in different microchips during a high-temperature stor-
age test. The growth kinetics is modeled employing a continuum model based
on the chemical affinity tensor concept. By combining experiment, theory, and
a comparison of experimental data and theoretical predictions, the values of
the diffusion coefficient and an estimate for the chemical reaction constants are
obtained. Each chapter is finalized by the conclusions section, and the final
Chapter 5 gives the general conclusions and the outlook.
Chapter 1. Introduction
Numerical and analytical studies of the chemical reaction front kinetics in solids 9
2 Chemical reaction front kinetics and
stability
2.1 Chemical affinity tensor and general problem
statement
Consider a solid body undergoing a chemical transformation caused by a
chemical reaction localized at a sharp interface. In general, the equation of a
chemical reaction can be written as follows
πβπ΅β+π*π΅*βπ+π΅+,(2.1)
where
πβ
,
π*
and
π+
are the stoichiometric coefficients,
π΅β
,
π΅*
and
π΅+
are
the species involved in the reaction: initial material, diffusing reactant, and the
reaction product, respectively. The reaction between
π΅β
and
π΅*
is localized
at the reaction front Ξand is supported by the diffusion of
π΅*
through
π΅+
as
schematically shown in Fig. 2.1. It is assumed that all delivered
π΅*
is consumed
by the reaction. The reaction is accompanied by the change of volume at the
reaction front. This deformation produces mechanical stresses at the front
which in turn may affect the front propagation. Here and further in the text
of the manuscript index β
β
β denotes the initial material, β
*
β states for the
diffusing constituent, and β+β denotes the transformed material, which is the
product of the chemical reaction.
The notion of the chemical affinity arises from fundamental results by Gibbs
and De Donde (see, e.g., Prigogine and Defay [1954]). It was shown that in the
case of a chemical reaction the factor conjugate to the reaction rate
π
in the
expression of the entropy production
π
multiplied by temperature
π
was equal
to the combination of the chemical potentials of the reaction constituents:
ππ =π΄π, π΄ =ββοΈππππππ,(2.2)
where
ππ
is the chemical potential of the
π
-th constituent (per unit mass),
ππ
is the molar mass of
π
-th constituent, the stoichiometric coefficient
ππ
contributes to the sum with a positive sign β+β if the
π
-th constituent is
produced in the reaction and with a negative sign β
β
β if the
π
-th constituent
10
Ξπ΅βπ΅+
π΅*
t0
u0
Ξ©+
Ξ©β
Fig. 2.1: A schematic representation of a localized chemical reaction in solids.
is consumed. For the reaction (2.1) the classical chemical affinity would be
π΄=πβπβπβ+π*π*π*βπ+π+π+,(2.3)
where
π+, πβ
and
π*
are the mass densities of the chemical potentials of the
constituents
π΅+, π΅β
and
π΅*
, respectively. The factor
π΄
was called the chemical
affinity of the reaction. Then a kinetic equation in a form of the dependence
of the reaction rate on the chemical affinity can be formulated. One of the
accepted dependencies (see Glansdorff and Prigogine [1971]) for the reaction
rate is the equation
π=π[οΈ1βexp (οΈβπ΄
π
π )οΈ]οΈ,(2.4)
where
π
is the so-called partial rate of a direct reaction, which is defined by
the concentrations of the reactants, π
is the universal gas constant.
The relationships
(2.2)
β
(2.4)
were formulated for reactions in systems like
gases or liquids for which chemical potentials and, therefore, chemical affinity
can be presented as scalar values. The observation that the phase equilibrium
and chemical reaction take place not just in a point but at an oriented area
element passing through the point, led to the idea of a tensorial nature of the
chemical potential and chemical affinity (see the discussion on the tensorial
nature of a chemical potential in Grinfeld [1991] and on the tensorial nature
of the chemical affinity in Freidin and Vilchevskaya [2020]). One can see here
an analogy to the concept of stress. A stress state is determined by a scalar
pressure acting in a point in the case of liquids and gases. However, in the case
of a solid, instead of pressure, a stress tensor determines a traction acting at
the oriented area element passing through a point.
It was shown that the reaction rate is determined by the normal component
π΄ππ
of the chemical affinity tensor
A
, see the derivations in Freidin [2013];
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 11
Freidin, Vilchevskaya, and Korolev [2014]; Freidin [2015], and the chemical
affinity tensor has the same mathematical form as the scalar affinity
(2.3)
but
with tensorial chemical potentials instead of scalar ones. A kinetic equation
similar to the equation
(2.4)
can be taken for the reaction rate
ππ
at the area
element with the normal n:
ππ=π*π[οΈ1βexp (οΈβπ΄ππ
π
π )οΈ]οΈ,(2.5)
where it is taken into account that in the case of the reaction between solid and
gaseous constituents
π
can be taken in a form
π
=
π*π
where
π
is the partial
molar concentration of the diffusing reactant. Mechanical stresses affect the
reaction rate via the normal component of the affinity tensor. The normal
component of the reaction front velocity
ππ
can be expressed in terms of the
reaction rate from the mass balance at the reaction front:
ππ=πβπβ
πβ
ππ.(2.6)
In order to handle further analytical calculations and to obtain the unknown
parameters of the model, linear elastic solid reactants are considered. In this
case the constitutive equations of π΅βand π΅+are
πβ=Cβ:πβ,
π+=C+:(π+βπtr),(2.7)
and the Helmholtz free energies of the solid reactants are represented by
πβ(π) = πβ+1
2πβ:Cβ:πβ,
π+(π) = π++1
2(π+βπtr):C+:(π+βπtr),
(2.8)
where
CΒ±
are stiffness tensors,
πΒ±
are temperature dependent chemical energies
of the reactants
π΅Β±
. An isotropic transformation strain tensor can be considered
πtr =πtrI,(2.9)
where principal strains
πtr
are the same at all points of the domain occupied
by
π΅+
, and differences related to the deviations of the concentration of the
diffusing reactant are neglected.
Two sources of the volume change are considered: deformation due to
the chemical transformation itself and deformation due to the diffusion of
the diffusing reactant (see, e.g., Freidin [2013]; MΓΌller, Vilchevskaya, and
Freidin [2015]; Freidin and Vilchevskaya [2020]). The total deformation can
Section 2.1. Chemical affinity tensor and general problem statement
12
be estimated as follows. A material
π΅β
is placed on one side of the reaction
front and the mix of the materials
π΅+
and
π΅*
is on the other side. Due to the
chemical reaction
(2.1)
, the elementary volume d
πβ
=
πβπβ/πβ
transforms
into the volume d
π+
=
π+π+/π+
where
πΒ±
are the reference mass densities
of solid reactants π΅Β±. The ratio of these volumes is
π½ch =π+π+/π+
πβπβ/πβ
.(2.10)
The diffusion may produce additional volume
ππ*π*/π*
inside material
π΅+
,
where
π*
is the reference mass density of the diffusing reactant
π΅*
. The
parameter
π
reflects the deformational interaction between the reactants
π΅*
and
π΅+
. In fact,
π΅β
βtransformsβ into the mix of
π΅+
and
π΅*
at the interface Ξ.
Then the ratio of stress-free volumes of materials coexisting across the reaction
front is given by
π½tr =π+π+/π++ππ*π*/π*
πβπβ/πβ
.(2.11)
The case
π
= 0 in
(2.11)
corresponds to a solid skeleton approach, according to
which the diffusion of
π΅*
is not accompanied by volume expansion of
π΅+
. The
case
π
= 1 corresponds to adding the volumes of
π΅+
and
π΅*
. The value of the
parameter
π
depends on the mechanism of the diffusion of
π΅*
and saturations
ability of
π΅+
with respect to diffusing
π΅*
. The diffusion of
π΅*
through
π΅+
may occur by a vacancy mechanism involving two counter fluxes: of
π΅*
and
of vacancies. Then at the reaction front atoms of
π΅*
take places of these
vacancies, and this may also lead to shrinkage of
π΅+
, which corresponds to
the case π < 0.
In the case of small deformations, the transformation strain can be calculated
as
πtr
= (
π½1/3
tr β
1). To keep strain compatibility conditions (i.e., continuity of a
body) across the reaction front in the presence of volume change
π½tr
, additional
strains appear and produce internal stresses. These internal stresses together
with external loading influence the reaction front velocity.
Further, it is assumed that the chemical potential of the diffusing reactant
is given by
π*π*=π*+π
π ln π
π*
,(2.12)
where
π*
is the chemical energy of the diffusing constituent
π΅*
and
π*
is a
reference concentration of π΅*.
In quasi static approach, neglecting the pressure terms acting on the diffusing
constituent (the solid skeleton approach) it can be shown (see, e.g., Freidin,
Vilchevskaya, and Korolev [2014]) that from
(2.7)
β
(2.12)
the normal component
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 13
of the chemical affinity tensor can be expressed as
π΄ππ =πβπβ
πβ
π+π*π
π ln π
π*
,(2.13)
where the contribution of mechanical and chemical energies is presented by
π=πΎ+1
2πβ:πββ1
2π+: (π+βπtr) + πΒ±:JπK.(2.14)
It should be noted that due to displacement and traction continuity conditions
πβ:JπK=π+:JπK. In the expression (2.14)
πΎ=πββπ++πβ
πβπβ
π*π*(2.15)
is the combination of the chemical energies of the reactants and square brackets
denote the jump of the value across the reaction front, JπK=π+βπβ.
In the works by Kubanov and Freidin [1988]; Freidin [1989]; Morozov and
Freidin [1998]; Freidin [2010] a phase transformation problem was considered
and expressions for the jump of the normal component of the Eshelby stress
tensor across the interface were derived. These results were used in later work by
Freidin, Vilchevskaya, and Korolev [2014]; Freidin [2015] where the expression
for the normal component of the chemical affinity tensor was derived through
strains at only one side of the chemical interface (e.g., here it is expressed by
the strains from the βββ side):
π΄ππ =πβπβ
πβ(οΈπΎβ1
2πtr :C+:πtr β1
2πβ: [[C]] : πβ+
πβ:C+:πtr +1
2qβ:K+(n) : qβ)οΈ+π*π
π ln π
π*
,
(2.16)
where
qΒ±= [[C]]:πΒ±βC+:πtr,
Kβ(n)={n Gβ(n)n}π ,Gβ(n)=(nΒ·CβΒ·n)β1.(2.17)
In
(2.17)
the upper and lower subscripts β+β and β
β
β in the relationships corre-
spond to each other,
G
(
n
)is the Fourier transform of the Green tensor (the in-
verse of the acoustic tensor), and
π
means symmetrization:
πΎππππ =π(ππΊπ)(πππ)
.
If the tensor [[
C
]]
β1
exists, then the mechanical contribution in the expression
for
π΄ππ
can be rewritten in a similar manner as it was shown in Kubanov and
Section 2.1. Chemical affinity tensor and general problem statement
14
Freidin [1988] for the case of phase transformations:
π΄ππ =πβπβ
πβ(οΈπΎ*β1
2qβ:(οΈ[[C]]β1βK+(n))οΈ:qβ)οΈ+π*π
π ln π
π*
,
where πΎ*=πΎβ1
2πtr : [[B]]β1:πtr,BΒ±=Cβ1
Β±.
(2.18)
One should note that the latter expression is suitable for analytical solutions and
for the following linear stability analysis. However, expressions
(2.13)
β
(2.14)
should be utilized in the numerical simulations.
The mechanical stresses
π
and the concentration
π
of the diffusing reac-
tant
π΅*
can be found from the solution of the system of equations which include:
(i) Mechanical equilibrium equations:
β Β· π=0(2.19)
with boundary and interface conditions:
u|Ξ©1=u0,πΒ·n|Ξ©2=t0,
JuK|Ξ=0,JπKΒ·n|Ξ=0,(2.20)
where
u0
and
t0
are the given displacement and traction vectors at the
corresponding outer boundaries Ξ©1and Ξ©2of the body;
(ii)The constitutive equations (2.7) of the solid reactants;
(iii) The diffusion equation defined over the domain occupied by the material
B+
Ξc = 0 ,(2.21)
with the boundary and interface conditions:
π·nΒ· βπβπ(π*βπ)=0 at Ξ©+,
π·nΒ· βπ+π*ππ= 0 at Ξ,(2.22)
where
π·
is the diffusivity,
π
is the mass transfer coefficient, or the dissolution
constant (Xin et al. [2003]; Lin et al. [2017]) at the outer boundary of the body,
n
is the normal directed outward the domain
π΅+
. For the sake of simplicity,
the steady-state diffusion is considered with stationary diffusion equation. One
should note, that if the characteristic time of the reaction is much smaller
than the relative time of the diffusion, then the chemical reaction is diffusion
controlled. For this case a non-stationary diffusion equation should be utilized.
The first boundary condition in
(2.22)
defines the flux of
π΅*
through the
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 15
external boundary in relation with the solubility of
π΅*
in
π΅+
and means that
the supply of
π΅*
through the boundary stops if the saturation
π*
is reached.
On the other hand, the Dirichlet boundary condition
π
=
π*
can be prescribed
at Ξ
+
, however, being less physically reasonable comparing to Robin boundary
condition (2.22)1.
The second boundary conditions in
(2.22)
follows from the mass balance
between the supplied and consumed
π΅*
at the reaction front and the assump-
tion that the velocity of diffusing particles is much greater than the velocity
of the reaction front. The reaction rate
ππ
in the boundary condition at the
chemical reaction front couples the elasticity and the diffusion problems. One
should note that starting from equation
(2.22)
, the coefficient
π*
is assumed
to be equal to one. This can be done by normalizing all of the stoichiometric
coefficients in
(2.1)
, i.e.,
πββπβ/π*, π+βπ+/π*
. This also leads to
avoiding a nonlinearity of the boundary condition at the moving interface
(2.22)2when it is rewritten in terms of concentrations (see, e.g., Freidin et al.
[2016]).
(iv) To close the system of equations one has to add the kinetic equation
(2.5)
for
ππ
which defines the reaction front velocity via
(2.6)
with the dependence
of π΄ππ on the stresses and concentration given by (2.13)β(2.15).
The reaction rate can be expressed in terms of the equilibrium concentration
πeq as shown, for instance, in Freidin, Vilchevskaya, and Korolev [2014]:
π=π*(πβπeq),where πeq =π*exp (οΈβπβπβ
πβ
π
π
π )οΈ,
where πis taken from (2.14), or can be rewritten as
π=πΎ(π)β1
2πtr :C+:πtr β1
2πβ:JCK:πβ+
πβ:C+:πtr +1
2qβ:K+(n) : qβ
=πΎ*β1
2qβ:(οΈ[[C]]β1βK+(n))οΈ:qβ,
(2.23)
similarly to (2.16) and (2.18).
It should be noted that if the concentration of the diffusing constituent at
the reaction front is equal to
πeq
then the driving force is equal to zero and the
thermodynamic equilibrium is reached.
Section 2.1. Chemical affinity tensor and general problem statement
16
2.2 On stability of a chemical interface
A special case when the driving force is equal to zero for all points of the interface
is considered as an equilibrium configuration. However, even if the interface
velocity is zero, an additional stability analysis of the βreaction blockingβ state
is required. The kinetic equation and the behavior of the perturbed solution on
the perturbed interface is analyzed in this section. To this end the procedure
described by Eremeev, Freidin, and Sharipova [2007] is followed, where the
authors analyzed the stability of the interface between two phases at the
mechanical equilibrium during a phase transformation. In the next sections it
will be shown that there are certain mathematical similarities between chemical
and phase transformation regarding the linear stability analysis. Therefore
phase transformation front kinetics and stability is studied as well as chemical
interface stability.
It is important to notice that within perturbation approach the configuration
is considered as stable if the equilibrium chemical reaction front with a normal
parallel to the tensile or compressive deformation is stable. In other words,
configurations, which tend to change their microstructure (e.g., optimal lami-
nates discussed in Antimonov, Cherkaev, and Freidin [2016]) are considered as
unstable.
The main idea behind perturbing the interface is that if the equilibrium
position is stable then the perturbed configuration will return to its unperturbed
state. As mentioned above, the equilibrium position of the interface is reached
when the thermodynamic force is equal to zero in all points of the interface.
Linear stability analysis considers small perturbations of the displacements
and position of the interface. Consequently, the displacement and the interface
position in the perturbed state are written as follows:
u=u0+w,R=R0+πn0,(2.24)
where
u0
and
R0
correspond to the equilibrium displacements and interface
position,
w
and
π
are the perturbations,
n0
is a normal to the unperturbed
interface (Fig. 2.2). Linearization of the Boundary Value Problem (BVP)
defined by
(2.19)
β
(2.20)
using the constitutive equation
(2.7)
and the additional
thermodynamic condition
π
= 0 leads to the following set of equations and
boundary conditions for
w
and
π
which were obtained in Eremeev, Freidin,
and Sharipova [2003]; Eremeev, Freidin, and Sharipova [2007]. In the domains
πΒ±the differential equation and the boundary conditions take the form:
β Β· πΒ±(w) = 0,πΒ±(w) = CΒ±:βw,
w|Ξ©1=0,nΒ·π+(w)|Ξ©2=0.(2.25)
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 17
Ξ0Ξ
R0
R
n0
n
Fig. 2.2: Unperturbed and perturbed interfaces Ξ0and Ξ.
Displacement and traction continuity conditions at the interface are
JwK=βπJnΒ· βu0K,nΒ·Jπ(w)K=βπΒ·Jπ0KβπnΒ·JnΒ· βπ0K.(2.26)
Analogously to the perturbed displacements, the perturbed solution for the
diffusion problem consists of the original solution c
0
and an additional small
perturbation term:
π=π0+π . (2.27)
A linearization of the boundary value problem defined by
(2.21)
β
(2.22)
provides
the following set of equations for the diffusion perturbation value
π
, which are
coupled with the displacement perturbation wby the chemical reaction rate
Ξπ = 0,
π·n0Β· βπ +π*πΏπ(w, π )=0 at the reaction front interface,
π·n0Β· βπ +πΌπ = 0 at the reactant supply surface.
(2.28)
In the following sections a detailed derivation of the boundary and interface
conditions is given for the perturbed boundary value problem written in
(2.25)
β
(2.28)
. For this purpose a variation of the normal vector, jump across the
interface conditions, the driving force and the chemical reaction rate are derived.
2.2.1 Variation of the normal
If
π1, π2
are the curvilinear coordinates at the interface, then unperturbed and
perturbed interfaces are presented by the vectorial functions
R
=
R0
Ξ
(
π1, π2
)
and
R
=
RΞ
(
π1, π2
), respectively. Then the basis vectors on the unperturbed
Section 2.2. On stability of a chemical interface
18
and perturbed interfaces are defined by
R0
πΌ=πR0
Ξ
πππΌ,RπΌ=R0
πΌ+π(πn)
πππΌ, πΌ = 1,2.(2.29)
Due to the perturbation, the normal to the interface becomes
n=n0+πΏn.(2.30)
Then
nΒ·RπΌ= 0 = (n0+πΏn)Β·R0
πΌ+ (n0+πΏn)Β·(οΈπn0
πππΌπ+n0ππ
πππΌ)οΈ.(2.31)
Neglecting second order terms and taking into account relationships
n0Β·πΏn
= 0
and n0Β·πn0
πππΌ= 0, the following expression can be derived
πΏnΒ·R0
πΌ=βππ
πππΌ.(2.32)
Since
πΏnΒ·R0
πΌ
=
πΏππΌ
is the covariant component of the vector
πΏn
in the dual
basis RπΌ0, the above equality takes the vectorial form
πΏn=βξ₯
βπ, where ξ₯
β=RπΌ0π
πππΌ(2.33)
stands for a 2D nabla-operator in π1-π2coordinate system.
2.2.2 Variation of a jump across the interface
To formulate displacement and traction continuity conditions across the per-
turbed interface, a general expressions of the variation of a jump due to
variation of the interface position and displacement field is derived.
Given displacement field u, consider a function
π(R|u) = {οΈπβ(R|u),Rβπ£β
π+(R|u),Rβπ£+
where the functions
πβ
(
R|u
)and
π+
(
R|u
)are smooth enough in domains
π£β
and
π£+
. Denote by
πΒ±
(
RΞ|u
)and
βπΒ±
(
RΞ|u
)the limit values of
π
and
βπ
at the corresponding sides of the interface. Referring to the Figure 2.2, consider
two configurations: (i) with the unperturbed interface
R0
Ξ
at displacement
field
u0
and (ii) with the perturbed interface
RΞ
=
R0
Ξ
+
πΏRΞ
at perturbed
displacement
u
=
u0
+
w
. Then the jumps of
π(RΞ,u)
across perturbed
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 19
interface, the jump of
π(οΈR0
Ξ,u0)οΈ
across unperturbed interface and the jump of
the variations
πΏπ’π(οΈR0
Ξ|u0)οΈ
of
π
due to displacement perturbations are related
as
Jπ(RΞ|u)K=Jπ(οΈR0
Ξ|u0)οΈK+JπΏπ’π(οΈR0
Ξ|u0)οΈK+πΏRΞΒ·Jβπ(οΈR0
Ξ|u0)οΈK.(2.34)
To prove
(2.34)
note that the functions
π+
(
R|u
)and
πβ
(
R|u
)are defined
only in the domains
π£+
and
π£β
, respectively, and the interface position is given
by
R0
Ξ
if the displacement is
u0
, and by
RΞ
if the displacement is
u
. Therefore
π+
is defined in points
R0
Ξ
and
RΞ
at displacement
u
and in points
R0
Ξ
at both
displacements
u0
and
u
. Then, neglecting second (and higher) order terms
with respect to perturbations one can obtain the following expressions
π+(RΞ|u) = π+(οΈR0
Ξ,u)οΈ+πΏRΞΒ· βπ+(οΈR0
Ξ|u)οΈ
π+(οΈR0
Ξ|u)οΈ=π+(οΈR0
Ξ|u0)οΈ+πΏπ’π+(οΈR0
Ξ|u0)οΈ(2.35)
and, thus, with the same accuracy
π+(RΞ|u) = π+(οΈR0
Ξ|u0)οΈ+πΏπ’π+(οΈR0
Ξ|u0)οΈ+πΏRΞΒ· βπ+(οΈR0
Ξ|u0)οΈ.(2.36)
Analogously, since function
πβ
is defined in points
RΞ
at displacement
u
and
at displacement
u0
and in points
RΞ
and
R0
Ξ
at displacement
u0
, one can write
πβ(RΞ|u) = πβ(οΈR0
Ξ|u0)οΈ+πΏπ’πβ(οΈR0
Ξ|u0)οΈ+πΏRΞΒ· βπβ(οΈR0
Ξ|u0)οΈ.(2.37)
Subtracting (2.37) from (2.36) gives the formula (2.34).
If quantity
π
denotes the displacement vector field
u
then
πΏπ’π
denotes
the displacement perturbation
w
. Since, due to displacement continuity,
Ju0(R0
Ξ)K=Ju(RΞ)K= 0, from (2.34) it follows that
JwK+πΏRΞΒ·Jβu0K= 0.(2.38)
From the displacement continuity it follows that
Jβu0K
=
n0a
, where
a
=
n0Β·Jβu0K
is the amplitude of jump. Then only normal component
π
of
πΏRΞ
remains in the displacement continuity condition that becomes so-called second
order compatibility condition (2.26)1:
JwK+πn0Β·Jβu0K= 0.(2.39)
Further, without loss of generality, we assume that πΏRΞ=πn0.
Section 2.2. On stability of a chemical interface
20
For the jump of stresses (2.34) can be rewritten as follows
Jπ(RΞ,u)K=Jπ(οΈR0
Ξ,u0)οΈK+Jπ(οΈR0
Ξ,w)οΈK+πn0Β·Jβπ(οΈR0
Ξ,u0)οΈK(2.40)
where
π(οΈR,u0)οΈ
and
πΒ±
(
R,u
)are given by
(2.7)
with strain tensors
π0
=
βu0
and
π
=
βu
, respectively,
πΒ±(w) = CΒ±:βw
. Then at the perturbed interface
nΒ·Jπ(RΞ,u)K= (n0+πΏn)Β·(οΈJπ(οΈR0
Ξ,u0)οΈK+Jπ(οΈR0
Ξ,w)οΈK+
πn0Β·Jβπ(οΈR0
Ξ,u0)οΈK)οΈ=0.
(2.41)
Taking into account that
n0Β·Jπ(οΈR0
Ξ,u0)οΈK
= 0, using the formula
(2.33)
for
πΏn
and neglecting second order terms we come to the final expression
(2.26)2
for
the linearized traction continuity condition at the perturbed interface
n0Β·Jπ(w)K=ξ₯
βπΒ·Jπ0Kβπn0Β·Jβπ0KΒ·n0,(2.42)
where π0
Β±=πΒ±(οΈR0
Ξ,u0)οΈand πΒ±(w) = πΒ±(οΈR0
Ξ,w)οΈ.
2.2.3 Linearized kinetic equation for perturbed phase interfaces
By (2.18), the driving forces at the perturbed and unperturbed interface are
π|Ξ=πΎ*β1
2(οΈq0
β+πΏqβ)οΈ:(οΈJCKβ1βK+(n0+πΏn))οΈ:(οΈq0
β+πΏqβ)οΈ,
π|Ξ0=πΎ*β1
2q0
β:(οΈJCKβ1βK+(n0))οΈ:q0
β,
(2.43)
where q0
β= [[C]]:π0
ββC+:πch,πΏqβis calculated by Eq. (2.37):
πΏqβ=πΏπ’q0
β+πn0Β· βq0
β=JCK:βwβ+πn0Β·(οΈβπ0
β:JCK)οΈ.(2.44)
Then, with
K+
(
n0
+
πΏn
) =
K
(
n0
) +
πΏK
(
n0
), the expression of the driving
force variation takes the form
πΏπ =βq0
β:(οΈJCKβ1βK+(n0))οΈ:πΏqβ+1
2q0
β:πΏK(n0) : q0
β.(2.45)
In the case when the normal
n0
is an eigenvector of
q0
β
and material β+β
is isotropic, the last term in
(2.45)
is equal to zero. Indeed, the fourth rank
tensor
K+(n)
in the case of an isotropic material can be presented with the
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 21
use of in special tensorial basis (see Kunin [1983]) as
K+(n) = 1
π+(οΈE5βΞΊ+E6)οΈ,(2.46)
where
E5= (neπeπn)π =1
4(neπeπn+eπneπn+neπneπ+eπnneπ),
E6=nnnn,
ΞΊ+=π++π+
π++ 2π+
=1
2(1 βπ+),
upper script
π
denotes here the symmetrization of a fourth rank tensor with
respect to permutation of indices within the first and second pairs, vectors
eπ
(
π
= 1
,
2
,
3) form an orthonormal basis, and
eπeπ
is a second rank unit
tensor. Note that
E5
is also symmetric with respect to the permutation of the
pairs of indices. The variation of E5equals to
πΏE5=1
4((πΏn)eπeπn+neπeππΏn+eπ(πΏn)eπn+eπneππΏn+
+ (πΏn)eπneπ+neπ(πΏn)eπ+eπ(πΏn)neπ+eπn(πΏn)eπ)
Then, since nΒ·πΏn= 0, it is clear that if nis an eigenvector of q,
q:πΏE5:q= 0 (2.47)
and, analogously, q:πΏE6:q= 0, and finally
q:πΏK:q= 0.
By
(2.44)
and
(2.45)
, with this additional assumption about the eigenvector of
q0, the expression for the variation of the driving force takes the form
πΏπ =βq0
β:(οΈJCKβ1βK+(n0))οΈ:JCK:πβ(wβ)β
πq0
β:(οΈJCKβ1βK+(n0))οΈ:(οΈn0Β·(οΈβπ0
β:JCK)οΈ)οΈ.
(2.48)
The variation can also be expressed in terms of strains on the β+β side of the
interface:
πΏπ =βq0
+:(οΈJCKβ1+Kβ(n0))οΈ:JCK:π+(w+)β
πq0
+:(οΈJCKβ1+Kβ(n0))οΈ:(οΈn0Β·(οΈβπ0
+:JCK)οΈ)οΈ.
(2.49)
Section 2.2. On stability of a chemical interface
22
Note that the last terms in Eqs.
(2.48)
and
(2.49)
are zero if the normal
derivatives
n0Β· βπ0
β
of the strains
π0
β
are zero. This is, for example, the case
of piece-wise homogeneous two phase deformations with plane interfaces, or
in the case of spherically or axially symmetric two-phase deformations in a
solid cylinder or sphere when the strain is uniform in the inner domain. In the
last case one can take representation
(2.48)
or
(2.49)
, depending on what the
phase, βββ or β+β, occupies the inner domain.
The normal components
ππ
and
π0
π
of the velocities
V
and
V0
of the
perturbed and unperturbed interfaces with an accuracy of the second order of
smallness are related as
ππ=πRΞ
ππ‘ Β·n=π(οΈR0
Ξ+πn0)οΈ
ππ‘ Β·(n0+πΏn) = π0
π+V0Β·πΏn+ππ
ππ‘ .(2.50)
Since only normal component of the interface velocity is essential, one may
accept V0=π0
πn0. Then
πΏππ=ππβπ0
π=dπ
dπ‘.(2.51)
The kinetic equation for the phase transformation front is
πph
π
=
π
phπ
. Then
from
(2.48)
,
(2.49)
and
(2.51)
it follows that perturbations evolve according to
the kinetic equation that can be presented in the following forms:
1
π
ph
dπ
dπ‘=βq0
β:(οΈJCKβ1βK+(n0))οΈ:JCK:πβ(wβ)β
πq0
β:(οΈJCKβ1βK+(n0))οΈ:(οΈn0Β·(οΈβπ0
β:JCK)οΈ)οΈ
=q0
+:(οΈJCKβ1+Kβ(n0))οΈ:JCK:π+(w+)
βπq0
+:(οΈJCKβ1+Kβ(n0))οΈ:(οΈn0Β·(οΈβπ0
+:JCK)οΈ)οΈ.
(2.52)
2.2.4 Linearized kinetic equation for perturbed chemical interfaces
The variation of thr reaction rate
πΏπ
in
(2.28)
depends on the displacements
perturbation vector
w
and the diffusion perturbation
π
. Keeping in mind that
π΄ππ β
0and
π*
=
π0
=
πeq
near equilibrium (for the details we refer to Freidin,
Vilchevskaya, and Korolev [2014]), the reaction rate can be expressed by using
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 23
Eq. (2.5):
π(w, π ) = π*π(οΈ1βexp (οΈβπ΄ππ(w, π )
π
π )οΈ)οΈ=π*(π0+π )π΄ππ(w, π )
π
π
=π*π0
π
π (οΈπβπβ
πβ
πΏπ +π
π ln (οΈπ0+π
π*)οΈ)οΈ+Β· Β· Β·
=π*π0
π
π (οΈπβπβ
πβ
πΏπ +π
π π
π*)οΈ+Β· Β· Β· ,
(2.53)
where
πΏπ
is a variation of the mechanical part of the chemical affinity tensor,
which has the same form as the variation of the driving force in the case of
phase transformations
(2.48)
. The variation of the diffusion related logarithmic
term in the expression for π΄ππ can be obtained as follows
πΏ(οΈln π(R)
π*)οΈ= ln π
π*ξξξξR0
+πΏπ
πξξξξR0
+πn0Β· β ln π
π*ξξξξR0
+Β· Β· Β· .(2.54)
At the thermodynamical equilibrium the concentration of the diffusing com-
ponent in the reactant material is constant (as can be seen from the solution
in Freidin, Vilchevskaya, and Korolev [2014]) and equal to
π*
. That means
that the first and the last terms in Eq (2.54) are equal to zero. According to
the notation in Eq (2.27), the variation of the concentration at the interface is
equal to π . Hence,
πΏ(οΈln π(r)
π*)οΈ=π
π*
.(2.55)
Finally, the linearized kinetic equation is obtained by substituting the linearized
reaction rate (2.53) into the kinetic equation (2.6):
1
π
ch
dπ
dπ‘=πβπβ
πβ[οΈβq0
β:(οΈJCKβ1βK+(n0))οΈ:JCK:πβ(wβ)β
πq0
β:(οΈJCKβ1βK+(n0))οΈ:(οΈn0Β·(οΈβπ0
β:JCK)οΈ)οΈ]οΈ+π
π π
π*
(2.56)
where π
ch is a positive constant.
2.3 Stability of a planar chemical interface
As a first example, a stability of the planar interface is analyzed using the
aforementioned procedure. A chemical reaction in the layer of thickness
π»
is
considered, as illustrated in Fig. 2.3. Lower boundary is fixed and displacement
π’0
in
π¦
-direction is applied at the upper boundary. This external load is further
considered as a parameter for controlling the resulting equilibrium position.
Section 2.3. Stability of a planar chemical interface
24
Diffusing constituent is supplied through the lower boundary and the reaction
takes place at
π¦
=
β
. This gives the following boundary condition for the
mechanical equilibrium equation, (2.19):
u|π¦=0 =0,u|π¦=π»=π’0eπ¦,
π’π₯= 0, ππ₯π¦ = 0,at π₯= 0 and π₯=πΏ. (2.57)
and for the diffusion equation, (2.21):
π·nΒ· βπ+πΌ(πβπ*)=0,at π¦= 0,
nΒ· βπ= 0,at π₯= 0 and π₯=πΏ. (2.58)
Linear elastic solid reactants are considered with the constitutive equations
(2.7) and elasticity tensors taken in a form
CΒ±=πΒ±II + 2πΒ±4I,(2.59)
where
πΒ±
and
πΒ±
are LamΓ© parameters and
4I
is the forth rank unit tensor.
For the sake of simplicity of further derivations, the chemical transformation
strains assumed to be planar, namely
πtr
=
πtr
(
eπ₯eπ₯
+
eπ¦eπ¦
). One should note
that the chosen boundary conditions (for both mechanical equilibrium and
diffusion problem) on left and right boundaries of the layer may be considered
as a way to model an infinite layer.
π₯
π¦
β
π»
πΏ0
-
β
+
β
Fig. 2.3: A schematic representation of the planar chemical reaction front.
Assuming that the parameters of the problem are such that there exists an
equilibrium configuration of the reaction front inside the layer, the mathematical
problem is to determine whether this configuration is physically stable or
unstable. Furthermore, only planar equilibrium configurations of the front are
considered, i.e., all points of the front belong to the line
π¦
=
β
, where
β
is a
constant. At this state, the solution of the balance equation and the diffusion
equation are displacements
u0
and concentration
π0
=
π*
, the stresses emerging
in the body are denoted as
π0
. The closed form expressions for these quantities
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 25
are shown in Appendix A.1. From this point, the superscript β0β refers to the
variable taken from the solution of the unperturbed BVP.
According to the utilized approach, as described in Section 2.2, the equilib-
rium configuration of the front is perturbed,
R=R0+π(π₯, π‘)n0,(2.60)
where
R
and
R0
are positions of the points of the front at the perturbed and
the unperturbed states, respectively,
n0
=
eπ¦
at the equilibrium,
π(π₯, π‘)
is the
amplitude of the perturbation. The perturbation of the front necessary leads
to the perturbation of the solutions of the PDEs:
u(π₯, π¦, π‘) = u0(π¦) + w(π₯, π¦, π‘),
w(π₯, π¦, π‘) = π’(π₯, π¦, π‘)eπ₯+π£(π₯, π¦, π‘)eπ¦,
π(π₯, π¦, π‘) = π0(π¦) + π (π₯, π¦, π‘),
(2.61)
where
w
and
π
are the perturbation functions for the displacements and
concentration respectively. From this point, the superscript β0β refers to the
variable taken from the solution of the unperturbed BVP. Linearization of the
balance equation gives
β Β· πΒ±(w) = 0,πΒ±(w) = CΒ±:βw.(2.62)
The boundary conditions from (2.57) are:
w=0,at π¦= 0 and π¦=π»,
wΒ·eπ₯= 0,πΒ±(w) : eπ₯eπ¦= 0,at π₯= 0 and π₯=πΏ, (2.63)
and interface conditions at π¦=βare defined in (2.26).
The linearized diffusion equation reads:
Ξπ = 0,(2.64)
and the linearized boundary conditions (2.58) are
π·nΒ· βπ +πΌπ = 0,at π¦= 0,
π·n0Β· βπ +π(w, π )=0,at π¦=β,
nΒ· βπ = 0,at π₯= 0 and π₯=πΏ.
(2.65)
Variations of energy
πΏπ
(
w
)and the reaction rate
πΏπ (w, π )
have form
(2.48)
and (2.53), respectively.
Solution of the perturbed system of equations
(2.62)
and
(2.64)
defines the
dependency of displacements
u
and concentration
π
(and therefore the reaction
Section 2.3. Stability of a planar chemical interface
26
front velocity) on the amplitude of the perturbation
π
. For the case of the
planar interface, these dependencies can be obtained analytically in a series
form. To satisfy the boundary conditions at
π₯
= 0 and
π₯
=
πΏ
, the solution is
found as series
π’(π₯, π¦, π‘) =
β
βοΈ
π=1
ππ(π¦, π‘) sin (πππ₯), π£ (π₯, π¦, π‘) =
β
βοΈ
π=1
ππ(π¦, π‘) cos (πππ₯),
π(π₯, π‘) =
β
βοΈ
π=1
ππ(π‘) cos (πππ₯), π (π₯, π¦, π‘) =
β
βοΈ
π=1
ππ(π¦, π‘) cos (πππ₯), ππ=ππ
πΏ.
One can show that functions
ππ(π¦, π‘)
and
ππ(π¦, π‘)
which satisfy equation
(2.62)
are
ππ(π¦, π‘) = π΄π(π‘) exp (πππ¦) + π΅π(π‘)π¦exp (πππ¦) +
πΆπ(π‘) exp (βπππ¦) + π·π(π‘)π¦exp (βπππ¦),
ππ(π¦, π‘) = (οΈβπ΄π(π‘) + π΅π(π‘)π+π+ 2
π+π
1
ππ)οΈexp (πππ¦)βπ΅π(π‘)π¦exp (πππ¦) +
(οΈπΆπ(π‘) + π·π(π‘)π+π+ 2
π+π
1
ππ)οΈexp (βπππ¦) + π·π(π‘)π¦exp (βπππ¦).
Since the material properties are different for domains
πΊ+
and
πΊβ
, there are
two functions
ππ
. The first is defined in
πΊ+
with LamΓ© parameters
π+
and
π+
,
the second is defined in
πΊβ
with LamΓ© parameters
πβ
and
πβ
. Analogously,
there are two functions
ππ
. Functions
π΄π
,
π΅π
,
πΆπ
,
π·π
are then found from
the boundary conditions
(2.63)
. Time dependency in these functions comes
from dependency of the boundary conditions on
π(π₯, π‘)
. It is easy to see from
the boundary conditions that each function
π΄π
,
π΅π
,
πΆπ
,
π·π
is proportional to
ππ(π‘)
, e.g.,
π΄π(π‘)
=
ππππ(π‘)
, where
ππ
is some constant. This allows expressing
variation πΏπ as a function of ππ:
πΏπ (π‘) =
β
βοΈ
π=1
πΏph
πππ(π‘) cos (πππ₯),(2.66)
where
πΏph
π
is a constant. Here the superscript β
ph
β refers to the phase transfor-
mation, since the expression
(2.66)
is to be used in the kinetic equation
(2.52)
if the phase transformation interface stability is considered.
Further, the perturbed diffusion equation
(2.64)
has to be solved. Function
ππ(π¦, π‘)can be found in a form
ππ(π¦, π‘) = πΈπ(π‘) exp (πππ¦) + πΉπ(π‘) exp (βπππ¦),(2.67)
where functions
πΈπ
,
πΉπ
are found from the boundary conditions. Furthermore,
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 27
as seen from the boundary conditions, each function
πΈπ
,
πΉπ
is proportional
to
ππ(π‘)
, e.g.,
πΉπ(π‘)
=
πππΏph
πππ(π‘)
, where
ππ
is a constant. Thus, from the
boundary conditions, it is easy to show that
πΈπ(π‘) = πΏph
πππ(π‘)βπ
ππ
, πΉπ(π‘) = πΏph
πππ(π‘)βπ
ππ
π·ππβπΌ
π·ππ+πΌ,(2.68)
where the following notation is introduced:
ππ=π·ππexp (ππβ)βπ·ππ
π·ππβπΌ
π·ππ+πΌexp (βππβ) +
π*exp (ππβ) + π*
π·ππβπΌ
π·ππ+πΌexp (βππβ),
π
=π*π*
1
π
π
πβπβ
πβ
.
When all the constants are found, and the corresponding perturbed functions
w
and
π
are obtained, one can substitute them into the kinetic equation
(2.56)
.
Then the following equation for the evolution of the perturbation
ππ
is obtained:
1
π
πβ
πβπβ
dππ
dπ‘=πΏch
πππ,
where
πΏch
π=πΏph
πβ
β
β
β1β
π*exp (ππβ) + π*
π·ππβπΌ
π·ππ+πΌexp (βππβ)
ππβ
β
β
β .(2.69)
Sign of
πΏch
π
defines the behavior of the solution: when the sign is negative for
each
π
, the perturbation decays in time, and the reaction front moves back to
its equilibrium position. If at least one
πΏπ
is positive, than the perturbation
with this wave number grows exponentially, which leads to the instability of
the interface.
The specific form of the obtained expression for
πΏch
π
is very useful. One
can show that the expression in parentheses in
(2.69)
is always positive. This
means that sings of
πΏch
π
and
πΏph
π
always coincide. If a stress-induced phase
transformation problem is considered, for which the velocity is defined by
equation (2.52), similar analysis of a perturbed boundary gives
1
πph
dππ
dπ‘=πΏph
πππ.
Since
πΏch
π
and
πΏph
π
have identical signs, the physical stability of a phase boundary
Section 2.3. Stability of a planar chemical interface
28
(π)
πΏch
π
βeq
(π)
πΏch
π
βeq
Fig. 2.4: πΏπ
-lines for the stability analysis of the planar chemical reaction front for
stable (π)and unstable (π)sets of elastic parameters.
in the phase-transformation problem necessarily results in the physical stability
of a chemical reaction front in the chemo-mechanical problem. Therefore, the
diffusing reactant does not influence stability of the front. The diffusion process
cannot stabilize the physically unstable front and also otherwise.
In Fig. 2.4 different possible scenarios are shown. The plots (a) and (b)
correspond to stable and unstable sets of material parameters for the planar
interface, respectively. For the stable configuration
πΏch
π
is negative for each
π
and all possible positions of the equilibrium position (which, as mentioned
earlier, is controlled by the external load
π’0
). For the unstable configuration
the solution shows positive values of the πΏch
πβs .
It should be noted that negative values of the
πΏπ
βs is only a necessary but not a
sufficient stability condition. For the problem of phase transformation, which is
mathematically pretty similar to the problem of chemical transformation, there
is another necessary stability criteria which is based on the analysis of so called
phase transition zones (PTZ). It was proved in Grabovsky and Truskinovsky
[2011]; Grabovsky and Truskinovsky [2013] that if strains at the equilibrium
interface belong to the boundaries of the phase transition zones, then this
configuration is stable. This approach cannot be explicitly extrapolated to the
case of the chemical reaction, therefore PTZ and perturbation approaches are
later compared for the phase transformation problem.
Another weakness of the perturbation method is as follows. Assuming
the aforementioned configuration with the planar chemical interface with the
given material properties, one can obtain values for
πΏch
π
for each positive
integer
π
. However, with growing
π
, system of equations for functions
π΄π
,
π΅π
,
πΆπ
,
π·π
becomes ill-conditioned. It means, that for higher wave numbers
one cannot guarantee an accurate solution. However, for several observed
cases, the absolute value of πΏch
πincreases with πfor both stable and unstable
configurations, as shown in Fig. 2.4.
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 29
One should note, that the material parameters for the following problems
are chosen to be dimensionless because of the following reasons. The purpose
of Sections 2.3β2.5 is the qualitative analysis of the stability behavior, that is
validated by the numerical results in the following Chapter 3. When considering
stability analysis, the material parameters are varied to consider different
physical scenarios, therefore, there is no reference to the real engineering
problem. Also, during the development of the numerical procedures, it is a
common practice to use dimensionless units. The units are taken into account
in the last Chapter 4, where the real engineering application is considered.
Having all this in mind, the following procedure was utilized in order to
obtain a stability regions in the material parametric space. It was shown,
that diffusion parameters of the model do not influence the stability behavior,
so they are kept fixed. In addition the Poissonβs ratios are fixed so that
π+
=
πβ
= 0
.
25, the elasticity modulus for one of the materials is selected to be
πΈβ
= 60, external load was set to
π’0
= 0
.
02. With this, the parametric space
for the following study is defined by ratio
πΈ+/πΈβ
, the value of the chemical
transformation
πtr
and the energy parameter
πΎ
. In this work for the sake of
convenience, the chemical energy parameter
πΎ
was adjusted in order to get
the equilibrium position at the middle (
βππ
= 0
.
5) of the layer. From physical
point of view, variation of πΎmay correspond to the variation of temperature.
Varying the selected parameters, one can obtain the areas, where the values
for πΏ1are positive or negative, Fig. 2.5.
πtr
πΈ+
πΈβ
(π)
πtr
πΈ+
πΈβ
(π)
Fig. 2.5:
The contour plots of: (
π
)β
πΎ
for given external load in material parameter
space; (
π
)β
πΏ1
for given external load in material parameter space at corresponding
πΎ
from figure (a). Black solid lines in (b) represent the πΏ1= 0 level.
The stability regions in Fig. 2.5 prescribe only a necessary condition, meaning
that if the region says βunstableβ (
πΏ1>
0) then the system is physically unstable.
However, if the region states βstableβ (
πΏ1<
0), an additional check for
πβ₯
2is
required. Using these plots the stable (with the additional check) and unstable
configurations are selected for further study. The sets of elasticity parameters
at which the equilibrium planar chemical reaction front with a normal parallel
Section 2.3. Stability of a planar chemical interface
30
π+πβπ+πβπtr
Stable 40 10 34 26 0.005
Unstable 66 34 50 66
Tab. 2.1:
Stable and unstable sets of material properties and parameters used in
analytical study and numerical simulation for the planar interface problem.
to the tensile or compressive deformation is stable or unstable are listed in
the Tab. 2.1. Later in the text set of parameters for which the equilibrium
position of the interface is stable or unstable is reffered to as βstable set of
parametersβ or βunstable set of parametersβ, respectively.
It has been shown that chemical and diffusion parameters do not affect the
stability of the planar interface in the equilibrium position. Therefore only one
set of the following chemical and diffusion parameters are used and they are
listed in Tab. 2.2.
Parameter πβπβ
πβ
π
π π*π· π*πΌ
Value 43.2 2434.8 0.1 0.1 0.01 0.2
Tab. 2.2:
Diffusion parameters used in analytical study and numerical simulation for
the planar interface problem.
These parameters are used in the following Chapter 3 as a reference values
to validate numerical procedures.
2.4 Stability of a cylindrical chemical interface
Consider a hollow cylinder undergoing a a chemical reaction, Fig. 2.6. Due to
symmetry only a quarter of the cylindrical cross-section will be shown in all
further figures.
The problem is axially symmetric, which means that
u
=
π’
(
π
)
eπ
. As it was
done in the case of planar interface, the transformation strains are assumed to
be planar:
πtr =πtr(eπeπ+eπeπ),(2.70)
where
eπ
and
eπ
are unit vectors of a cylindrical coordinate system. The
boundary conditions refer to the quantities shown in Fig. 2.6, where
π
is an
interface radius,
π
and
π
are inner and outer radii of the cylinder, respectively,
can be written as
ππ(π)=0, π’(π) = π’0,
Jπ’(π)K= 0,Jππ(π)K= 0.
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 31
Ξ
-
β+
β
π’0
ππ
π
Fig. 2.6: A hollow cylinder undergoing a phase transformation.
Diffusing reactant for the chemical reaction is supplied through the outer
surface of the cylinder. As it was done in previous sections a stationary
diffusion is considered. The corresponding diffusion equation in the case of
axial symmetry becomes:
d
dπ(οΈ1
π
dπ
dπ)οΈ= 0, π β[π, π].(2.71)
Boundary and interface conditions according to the notation of Figure 2.6 can
be written as follows (as described in Section 2.1):
π·dπ
dπξξξξπ=π
βπΌ(π*βπ(π)) = 0, π· dπ
dπξξξξπ=π
βπ*(π(π)βπeq(π)) = 0.(2.72)
The solution to this particular problem can be found analytically as it
was shown in details, e.g., in Vilchevskaya and Freidin [2007] for the phase
transformation front kinetics problem. This solution and its extension for the
chemical reaction front propagation are given in the Appendix A.2.
If a cylinder undergoing a phase transformation,the phase transition front
velocity is given by the following expression (see Freidin [2007])
πph
π=βπΛ = π
phπ(π), π
ph >0,(2.73)
where
π
is the same as in the expression of the chemical affinity tensor
π΄ππ
,
(2.14), but without the chemical energy term in the expression for the energy
parameter
πΎ
,
(2.15)
. Parameter
π
ph
is a positive constant. One should note
that in current problem formulation
n
=
βeπ
. For a given load
π’0
and
choice of material parameters the dependence of the driving force
π
on the
interface radius can be plotted. This driving force may have a zero value for
Section 2.4. Stability of a cylindrical chemical interface
32
a certain radius as shown in Fig. 2.7. For this radius the normal component
of the velocity of the interface is zero, so this radius is then referred to as an
equilibrium position.
ππ
Fig. 2.7: Dependence of the driving force on the interface radius.
Note that if the initial radius of the interface is greater than the equilibrium
one, the sign of the driving force
π
is positive. Then according to
(2.73)
the
normal component of the velocity is positive. Keeping in mind that n=βeπ,
the phase transformation front is forced to move toward the equilibrium position.
The same happens if the initial radius of the interface is less then equilibrium
radius. This means that the solution is stable at the equilibrium radius with
respect to the radial perturbations.
The external load
π’0
can be considered as a parameter for controlling the
resulting equilibrium radius. Examples of the dependencies of the position of
the equilibrium radius on the external load are shown in Fig. 2.8 (a) and (b)
for a solid and for a hollow cylinder, respectively. Note that an equilibrium
radius may not exist within the range between internal and external radii of
the cylinder for some values of
π’0
. In the case of a hollow sphere or a hollow
cylinder two equilibrium radii may correspond to one external load. A more
detailed discussion of this phenomenon for the spherical problem can be found
in Eremeev, Freidin, and Sharipova [2007]. Moreover, the case of an external
load at which there is only one equilibrium radius and the interface is stable
with respect to radial perturbations is considered in that reference.
As mentioned above, the equilibrium position of the interface is reached
when the thermodynamic force is equal to zero in all points of the interface.
Linear stability analysis considers small perturbations of the displacements
and position of the interface. Consequently, similarly to the planar interface
problem in Section 2.3, the displacement, concentration and the interface
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 33
(π) (π)
u0u0
πeq πeq
Fig. 2.8:
Dependence of thermodynamically equilibrium interface radius on the exterξΏ
nal load for (
π
)β solid and (
π
)β hollow cylinders. Dashed line in (
π
)corresponds to
the inner radius of the hollow cylinder.
position in a cylindrical coordinate can be written as
u=u0(π) + w(π, π),R=R0+π(π)n0,
w(π, π) = π’(π, π)eπ+π£(π, π)eπ,
π=π0(π) + π (π, π).
(2.74)
For the case of the hollow cylinder equilibrium equation
(2.25)
can be
rewritten as shown in Vilchevskaya and Freidin [2007]:
(π+ 2π)ππ
ππ β2π
π
ππ
ππ = 0,(π+ 2π)1
π
ππ
ππ β2πππ
ππ = 0,(2.75)
where
π=ππ’
ππ +1
π(οΈπ’+ππ£
ππ)οΈ, π =ππ£
ππ +1
π(οΈπ£βππ’
ππ)οΈ.
Then the boundary conditions on the free inner surface (π=π) are
ππ + 2πππ’
ππ = 0,ππ£
ππ = 0.
On the outer surface (
π
=
π
), where Dirichlet boundary conditions are applied,
the perturbation function has to vanish. Therefore
π’= 0, π£ = 0.
The interface conditions at
π
=
π
from
(2.26)
can be written for the new
Section 2.4. Stability of a cylindrical chemical interface
34
variables π’and π£as
Jπ’K=βπ[οΈ[οΈdπ’0
dπ]οΈ]οΈ,Jπ£K= 0,
[οΈ[οΈππ + 2πππ’
ππ ]οΈ]οΈ=βπ[οΈ[οΈdπ0
π
dπ]οΈ]οΈ,[οΈ[οΈ2πππ£
ππ ]οΈ]οΈ=dπ
dπJπ0
πK.
(2.76)
Linearized diffusion equation
(2.28)
for the perturbed concentration in cylin-
drical coordinates takes the form
Ξπ (π, π) = π2π
ππ2+1
π
ππ
ππ +1
π2
π2π
ππ2= 0 (2.77)
with boundary conditions
βπ·ππ
ππ +π*(οΈπ +π*
πβπβ
πβ
πΏπ
π
π )οΈ= 0 at the reaction front,
π·ππ
ππ +πΌπ = 0 at the outer surface.
(2.78)
Together with Eq. (2.75) and its boundary conditions, which are related to
the perturbed mechanical problem, equations
(2.77)
and
(2.78)
form a set of
equations for
π’
,
π£
and
π
. The aforementioned system of equations with the
boundary and interface conditions can be uniquely solved with respect to the
unknown functions. To do so, one has to follow the method proposed in Ere-
meev, Freidin, and Sharipova [2007] and the solution proposed by Vilchevskaya
and Freidin [2007] for the phase transformation problem. One can seek the
solution in a series form:
π’(π, π) =
β
βοΈ
π=2
ππ(π) cos(ππ), π£(π, π) =
β
βοΈ
π=2
ππ(π) sin(ππ),
π(π) =
β
βοΈ
π=2
ππcos(ππ).
(2.79)
Then it is natural to seek the solution for π also as a series:
π (π, π) =
β
βοΈ
π=2
ππ(π) cos(ππ).(2.80)
Substitution of these series into the BVP provides a set of equations for the
amplitude functions
ππ
(
π
),
ππ
(
π
)and
ππ
(
π
)for each integer number
πβ₯
2,
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 35
which can be found uniquely. This means, that they can be expressed as
ππ(π) = π΄π
π(π)ππ, ππ(π) = π΄π
π(π)ππ, ππ(π) = π΄π
π(π)ππ(2.81)
where
π΄π
π, π΄π
π
and
π΄π
π
are certain linear operators. After the solution of these
equations is obtained and substituted into the kinetic equation
(2.56)
, the
latter takes the following form:
1
π
ch
dππ
dπ‘=πΏch
πππ,(2.82)
where
πΏch
π
is a linear integro-differential operator. Similarly to the planar
interface solution given in Section 2.3, the solution of this equation is an
exponential one, and its power is defined by
πΏch
π
. The character of the solution
is then governed by the sign of
πΏch
π
. If the sign for each
π
is negative then the
power of the exponent is also negative and the coefficients
ππ
tend to zero in
time. This means that the function
π
,
(2.79)3
, vanishes in time and it can be
concluded that the perturbed interface moves toward its original equilibrium
position. If for at least one value of
π
the sign of the
πΏch
π
is positive, then
the amplitude
ππ
will grow exponentially in time. Then the solution will be
unstable and the lowest number
π
for which
πΏch
π
is positive will define the
mode of the stability loss.
Analogously, if the stability analysis is performed for a cylinder undergoing
a phase transformation, one can obtain values for
πΏph
π
from the solution of the
perturbed BVP
(2.75)
β
(2.76)
and the variation of the kinetic equation
(2.52)
.
As a result for each number
π
an equation for the amplitude of perturbation
ππcan be written as 1
π
ph
dππ
dπ‘=πΏph
πππ,(2.83)
where
πΏph
π
is a linear integro-differential operator. As in the problem with
chemical interface, the character of the solution is governed by the sign of
the
πΏph
π
. If the sign for each
π
is negative, then the power of the exponent is
also negative and the coefficients
ππ
tend to zero in time. This means that
the function
π
, which describes the divergence from the equilibrium position
of the transformation front, vanishes in time and one can conclude that the
solution is stable. If at least for one
π
the sign of the
πΏph
π
is positive, then the
amplitude ππwill grow exponentially in time, the solution will be unstable.
In the case of cylindrical problem, the expressions for the
πΏπ
βs are very long,
which makes it impossible to relate
πΏch
π
and
πΏph
π
as it was done for the planar
interface in
(2.69)
. Therefore for selected number of material parameters both
stability analyses are done: for phase transformation and for chemical reaction.
This is done to check whether the diffusion can impact the stability.
Section 2.4. Stability of a cylindrical chemical interface
36
(π)
πΏph
π
πeq
(π)
πΏph
π
πeq
(π)
πΏph
π
πeq
(π)
πΏph
π
πeq
Fig. 2.9: πΏπ
-lines for the stability analysis of the cylindrical phase interface for: (
π
)
- stable set of elastic parameters for the solid cylinder, (
π
)- unstable set of elastic
parameters for the solid cylinder, (
π
)- stable set of elastic parameters for the hollow
cylinder with internal radius
π
= 0
.
1
π
,(
π
)- stable set of elastic parameters for the
hollow cylinder with internal radius π= 0.5π
.
Since expressions for
πΏπ
βs do not add any value to the results they are
omitted in the manuscript. However, these complicated equations still can be
solved analytically without numerical methods. To handle them the symbolic
mathematics software Mathematica was used.
In Fig. 2.9 different possible phase transformation front stability scenarios are
shown. The plots (a) and (b) correspond to stable and unstable sets of material
parameters for the solid cylinder, respectively. For a stable configuration
πΏph
π
is negative for each
πβ₯
2and all possible positions of the equilibrium radius.
For the unstable configuration the solution shows positive
πΏph
π
βs. Plots (c)
and (d) were calculated only for βstableβ material parameters and show the
influence of the inner diameter of the hollow cylinder on the stability of the
interface. Small holes, Fig. 2.9(c), lead to an unstable solution only for the
configuration where the equilibrium radius is close to the hole. The limit case,
when the hole diameter tends to zero, gives the solution for the solid cylinder,
as expected. Fig. 2.9(d) shows that for large inner diameters of the cylinder
the solution is unstable even for a stable set of parameters.
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 37
(π)
πΏπ
πeq
(π)
πΏπ
πeq
Fig. 2.10: πΏπ
-lines for the stability analysis of the cylindrical chemical reaction front
for: (
π
)β stable set of elastic parameters for the hollow cylinder with internal radius
π= 0.1π
,(π)β unstable set of elastic parameters for the solid cylinder.
π+πβπ+πβ
Stable configuration 25 45 10 25
Unstable configuration 45 25 25 10
Tab. 2.3:
Material properties and parameters used in analytical and numerical simuξΏ
lation. Material β+β refers to the outer, material β
β
β to the inner region, respectively.
In Fig. (2.10) results are shown for the stable configuration of material
properties from Tab. 2.3, diffusion properties from Tab. 2.2 and energy
parameter
πΎ
= 0
.
15 for solid (a) and hollow (b) cylinders. For all configurations,
the presence of the diffusion in the model does not change the system behavior
compared to the phase transformation problem. The presence of inner hole
has the same effect on the stability behavior as in the phase transformation: a
small inner diameter affects only configurations where the equilibrium radius
is close to inner surface and in the limit of a vanishing radius results in the
solution for the solid cylinder. Large holes make the system unstable.
Using these results, the stable and unstable configuration is selected for
further studies. These values are listed in the Tab. 2.3. These parameters are
used in the following Chapter 3 as a reference values to validate numerical
procedures.
2.5 Remark on phase transition zones
For the problem of phase transformation, which is mathematically pretty
similar to the problem of chemical transformation, there is another necessary
stability criteria. The latter is based on the analysis of so called phase transition
zones. It was proved in Grabovsky and Truskinovsky [2011]; Grabovsky and
Truskinovsky [2013] that if strains at the equilibrium interface belong to the
Section 2.5. Remark on phase transition zones
38
boundaries of the phase transition zones, then this configuration is stable.
However, this approach cannot be explicitly extrapolated to the case of the
chemical reaction.
The procedure of obtaining the process zones was discussed in many papers,
starting from Freidin and Chiskis [1994a]; Freidin and Chiskis [1994b] (more
general nonlinear elastic cases), Morozov and Freidin [1998] (linear elastic case)
and recent works about microstructures and optimal composites Antimonov,
Cherkaev, and Freidin [2016]; Freidin and Sharipova [2019]. In this section,
a short overview of this procedure is presented. Then, the PTZ approach is
utilized in order to validate results obtained with the perturbation method.
For the phase transformation problems, to construct the PTZs one has to
analyze the driving force
π
given in a form
(2.23)
. With this expression, only
the last term depends on the normal vector
n
. For the case of linear isotropic
material, this term is a quadratic form and reads
π¦β(n) = q+:Kβ(n) : q+=1
πβ(οΈπ2βπβ+πβ
πβ+ 2πβ
π2
1)οΈ,(2.84)
where
π1=nΒ·q+Β·n, π2=nΒ·q2
+Β·n(2.85)
are so-called orientation invariants. It was shown that minimization of energy
with respect to the microstructure geometry (i.e., normal
n
direction) leads to
extrema condition of
(2.84)
at given
q+
. The maximum of
π¦β
(
n
)corresponds
to minimum of energy.
For the case of two phase configuration with only one inter-phase boundary,
the direction of the maximizing normal
n*
(
q+
)can be expressed through
the principal values of tensor
q+
. Let
ππππ₯, ππππ
and
ππππ
be the maximum,
intermediate and minimum eigenvalue of
q+
respectively with corresponding
eigenvectors
eπππ₯,eπππ
and
eπππ
;
|π|πππ₯
and
|π|πππ
be the maximum and
minimum absolute value of q+eigenvalues. Then if
ππππππππ₯ <0,or {οΈππππππππ₯ >0
(1 βπβ)|π|πππ < πβ|π|πππ₯
(2.86)
the maximizing normal n*=ππππ₯eπππ₯ +ππππeπππ and
π2
πππ₯ =(1 βπβ)ππππ₯ βπβππππ
ππππ₯ βππππ
, π2
πππ =πβππππ₯ β(1 βπβ)ππππ
ππππ₯ βππππ
.(2.87)
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 39
Otherwise, namely
{οΈππππππππ₯ >0
(1 βπβ)|π|πππ > πβ|π|πππ₯,(2.88)
the maximizing normal has only one component the basis of eigenvectors of
q+
and n*=ππππ₯eπππ₯.
Recalling the case of planar interface described in Section 2.3, the equilibrium
configuration is a horizontal interface with the normal
n
=
eπ¦
. In this problem,
the eigenvectors of
q+
correspond to the Cartesian basis vectors, and
q+
=
π+
π₯eπ₯eπ₯
+
π+
π¦eπ¦eπ¦
+
π+
π§eπ§eπ§
. With this, the horizontal interface at the equilibrium
position is stable only if the condition
(2.88)
is satisfied and
ππ¦
is the maximum
eigenvalue of
q+
. Analogously, for the cylindrical interface from Section 2.4,
only condition
(2.88)
leads to the stable equilibrium interface with the normal
n=eπ.
Phase transition zones built with the use of
(2.86)
and
(2.88)
are shown
in Figs. 2.11 and 2.12 for the planar and cylindrical problems, respectively.
These plots are built for the stable cases of phase transformation equilibrium,
and material parameters are given in Tabs. 2.1 and 2.3. Solid dots in both
figures represent the actual strains at the interface in the equilibrium position.
Dashed lines correspond to the condition
(2.86)
and solid lines β to condition
(2.88).
For selected material parameters strains at the interface belong to the phase
transition zones in the case of planar interface and the minimizing normal
vector has only one component
eπ¦
, so the configuration is stable. It is shown
earlier in the discussions about
(2.69)
that the diffusion does not affect the
stability of planar interface. Therefore, one can conclude that the chemical
interface is also stable.
For selected material parameters strains at the interface belong to the phase
transition zones in the case of cylindrical interface and the minimizing normal
vector has only one component
eπ
, so the configuration is stable. It is not
possible to study analytically the influence of the diffusion on the stability of
the cylindrical interface (see Section 2.4), therefore stability of the cylindrical
chemical interface cannot be proved.
2.6 Conclusions
Problem of a chemical interface propagation in linear elastic solids is formulated
based on the concept of chemical affinity tensor from works by Freidin [2013];
Freidin, Vilchevskaya, and Korolev [2014]; Freidin [2015]; Freidin et al. [2016].
Section 2.6. Conclusions
40
β0.04 β0.02 0.00 0.02 0.04
β0.04
β0.02
0.00
0.02
0.04
ππ¦π¦
ππ₯π₯
Fig. 2.11:
Phase transition zones for the stable set of material parameters from Tab.
2.1. Blue and red dots correspond to the strain state for planar interface configuration.
Dashed lines correspond to the condition (2.86), solid lines β to condition (2.88).
β0.3 β0.2 β0.1 0.0 0.1
β0.3
β0.2
β0.1
0.0
0.1
πππ
πππ
Fig. 2.12:
Phase transition zones for the stable set of material parameters from
Table 2.3. Blue and red dots correspond to the strain state for cylindrical interface
configuration. Dashed lines correspond to the condition
(2.86)
, solid lines β to condition
(2.88).
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 41
An estimation for the total deformation
π½tr
due to the chemical transformation
and the diffusion model is proposed.
A linear stability analysis problem is formulated for the phase transformation
interface in the equilibrium position (based on works by Eremeev, Freidin, and
Sharipova [2003]; Eremeev, Freidin, and Sharipova [2007]) and extended to the
case of the chemical interface. A detailed derivation for of the boundary and
interface conditions, as well as the kinetic equation, is given for the perturbed
phase transformation front and chemical interface.
In this work, the stability problems of a planar and cylindrical chemical
interface were considered analytically using the perturbed kinetic equation
approach. For the planar interface, it was shown that the diffusion does
not affect the stability of the interface. Stability regions were constructed
in material parametric space for that problem. For a cylindrical chemical
interface, two cases were studied: hollow and solid cylinders undergoing a
chemical reaction. It was shown that a small inner diameter does not affect
the stability behavior of the interface when the equilibrium position is far
enough from the hole, while large inner diameters with prescribed boundary
conditions on them make the system always unstable. Due to the complexity
of the perturbed kinetic equation, it was not possible to prove analytically the
influence of the diffusion parameters on the stability of the interface, when
compared to the mathematically similar problem for the phase transformation.
However, the results of the analyses show that for the selected choice of material
parameters diffusion did not change the stability behavior.
Even for the simple geometries and linear elastic materials, the perturbed
stability analysis is quite complicated and cumbersome. Therefore, in the next
chapter, it is studied whether the developed numerical approach for simulating
the chemical reaction (and the phase transformation) front propagation can
reveal stability or instability of the chemical interface.
Based on the results of the stability analysis, two sets of material parameters
are selected for further numerical studies of planar and cylindrical interfaces.
These sets are further referred as βstableβ and βunstableβ for stable and unstable
interface equilibrium position, respectively. An additional check was done for
the stable set of parameters using the PTZ approach. The latter, however, is
only relevant for the phase transformation problem. Its extension to the case
of the chemical reaction is not straightforward and was out of the scope of this
work. The PTZ criteria confirmed the stability of the phase transformation
front for the selected material parameters. As mentioned earlier, for all studied
cases the diffusion did not affect the stability behavior, therefore, it is assumed
that considered stable configurations are stable indeed.
All methods in this chapter considered the stability of the chemical reaction
front at the equilibrium position. The next chapter considers numerical studies
Section 2.6. Conclusions
42
of the chemical interface approaching this equilibrium position for both: stable
and unstable configurations.
Chapter 2. Chemical reaction front kinetics and stability
Numerical and analytical studies of the chemical reaction front kinetics in solids 43
3 Numerical simulation and study of
the chemical reaction front kinetics
There is a large variety of numerical approaches that can handle problems
with moving interfaces. All of them can be divided into two groups relying on
a smooth and a sharp representation of the interface. A typical example of
the smooth-interface approach is the phase-field method, see e.g. Svendsen,
Shanthraj, and Raabe [2018]; Weinberg, Werner, and Anders [2018] in applica-
tion to chemo-mechanics and Schneider et al. [2017]; Schneider et al. [2018] in
application to the phase transformations. These methods require an additional
equation which governs the evolution of the phase field variable. Since the
kinetic equation for the interface is already chosen in the previous chapter, it
might be not convenient to utilize the phase field method to solve considered
type of problems. Moreover, given kinetic equation based on chemical affinity
tensor concept uses jump conditions across the interface. They are much easier
to satisfy for the sharp interface methods, where the thickness of the interface
is neglected. Therefore, the latter was chosen to be implemented using the
finite element method.
All sharp-interface finite-element based methods can be divided into three
subcategories. The first subcategory is when the interface coincides with the
element edges (or faces in 3D case) and the geometry is completely remeshed
each time the interface moves, e.g., Mueller and Gross [1998]; Mueller and
Gross [1999]; Gross, Mueller, and Kolling [2002]; Mueller, Gross, and Lupascu
[2006] for the case of phase transformations and, e.g., Freidin et al. [2016] in
application to chemo-mechanics. The second subcategory also relies on the
interface coinciding with the element edges/faces, however, the mesh in only
distorted as the interface moves, i.e., the nodes are moved, without changing
neither the number of the nodes nor the mesh connectivity, e.g., Morozov
et al. [2018a] in application to chemo-mechanics where such approach has been
implemented using the isogeometric method. The third subcategory unites
approaches where the interface cuts through the finite-element mesh in an
arbitrary way and moves independently of the mesh, which is unchanged from
one time increment to another. Typical examples of such approaches are the
level set method (e.g., Duddu et al. [2011]; Zhao, Bordas, and Qu [2015];
Moghadam and Voorhees [2016]) and the CutFEM, e.g., Hansbo, Larson, and
Larsson [2017]; Burman et al. [2018]; Poluektov and Figiel [2019]. In current
44
work, the remeshing method is utilized for solving problems of moving chemical
interfaces. The used procedure is compared with the isogeometric analysis
and with the CutFEM method, therefore all free sharp-interface finite-element
based methods are covered in this study.
3.1 Description of the numerical procedure
Numerical procedure in the case of small strains and linear elastic materials
explicitly follows the analytical procedure showed in the Section 2.1. All the
assumptions from the previous chapter still hold:
β’A plain strain formulation is considered;
β’Both materials are assumed to be linear elastic;
β’
Diffusion of the particles does not introduce additional stresses (the solid
skeleton approach);
β’
Diffusion process is assumed to be much faster than the chemical reaction
front propagation;
β’
Some initial position of the chemical interface is introduced in the geom-
etry.
Necessity of existence of some thin initial layer of the new material is one of
the weaknesses of this procedure. In order to obtain the solution for stresses,
strains and concentration, one have to get at least one layer of finite elements
for each material domain. This issue is much less strict in the meshless methods
like aforementioned CutFEM. However it has its own challenges, which will be
discussed later in this chapter.
With the given geometry, interface position and boundary conditions, one can
obtain all the required values to calculate the normal component of the affinity
tensor,
π΄ππ
. Indeed, since in the considered class of problems (as described
in Chapter 2) the mechanical and diffusion problems are decoupled, one can
solve them separately one after another. Solution of the first problem gives the
stresses and strains at the reaction front. These values are then substituted
into the boundary condition at the moving interface for the diffusion problem.
Solution of the latter gives the concentration of the diffusing constituent at the
chemical reaction front. Now, the
π΄ππ
can be calculated, and substituted into
the kinetic equation, providing the normal component of the reaction front
velocity. This kinetic equation is solved with the explicit Euler scheme and with
the given time increment one can calculate the normal displacement for each
point of the interface. With the new position of the chemical reaction front,
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 45
Start
Input geometry and material
properties. Set initial position of the
interface (R0), Ξπ‘and π= 0
Generate the FE mesh for this
iteration and find normals nπat
each node of the interface
Calculate ππand ππfor each node
of the interface, obtain ππ(ππ,ππ)
Calculate ππfor each node at the
interface, obtain π΄π
ππ(ππ,ππ, ππ)
Ξπ’π=π*ππ(οΈ1βexp (οΈβπ΄π
ππ
π
π )οΈ)οΈ,
Rπ+1 =Rπ+ Ξπ’πnπ
Finish
New iteration?
π=π+ 1
yes
no
Fig. 3.1: The flowchart for the numerical algorithm.
the procedure starts over with the next iteration. This process is summarized
in the flowchart shown in Fig. 3.1.
Since the used mathematical model provides the expression for the normal
Section 3.1. Description of the numerical procedure
46
component of the reaction front velocity, one has to accurately calculate the
normal vector
n
to the reaction front. This is one of the main reasons why on
the first stages of this research, the isogeometric FEM procedures (Hughes,
Cottrell, and Bazilevs [2005]) were utilized. Problem formulation and main
results can be found in Morozov et al. [2018a]. By using Non-Uniform Rational
B-Splines (NURBS) as basis functions, one can calculate with high accuracy
the direction of the normal vector for any external or internal edge (curve)
of a geometrical domain, and particularly for the chemical reaction front. In
addition, IGA is more accurate with less computational effort compared to
standard FEA, cf. Morganti et al. [2015]. Comparison of the IGA and classical
FE models for the cylinder problem is shown in Fig. 3.2.
(a) (b)
Fig. 3.2:
Geometry described by NURBS curves which is directly used for the numerical
simulation (a). Classical FE model for the numerical simulation (b).
Since IGA uses displacements of control points as degrees of freedom which
are not interpolative inside a domain, one cannot move the internal front
directly in general case (not for simple translation and rotation). Therefore,
a special procedure should be introduced to apply the displacement to the
interface and βre-meshβ the domain appropriately. To do so, an additional
elasticity problem is solved with Ξ
π’
applied at the interface boundary as
non-homogeneous Dirichlet boundary conditions. Nitscheβs method is used for
imposing such kind of boundary conditions weakly (Hansbo [2005], see also
Juntunen and Stenberg [2009] and references therein). Doing that, the elements
are simply distorted without changing the topology or mesh connectivity. This
new distorted mesh defines the new position of all the nodes (the whole
isogeometric model) and can be used as an input to the next iteration. Number
of elements remain the same during all analysis. This can be handled without
loss of the numerical accuracy because higher distortion of the elements is
allowed in the isogeometric method when compared to the standard FEM.
Despite all advantages, iterative update of the interface position may become
very complicated when the geometry is defined by NURBS, especially when
the interface between the domains has to be extrapolated or trimmed after the
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 47
ππ
ππ
nπ
nπ
n
Fig. 3.3:
Definition of the normal vector to the piecewise linear interface on the FE
mesh
movement. It was decided that working with more complex geometries and
loading scenarios is more important than computational time, therefore, for
the later studies the classical FEM procedure was utilized.
Traditional finite elements defines the interface as a piecewise linear curve
which makes it impossible to determine the normal uniquely. Usually, the
normal vector at a node is defined as a weighted sum of normal vectors to the
element edges, Fig. 3.3,
n=1
βοΈπ2
π+π2
π+ 2ππππnπΒ·nπ
(ππnπ+ππnπ),(3.1)
where
nπ
and
nπ
are the normal vectors to the element edges,
ππ
and
ππ
are the
lengths of these edges.
As mentioned earlier, in some configurations when moving the interface it
has to be adjusted, as shown in Fig. 3.4. The interface (the blue curve) may
move inside (left boundary) or outside (right boundary) the body domain. The
movement in normal direction is shown by black arrows and the new position
of the interface is indicated by a black dashed line.
In order to avoid this unphysical configuration, where the start and end of
the reaction front do not belong to the boundary of the body domain, a special
procedure was created. It cuts and linearly extrapolates the interface to adjust
the boundaries as shown by red curve in Fig. 3.4. In the analysis, this adjusted
curve is considered as a new position of the interface.
The numerical analysis was performed with the commercial finite element
program Abaqus. The problem was solved quasi-statically as described in
Freidin et al. [2016]; Morozov et al. [2018a]; Morozov et al. [2018b], where the
position of the interface defines the time dependence. A special Python script
Section 3.1. Description of the numerical procedure
48
Fig. 3.4: Automated adjustment of the interface position during its propagation.
for Abaqus was developed to automate the process flow and postprocess the
results.
3.2 Kinetics of the cylindrical chemical reaction front
approaching the blocking state
A cylinder undergoing the chemical reaction is considered, as described in
Section 2.4 and in Fig. 2.6. Due to symmetry, only a quarter of the cylinder is
modeled. A thin layer of transformed material was assumed in the initial state
for the analysis. Material parameters for stable and unstable configurations
were chosen for the finite element model in accordance with Tab. 2.3. In order
to compare the kinetics near the thermodynamic equilibrium, the external
boundary conditions
π’0
were adjusted so that the equilibrium radii are the
same for both problems.
Fig. 3.5 (b) shows how the radius of the interface changes in time. The
stable configuration converges smoothly to the analytically predicted equilib-
rium radius. For the unstable configuration, the blue line, change of radial
coordinates are plotted for twenty equally spaced reference points distributed
along the interface (highlighted in Fig. 3.5 (a)). Note that during the first
fifteen increments, the displacements of these points are equal and the phase
transformation front keeps its originally circular shape. However, when the
front approaches the equilibrium radius, it changes its smooth behavior and
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 49
(a)
πππ‘ππππ‘πππ
π
(b)
Fig. 3.5:
(a) β set of reference points on the interface for the output. (b) β results of
numerical simulations for stable (red curve) and unstable (blue curve) configurations
for the phase transformation front propagation problem.
Fig. 3.6:
Typical shape of the phase interface after losing stability near the equilibrium
radius during numerical simulation.
the loss of stability is observed. The typical shape of the interface right after
losing the stability is shown in Fig. 3.6. This mode of stability loss is governed
by the numerical accuracy and, therefore, depends on the mesh discretization.
The instability shape of the initially circular interface on the FE mesh is shown
in Fig. 3.7.
The stable and unstable configurations behave in the same manner for phase
and chemical transformation fronts. The interface approaches the equilibrium
position similarly. Fig. 3.8 shows how the radius of the interface changes in time
during the chemical reaction. For the unstable configuration, the blue line, the
Section 3.2. Kinetics of the cylindrical chemical reaction front approaching the blocking state
50
Fig. 3.7:
Shape of initially circular interface after the loss of stability on the FE mesh.
change of radial coordinates are plotted for twenty points equally distributed
along the interface. During the first twenty increments, the displacements
of these points are equal and the chemical reaction front keeps its originally
circular shape. As it approaches the predicted equilibrium radius, it changes
the smooth behavior and stability loss occurs. If the initial shape of the
reaction front is circular, the shape of the unstable interface is governed by the
numerical inaccuracy as it was in the case of phase transformations. One can
predefine the mode of the instability by introducing wave-type perturbations
with small amplitude and different frequency into the initial configuration.
In order to present kinetics of the interface approaching the equilibrium
radius from inside and from outside, the artificial case of so-called βbackward
reactionβ is analyzed. The reverse reaction mechanism might differ from the
direct reaction, but here it is assumed to be the same, just to test the numerical
procedure. Therefore, the initial chemical interface radius is set closer to the
center of the cross-section then the analytically predicted equilibrium radius.
For simplicity, it is assumed that during the simulation the gaseous constituent
diffuses from the chemical interface to the outer boundary of the body and the
material β+β transforms back to the material βββ.
Kinetics of the back and direct reactions with initially perturbed front is
shown in Fig. 3.9 for unstable configuration. Top row corresponds to the back
reaction and bottom - to the direct reaction. Left and right columns have
different frequency of the initial perturbation. In the case of stable material
parameters, the perturbations, as well as the ones introduced by numerical
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 51
πππ‘ππππ‘πππ
π
Fig. 3.8:
Results of numerical simulations for stable (red curve) and unstable (blue
curve) configurations for the chemical reaction front propagation.
inaccuracy, vanish and stability occurs at the equilibrium radius. But for
another set of parameters, the perturbations amplify when the reaction front
comes to the equilibrium position, so that one can see the predefined mode of
the instability as shown in Fig. 3.9. This effect appears when the amplitude of
the initial perturbation is greater then the accuracy of the numerical scheme.
3.3 A note about stresses caused by the interface
instability
The growing perturbations in the case of unstable configuration result in the
redistribution of the stress fields. Examples of the von Mises stress fields for
circular and perturbed interfaces are shown in Fig. 3.10. For given example,
stresses at the perturbed interface are much higher than on circular interface
and, therefore, may lead to the plasticity or to some failure scenarios, like
delamination or fracture.
Stresses caused by the instability amplitude growth were analyzed in Morozov,
Freidin, and MΓΌller [2019]. In this section, a static configuration is considered,
in which the chemical interface has a predefined sinusoidal-type shape near the
blocking state, as shown in Fig. 3.11. This position of the interface is similar
to the one shown in Fig. 3.9, e.g., top right. The unstable set of parameters
from Tab. 2.3 is used for the numerical analysis.
Section 3.3. A note about stresses caused by the interface instability
52
Fig. 3.9:
Kinetics of the initially perturbed interface with different predefined modes of
instability for the unstable set of material parameters. Top row - βbackward reactionβ,
bottom row - direct reaction.
(a) (b)
Fig. 3.10:
The von Mises stress distribution for stable (a) and unstable (a) configuraξΏ
tions.
Various amplitudes and frequencies for the interface shape were analyzed.
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 53
Fig. 3.11:
Unstable configuration of a solid cylinder with predefined shape of the
interface. Stress distributions are analyzed along radial lines OA and OB.
Fig. 3.12: Distributions of von Mises stress along corresponding directions
Stress distributions along the radial lines OA and OB (Fig. 3.11) are shown in
Figs. 3.12β3.14.
The maximum von Mises stress in OB direction increases with the amplitude
and frequency of the interface shape. As shown in the previous section, the
amplitude grows with the propagation of the front in the unstable configuration.
This may lead to intensive plastic deformations.
Section 3.3. A note about stresses caused by the interface instability
54
Fig. 3.13: Distributions of the hoop stress along corresponding directions
Fig. 3.14: Distributions of the radial stress along corresponding directions
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 55
One should note that negative hoop stress (Fig. 3.13) in both OA and OB
directions do not change their sign in the vicinity of the stress concentration
due to the amplitude growth of the reaction front perturbation. However, an
increasing magnitude of the radial stress in the OB direction (Fig. 3.14) may
lead to delamination. In this case, not only the stress distribution but also the
conditions for the chemical reaction will change.
3.4 Comparison of the CutFEM and remeshing results
3.4.1 CutFEM overview
CutFEM belongs to subcategory of the numerical approaches where the in-
terface cuts through the finite-element mesh in an arbitrary way and moves
independently of the mesh. The mesh remains unchanged from one time
increment to another. Similar numerical approaches based on the same ideas
are the combination of the extended finite-element method (XFEM) to solve
the PDEs and the level-set method to move the interface, e.g., Zhao, Bordas,
and Qu [2013]; Zhao et al. [2013] in application to the phase transformations
and Duddu et al. [2011]; Zhao, Bordas, and Qu [2015] in application to chemo-
mechanics. Recently, the CutFEM method has been adapted for problems of
mechanochemistry by Poluektov and Figiel [2019], where the numerical method
has been formulated for the most general finite-strain mechanochemical setting,
i.e., involving non-linear PDEs. The CutFEM method has been originally
formulated for linear problems Burman and Hansbo [2012] and later adapted
specifically to linear elasticity Burman et al. [2018]; Burman et al. [2019a];
Burman et al. [2019b].
The CutFEM method relies on two main features: enforcement of the in-
terface conditions weakly using the Nitsche method, which allows solving the
discretized PDEs with the interface cutting through the elements, and intro-
duction of an inter-element stabilization, which addresses the ill-conditionality
of the discrete problem related to the interface partitioning the elements into
highly unequal spatial fractions. In Poluektov and Figiel [2019], the weak
form for the finite-strain coupled mechanochemical problem has been derived
from the variational principle, where problem-specific interface conditions (e.g.,
the force equilibrium and the displacement continuity for mechanics) and
inter-element stabilization terms have been used.
In this section the comparison of CutFEM and remeshing methods is per-
formed for linear elasticity. Therefore, the CutFEM-based method of Poluektov
and Figiel [2019] has been simplified to a linear elastic setting. This step is
out of scope of this work, therefore, the details are omitted.
The CutFEM-based and the remeshing methods are compared in three differ-
ent examples, which are considered in the following sections. The first example
Section 3.4. Comparison of the CutFEM and remeshing results
56
is the propagation of a planar reaction front, with an initially introduced
perturbation. Two scenarios are considered: physically stable and physically
unstable behavior. The second example simulates the interface kinetics in the
specimen under shear load. The third example focuses on a different topology
of the reaction front, namely a closed curve. All mentioned examples are solved
on a 2D geometry with a plane-strain formulation, as it was done in previous
sections.
3.4.2 Stable and unstable configurations for planar chemical
interface
As mentioned above, to investigate the performance of the methods, physically
stable and physically unstable scenarios are considered. Geometry of the
problem is similar to one described in Section 2.3 and shown in Fig. 2.3.
However in this simulation the planar reaction front has an initial
cos
-type
perturbation.
For the mechanical problem, the following boundary conditions are used:
clamped bottom boundary, symmetry conditions on left and right boundaries,
prescribed displacements on top boundary (vertical displacement is
π’0
, hori-
zontal displacement is zero). For the diffusion problem, the following boundary
conditions are used: mixed boundary conditions on bottom boundary and
interface (according to (2.22)), zero flux at left and right boundaries.
For the geometry dimensions
π»
=
πΏ
= 1 are used. The initial position of
the perturbed interface is taken to be curve
π¦
= 0
.
1 + 0
.
002
cos (6ππ₯)
. For the
physically stable case, LamΓ© parameters of the materials and transformation
strain are given in Tab. 2.1. External displacement was set to
π’0
= 0
.
0453.
For the physically unstable case, the parameters are also taken from Tab. 2.1,
but the applied displacement is
π’0
=
β
0
.
0381. One should note, that external
loads are different to adjust the equilibrium position to be in the middle of the
layer at
β
= 0
.
5. As it was shown in Section 2.3, diffusion parameters do not
influence the stability of the interface. Therefore, only one set of the chemical
and diffusion parameters were used and they are listed in Tab. 2.2. Energy
parameter is taken πΎ= 0.05.
In Fig. 3.15, for the physically stable case, the time-evolution of the reaction
front is shown by plotting the
π¦
-coordinate of three different points of the
interface, with coordinates
π₯
= 1
/
3,
π₯
= 1
/
2,
π₯
= 2
/
3, for both CutFEM and
remeshing solutions. It can be seen that for the most part of the simulation
time, all six curves are indistinguishable. At that stage, the absolute difference
between
π¦
-coordinates of the points resulting from two methods is in the range
of 10
β4
to 10
β6
. However, when the velocity of the front becomes slow, near
the equilibrium position, the CutFEM-based approach produces an artifact:
the front rapidly aligns with the nearest element edges. This happens due
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 57
Time, π‘
β|π₯
Fig. 3.15:
Results of numerical simulations by CutFEM and by the remeshing proξΏ
cedures. Propagation of a planar chemical reaction front for stable configuration is
illustrated by the evolution of π¦-coordinate of three points of the front.
to numerical error of the stresses becoming dominant in the reaction front
driving force at slow front velocities and is the major disadvantage of using this
method for modeling an approach to the equilibrium. This issue arises from
the use of linear triangular elements and the calculation of the finite-element
stresses, which are constant within such elements. Therefore, to correct this
artifact in th CutFEM procedure, either a higher-order elements must be used,
or improved stress-calculation procedures must be involved.
In Fig. 3.16, for the physically unstable case, the time-evolution of the
reaction front is shown by snapshots of the front configuration at four different
moments of time. It can be seen that the discrepancy between the results
obtained by numerical methods accumulates with time and is mostly revealed
at the boundaries of the domain. The accumulation in time is related to
physical instability of the interface: as any perturbation of the front should
grow in time, a numerical perturbation (i.e. a numerical error) also grows.
The discrepancy at the edges of the domain is related to slightly different
extrapolation procedures in both approaches, the edge points of the interface
are allowed to move inside the domain (as shown in Fig. 3.4). In this case
the extrapolation is used to find new intersection points of the interface and
the boundary of the domain. For the CutFEM-based method, this has been
described in Poluektov and Figiel [2019] and for the remeshing it is shown in
Section 3.1.
For CutFEM, the mesh consists of linear elements in the form of isosceles
right triangle with the side of Ξ
π₯
. In the remeshing procedure, linear quads
with full integration were used, while the size of an element was approximately
Section 3.4. Comparison of the CutFEM and remeshing results
58
π‘= 9600 π‘= 13600
π‘= 17600 π‘= 21600
Fig. 3.16:
Results of numerical simulations by CutFEM and by the remeshing proceξΏ
dures. Kinetics of the initially perturbed interface for the unstable configuration is
illustrated by four snapshots of the reaction front at different times.
equal to Ξ
π₯
. For this example Ξ
π₯
= 1
/
64 was taken. Time steps of Ξ
π‘
=
80 and Ξ
π‘
= 160 were taken for the physically stable and unstable cases
respectively. For CutFEM, numerical parameters, which were denoted as
π
and
π
in Poluektov and Figiel [2019], were taken 10
4
and 10 respectively. Examples
of the meshes are shown in Fig. 3.17.
As shown in Fig. 3.15, in the case of stable equilibrium position of the
interface, the perturbation diminishes as the front approaches the equilibrium
position. When the selected elastic material constants correspond to equilibrium
being physically unstable, Fig. 3.16, the amplitude of the perturbation rapidly
(exponentially, as shown analytically in Section 2.3) grows, even before the front
approaches the equilibrium position. From this, one can conclude that both
methods reveal the instability effect and adequately simulate stable scenarios.
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 59
Fig. 3.17:
Cropped examples of the finite-element meshes used in CutFEM and
remeshing procedures are shown for the snapshot at π‘= 17600 in Fig. 3.16.
Therefore, both of them can be used for modeling the reaction front kinetics
approaching stable and unstable equilibrium position.
3.4.3 Stable configuration for the chemical interface under shear
The second considered case is similar to the first example, however, an additional
shear displacements are applied at the top boundary. Therefore, the model
setup and the boundary conditions are same, except prescribed displacements
on top boundary:
u
=
π’0eπ¦
+0
.
01
eπ₯
. This creates a shear stress state, therefore,
if the initial position of the interface is horizontal the interface should rotate,
as it approaches the stable equilibrium position. Exactly this is observed in
the results of numerical simulations shown in Fig. 3.18, where snapshots of
the front configuration at four different moments of time are plotted.
The time-evolution of the reaction front is shown in Fig. 3.19 by plotting the
π¦
-coordinate of three different points of the interface with coordinates
π₯
= 1
/
5,
π₯= 1/2,π₯= 4/5, for both CutFEM-based and remeshing approaches.
Since the reaction front rotates, point
π₯
= 1
/
5moves slower and lags behind
point
π₯
= 1
/
2, while point
π₯
= 4
/
5overtakes point
π₯
= 1
/
2. Also, as in the
previous example, during most of the simulation time, when the velocity of the
front is relatively large, both numerical approaches give indistinguishable curves.
However, as the velocity drops with the front approaching the equilibrium
position, the CutFEM-based method results in the front aligning with the
nearest element edges. As the equilibrium position of the front is an inclined
curve, while element edges are either vertical, horizontal or inclined by 45,
the front acquires a stair-like shape. The remeshing procedure gives a smooth
enough piecewise-linear representation of the equilibrium position. Although
this is a noticeable drawback of the existing CutFEM-based approach, it can
Section 3.4. Comparison of the CutFEM and remeshing results
60
π‘= 4000 π‘= 16000
π‘= 40000 π‘= 60000
Fig. 3.18:
Propagation of the initially planar chemical reaction front in a body under
shear loading. The kinetics is illustrated by four snapshots of the reaction front at
different times.
be clearly seen that the front position numerical error close to the equilibrium
position has an order of an element size.
3.4.4 Closed interface in the square domain
The third example focuses on a different topology of the reaction front, namely
a closed curve. In Section 2.4 cylindrical geometries were considered and
a set of parameters leading to stable configuration of the circular reaction
front have been established. In this section, the geometry is changed to a
square, which creates inhomogeneous stress distributions (with respect to the
polar angle in the polar coordinate system), and, therefore, leads to a more
complex equilibrium configuration. Diffusing reactant is supplied trough the
external boundaries and the initial position of the interface was a circle in the
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 61
Time, π‘
β|π₯
Fig. 3.19:
Propagation of the initially planar chemical reaction front in a body under
shear loading. The kinetics is illustrated by the evolution of
π¦
-coordinate of three
points of the front.
center of this square. To highlight the effect of stresses on the equilibrium
configuration of the reaction front, two different loading cases are considered:
biaxial stretching and piece-wise linear prescribed displacement. The last one
is chosen in a form
u0=π0
2π₯
πΏeπ₯+π0
2π¦
π»(οΈ1 + ξξξξ
2π₯
πΏξξξξ)οΈeπ¦,at π¦=Β±π»
2,
u0=2π₯
π»π0(οΈ2βξξξξ
2π¦
π»ξξξξ)οΈeπ₯+ 2π0
2π¦
π»eπ¦,at π₯=Β±πΏ
2,
(3.2)
where the origin of the coordinate system is assumed to be in the center of the
square. Graphical interpretation of the loading scenarios is shown in Fig. 3.20.
Before the analysis it was assumed that change of the main domain from
circular to square (to be precise, the solution domains are a square and
circle, but the entities they are modeling are a prism and cylinder) would not
change drastically neither position of the equilibrium interface nor its stability.
Therefore, for material parameters, a stable set for cylindrical problem from
Tab. 2.3 was used. The transformation strain was taken to be
π
= 0
.
05. The
same chemical and diffusion parameters were used as in the first example
(Tab. 2.2), but with another value for the energy parameter
πΎ
= 0
.
15. For
the diffusion problem, the mixed boundary conditions were prescribed on all
Section 3.4. Comparison of the CutFEM and remeshing results
62
(a) (b)
Fig. 3.20:
Considered loading cases: biaxial stretching (a) and piecewise-linear preξΏ
scribed displacement (b).
boundaries. For this example height and width of the square
π»
=
πΏ
= 2 are
taken. The radius of the interface initial position is set to 0
.
73. External loading
for biaxial stretching is applied to all the exterior edges and the amplitude is
set to
π’0
= 0
.
076. For the second loading case, the loading parameter
π0
in
(3.2) is set to 0.025.
The initially circular configuration of the front evolves then into two different
shapes for these two loading scenarios, which are Shown in Figs. 3.21 - 3.22.
As in the previous examples, during initial stage (fast kinetics) both methods
produced indistinguishable results, therefore, the evolution of the reaction front
is not shown. At the final stage, i.e. close to equilibrium configuration, the
CutFEM-based method aligns with the element edges. Here, again, it is clearly
seen that the discrepancy in the front position has an order of an element size.
For solving these problems time step of Ξ
π‘
= 50 is taken. For the CutFEM-
based approach, spatial step of Ξ
π₯
= 1
/
32 was taken, for the remeshing
approach the average size of element on the interface Ξ
π₯
= 0
.
0116 is taken for
the biaxial stretching case and Ξ
π₯
= 0
.
0077 β for the second loading scenario.
Cropped examples of the finite-element meshes are shown in Fig. 3.23 for the
second loading case at time π‘= 6100 (Fig. 3.22).
All the three simulated examples reveal that both numerical procedures show
qualitatively similar results. For the planar interface numerical simulations cor-
relate with the analytical predictions. Moreover, both CutFEM and remeshing
simulations shows same kinetic features for stable and unstable configurations.
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 63
π‘= 2000 π‘= 15000
Fig. 3.21: Intermediate and final configurations of a closed-curve reaction front in a
2D body under biaxial stretching.
3.5 Kinetics of a chemical reaction front in a body with
a pore or an inclusion
Another example where the influence of the mechanical stresses on the kinetics
of the chemical interface can be clearly seen is a problem with the interface
approaching some inclusion. In the current section, a chemical reaction in
a finite layer with a cylindrical pore (or void) and with a cylindrical rigid
inclusion is considered, Fig. 3.24.
The reactant is supplied through the lower boundary. Left and right edges
are traction free. External displacements
u
=
π’0eπ¦
are applied to the upper
boundary, and lower boundary is fixed. The material parameters used in the
following numerical simulation are given in Tab. 2.1, the stable configuration
is considered. Diffusion parameters are the same as in Section 3.4.2. For this
particular problem
π»
= 1,
πΏ
= 2, the radius of the pore (or of the inclusion) is
π
= 1 and the initial thickness of the transformed layer is
β
= 0
.
1. An external
displacement is set to π’0= 0.0453.
Due to symmetry, only half of the model is considered. For both problems
time step of Ξ
π‘
= 80 and spatial step of Ξ
π₯
= 1
/
64 are taken. The kinetics of
the initially planar interface approaching the circular pore and solid inclusion
Section 3.5. Kinetics of a chemical reaction front in a body with a pore or an inclusion
64
π‘= 3000 π‘= 6100
Fig. 3.22: Intermediate and final configurations of a closed-curve reaction front in a
2D body under complex loading.
(a) (b)
Fig. 3.23:
Cropped examples of the finite-element meshes for final configuration of a
of a closed-curve reaction front in a 2D body under complex loading. (a) β CutFEM
results, (b) β remeshing results.
is shown in Fig. 3.25, (a) and (b), respectively. Similar setup for a uniform
layer (without the pore and inclusion) is shown in Fig. 3.25 (c). Dashed lines
correspond to intermediate positions of the interface at regular time intervals,
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 65
π₯
π¦
β
0
π
πΏ
π»
Fig. 3.24: Model for the finite layer with circular pore or inclusion.
which are equal to 20Ξ
π‘
for all plots. One can see that the void accelerates
the front propagation, and the inclusion retards it.
(a) (b) (c)
Fig. 3.25:
Kinetics of the initially planar chemical interface in the vicinity of a
cylindrical pore (a), a rigid inclusion (b), and in a uniform layer (c).
The distributions of the von Mises stress in the layer with a pore and with a
rigid inclusion are shown in Figs. 3.26 and 3.27, respectively.
The experiments show that pores (or voids) or another material insertions
may significantly influence the kinetics of the chemical reaction front prop-
agation. As an outlook for further study of the interface passing through
the region with inclusions, an additional procedure for adjusting the moved
interface might be developed, similar to one described in Fig. 3.4.
Section 3.5. Kinetics of a chemical reaction front in a body with a pore or an inclusion
66
Fig. 3.26:
The von Mises stress distribution in the vicinity of a pore for the reaction
front position shown in the right inset.
Fig. 3.27:
The von Mises stress distribution in the vicinity of a solid inclusion for the
reaction front position shown in the right inset.
3.6 Conclusions
An approach to study the influence of mechanical stresses on the chemical
reaction front propagation based on the chemical affinity tensor concept was
implemented in a numerical procedure.
From the correspondence of the analytical predictions and the results of
finite element simulations, it becomes evident that the numerical procedure
can be reliably used for the simulation of reaction and phase transformation
front propagation and a stability check of the interface in thermodynamic
equilibrium. The instability of the reaction front propagating towards the
unstable equilibrium was simulated and analyzed for initial states with a
not perfectly circular (or planar) interface by using the proposed numerical
procedure. Also, it was shown numerically that growing perturbations in the
case of unstable configuration may lead to plasticity and failure.
For the cylindrical interface, numerical simulations show that for the stable
set of material parameters both round and perturbed initial interface smoothly
Chapter 3. Numerical simulation and study of the chemical reaction front kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids 67
converge to the circular equilibrium position. For the unstable set of parameters
two scenarios are realized:
(i).
An initially round interface keeps its shape while propagating toward
thermodynamic equilibrium. It can be concluded that the global kinetics
of the reaction front interface suppresses the process of stability loss.
When the interface approaches equilibrium, its velocity decreases, so that
it becomes comparable with the growth speed of the instability amplitude.
At this point, the shape of the unstable front becomes visible.
(ii).
The amplitude of the initial perturbations of the interface grows with
the front propagation toward equilibrium. Numerical simulations show
that this process accelerates while the interface approaches equilibrium.
Note that for the studied parameters of the initial front perturbations,
the interface keeps the frequency of the predefined perturbations during
its propagation.
The proposed procedure was compared with the CutFEM approach to model
numerically the reaction front propagation. Both CutFEM and remeshing
methods show the same kinetics of the reaction front when it is far from the
equilibrium position. However, when the interface approaches the reaction
blocking state, the CutFEM approach produces an artifact solution, in which
the energy minimizing configuration forces the interface to align with the
nearest element boundary. This issue can be resolved by introducing a more
accurate procedure of calculating stresses at the interface, or by using the
second-order finite elements. When compared with CutFEM, the remeshing
procedure requires tracking the position of the interface, writing additional
scripts to handle geometry- and self-intersections of the front. Nevertheless,
both studied methods can be used to model the chemical reaction front kinetics.
A possible outlook for the remeshing procedure might be upgrading the
method to the case of large deformations and non-linear material constitu-
tive relations. This is already realized in CutFEM (Poluektov and Figiel
[2019]), while the remeshing procedure is rather limited by the Abaqus software
capabilities.
Having the remeshing procedure approbated and validated, it can be used
to solve a real engineering problem.
Section 3.6. Conclusions
Numerical and analytical studies of the chemical reaction front kinetics in solids 69
4 Experimental and theoretical studies
of Cu-Sn intermetallic phase growth
In this Chapter the growth of Inter-Metallic Compound (IMC) layers is consid-
ered: after soldering an IMC layer appears and establishes a mechanical contact
between eutectic tin-silver solder bumps and Cu interconnects in microelectronic
components. Intermetallics are relatively brittle in comparison with copper
and tin. In addition, IMC formation is typically based on multi-component
diffusion, which may include vacancy migration leading to Kirkendall voiding.
Consequently, the rate of IMC growth has a strong implication on solder joint
reliability. Experiments show that the intermetallic layers grow considerably
when the structure is exposed to heat. Mechanical stresses may also affect
intermetallic growth behavior. These stresses arise not only from external
loadings but also from thermal mismatch of the materials constituting the
joint, and from the mismatch produced by the change in shape and volume
due to the chemical reactions of IMC formation. This explains why in this
work special attention is being paid to the influence of stresses on the kinetics
of the IMC growth.
This chapter starts with a report of experimental findings regarding the IMC
growth at the interface between copper pads and tin based solder alloys in
different microchips during a high temperature storage test. Then the growth
kinetics is analyzed by means of a continuum model. By combining experiment,
theory, and a comparison of experimental data and theoretical predictions the
values of the diffusion coefficient and an estimate for the chemical reaction
constant are found. A comparison with literature data is also performed. This
chapter contains an overview of the results obtained in Morozov et al. [2018b];
Morozov et al. [2018a]; Morozov et al. [2020].
4.1 Overview on intermetallic compound growth
The main technological process for creating an electrical contact between
components in a (micro-) electric circuit is soldering. Different eutectic tin-
based alloys containing various metallic chemical elements (Zn, Bi, Pb, Ag,
Cu, etc.) are used for this process. One of the most common lead free solder
alloys in microelectronics is eutectic SnAg3.7 (or ternary SnAg3.6Cu0.8). This
70
alloy is also used in Ball-Grid Array (BGA) components in the microelectronic
industry for solder bumps and paste (Lee and Mohamad [2013]).
Cu6Sn5
Cu3Sn
Fig. 4.1: Cu-Sn phase diagram, reprinted from FΓΌrtauer et al. [2013].
An intermetallic, also known as an intermetallic compound, intermetallic
alloy, ordered intermetallic alloy, and a long-range-ordered alloy is a type of
metallic alloy that forms a solid-state compound exhibiting defined stoichiome-
try and ordered crystal structure, see, e.g., Callister Jr. and Rethwisch [2010].
In electric circuits Cu is used as a conductive material for contacting the
surfaces of the electronic components. During soldering the solder material
melts and gets in contact with the copper substrate so that a thin layer of a
particular IMC forms at the interface between the Cu and the tin-based solder.
In the absence of an IMC layer the bond between the solder and the substrate
is weak, since there is hardly no interaction between the metals at the boundary.
In Sn-Cu and Sn-Ag-Cu eutectic alloys the electrical and mechanical contact is
established by one (Cu
6
Sn
5
) or two (Cu
6
Sn
5
and Cu
3
Sn) intermetallic phases,
respectively (e.g., Liashenko, Gusak, and Hodaj [2014]). Their formation is
determined by the temperature regime, according to the phase diagram shown
in Fig. 4.1. It involves three stages: dissolving of Cu in liquid Sn, a chemical
reaction between the components, and further crystallization. Further growth
of the IMC layer takes place in the solid state. The formation of IMCs occurs
according to the following chemical reactions between copper and tin
6 Cu + 5 Sn βCu6Sn5,
3 Cu + Sn βCu3Sn.(4.1)
Phase Cu6Sn5is also referred to as the πand Cu3Sn as the π-phase.
As mentioned earlier, the existence of the IMC is a necessary condition for
the electro-mechanical contact. However, when compared to pure copper or tin
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 71
the IMC is more brittle. This may lead to a decrease of the reliability of the
joint. In addition, intense diffusion of Cu (from the pad) or Sn atoms through
the IMC may lead to void formation due to the Kirkendall effect, see, e.g.,
Paul [2004]; Dybkov [2010].
There are many experimental studies of the IMC growth in solders of the
aforementioned two and other phases of more complex stoichiometry (Kim,
Huh, and Suganuma [2003]; Yu et al. [2005]; Min-Suk, Chan-Jin, and Hyuk-
Sang [2008]), such as microrelief mapping of the interface, the determination
of the chemical composition of the phases, or the influence of the soldering
time on the kinetics of the IMC formation, etc. In these systems the growth
velocities of the IMC interfaces are much higher than in homogeneously phased
solders.
A typical empirical kinetic equation for the IMC phase growth is based on
fitting the experimentally observed data and has the form (Cogan et al. [1984];
Dariavach et al. [2006]):
β=β0+ππ‘1/π,(4.2)
where
β
and
β0
are the current and initial new phase thicknesses,
π‘
is the
time,
π
and
π
are growth constants. Under the assumption that the major
contribution to IMC growth is bulk diffusion, the power law
(4.2)
takes the
form of a square root dependency (see, e.g., Gao et al. [2006]; Gao et al. [2019]).
In other words, any deviation from the square root behavior indicates that not
only bulk diffusion is governing the growth kinetics (Cogan et al. [1984]). In this
case the growth constant
π
may be expressed through the diffusion coefficients,
π·
, (Mei, Sunwoo, and Morris [1992]). The temperature dependence of the
latter can be expressed by an Arrhenius equation,
π·=π·0exp (οΈβπ
π
π )οΈ,(4.3)
where
π·0
is a pre-exponential factor,
π
is the reaction activation energy,
π
is the gas constant, and
π
is the absolute temperature in
πΎ
, see, e.g., Ross,
Vuorinen, and Paulasto-KrΓΆckel [2016] and the references therein. Based on
these ideas a model for intermetallic growth in thin Sn joints between Cu
substrates was proposed in Arafat et al. [2020] with application to solder
microjoints.
An understanding of the IMC formation process and predicting its kinetics
in various range of temperatures is essential for evaluating the structural
integrity of solder interconnects. The electrical current during operation of the
microelectronic device can heat up a solder bump to 100
Β°
C - 150
Β°
C. This heat
stimulates IMC growth and as a result reduces the lifetime of the joint.
Section 4.1. Overview on intermetallic compound growth
72
4.2 Experiment overview
4.2.1 Specimen preparation
Two different types of microchips (referred to as Series I and II in Fig. 4.2
and in what follows) with BGA packages were available. The solder balls were
approximately of the same diameter, 500
π
m. The frames of the microchips
were made of plastic or metal for Series I or II, respectively.
Series I Series II
1500πm3600πm
Fig. 4.2: Photos of the microprocessor boards.
For the commercially obtained packages the exact information about the
solder material, the package substrate material, and its surface finish materials,
as well the ball attachment process was not available. Clearly, all of this
affects the initial IMC formation, as well as the diffusion processes during
its growth. Therefore, in order to obtain a somewhat clearer picture, the
following was assumed. According to performed EDXS (Energy-Dispersive
X-ray Spectroscopy) analysis, no other elements were found in the IMC region
apart from Cu and Sn. Hence it is fair to assume that the solder material was
eutectic tin silver (AgSn), and the substrate had an OSP (Organic Surface
Protection) coating. Moreover, note that the experimental work was carried
out in order to validate the analytical model of IMC growth but not of its
formation. The thickness of the IMC layers was measured before the heat
treatment and used as an initial condition for the theoretical model.
In order to reduce the number of tests and still get a reliable statistical
data, the experimental procedure to determine the kinetics of the IMC phase
growth was organized as follows. Before the experiment each of two types of
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 73
microchips was dissected through an array of balls, as shown schematically in
Fig. 4.3, (a). The cross-section was polished with 3
π
m diamond polishing
suspension. In order to track the interface propagation, two micro-indenter
marks (Vickers method) were applied to each ball by a Buehler MicroMet
5103 microindentation hardness tester. These marks served as reference for
tracking the IMC growth interfaces in the neighborhood, Fig. 4.3, (b). After
the heat treatment and the EDXS analysis, the specimens were returned to
the oven. The same cross-sections of the same solder balls were examined
during the experiment. Note that only three balls per series were continuously
monitored during the heat treatment, which was already quite time-consuming
and operationally elaborate. In the curves below showing thickness of IMC
over time an average of the observed growth for the three balls per series is
displayed.
Cu
Sn
(a) (b)
Fig. 4.3: Cutting scheme applied to the microprocessor.
The experimental procedure does not allow to fix the specimen in an epoxy
casing. Indeed, the melting temperature of the plastic is higher than melting
temperature of the solder material. However, the epoxy curing temperature is
close to the experimental temperature of 150
o
C. Therefore, all polishing was
performed βby hand,β which resulted in a relatively rough polished surface.
Nevertheless, the quality is fair enough to obtain the chemical interface kinetics
data, which is the aim of this study. Figure 4.4 shows the initial cross-section
of a solder ball. Differences in contrast correspond to different chemical
composition.
4.2.2 Experimental procedure
In the current study the process of IMC growth is analyzed in microprocessors
with BGA packages during a high temperature storage test. According to
the JESD22-A103 specification this experiment was performed at 150
Β°
C in a
vacuum oven over 1100 hours. The Series I ball was examined after 0, 120, 240,
360, 680, and 1000 h, and the Series II ball at 0, 120, 240, 360, 480, 800, and
1120 h, respectively. The chemical compositions were obtained from EDXS
Section 4.2. Experiment overview
74
(a)
Cu
Sn
100πm(b)
Cu6Sn5
Cu
Sn
10πm
Fig. 4.4:
Cross-sections of a solder ball from Series I before heat treatment, general
view (a) and interface zone (b).
analysis performed on an Oxford Instruments INCA X-Max detector fitted on
a MIRA 3 (TESCAN) microscope at voltage of 20 kV. The thickness evaluation
was based on postprocessing of micrographs with Python scripts.
4.3 Experimental results
The results for the two series are qualitatively different even though the
corresponding samples were heated and cooled under the same conditions. It
is suspected that this is due to the different casings and the associated heat
conduction properties: the thickness of the microprocessor case for Series I
(plastic casing) was about 150
π
m, and for series II (metal casing) about 360
π
m. Due to differences in the case materials and their thicknesses for samples
of series I and II, the true cooling rate down to room temperature could differ.
Both series of specimens contain a plastic circuit board and BGA in their
assembly. However, specimens from Series II have a additional metal framed
microchip. Therefore, the cooling time for Series II specimens is longer than
for specimens from Series I. The formation of IMC phases is a function of two
competing processes, namely growth and dissolution, which in turn depend
on the cooling rate. This could lead to the formation of various phases in
microprocessors of various designs.
Fig. 4.4 (b) shows that an IMC layer forms at the copper-solder interface
right after the attachment of solder balls on the microprocessor board copper
substrates. A spectroscopic analysis showed that its composition corresponds
to the compound Cu
6
Sn
5
. Note that the shape of the Cu
6
Sn
5
layer has a
comb-like structure (also referred to as scallops in the literature, e.g., Kim
and Tu [1996]) for all of the studied samples. Its relief is governed by the
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 75
roughness of the copper substrate, and the initial average height in the samples
was in the range 1.9 - 3.3
π
m. In addition, the tetragonal crystal structure of
Sn will result in an anisotropy of the diffusivity of the Cu atoms (Han et al.
[2017]; Wang et al. [2020]). It was also shown that the relief depends on the
crystalline structure of the underlying copper pad (Sunwoo, Morris, and Lucey
[1992]; Sang, Du, and Ye [2009]; Suganuma [2003]), which can be fine-grained
polycrystalline or even single-crystal based.
4.3.1 Chemical composition of the compounds
An EDXS analysis was performed at each stage of the experiment in order
to obtain the composition of the materials observed through the microscope.
The set of examined points for specimens from Series I after 1000 h heating
is shown in Fig. 4.5 and the corresponding compositions are listed in Tab.
4.1. All results are given in at% (atom percent). One should note that for
the composition 54.5 at% Cu and 45.5 at% Sn corresponds to Cu
6
Sn
5
and
the composition 75 at% Cu and 25 at% Sn to Cu
3
Sn, respectively. From the
analyzed data it follows that during heat treatment (up to 1000 h) the chemical
composition (stoichiometric ratio) of the IMC remained constant and no other
intermetallic compounds (e.g., Cu10Sn3or Cu41Sn11) were detected.
.1
.2
.3.4
10πm
Fig. 4.5:
First set of points from the
interface region for EDXS analysis, a
ball from Series I at
π‘
= 1000 h. The
corresponding data is listed in Tab. 4.1.
Point Cu Sn
1 54.5 45.5
2 52.4 47.6
3 52.7 47.3
4 53.1 46.9
Average 53.2 46,8
Tab. 4.1:
Data obtained from the
EDXS analysis for the points in Fig. 4.5,
values are given in at%.
The set of points for the EDXS analysis of the solder material and the Cu
substrate of a specimen from Series I after 1000 h heating is shown in Fig. 4.6
and the results of the analysis are listed in Tab. 4.2. Points 1 and 2 in Fig. 4.6
are located in the copper substrate area. However, they contain 2.6 and 5.6
at% of Sn, respectively. Moreover, during separation of Sn atoms from Cu
6
Sn
5
Section 4.3. Experimental results
76
.1.2
.3
.4
.5
.6
10πm
Fig. 4.6:
Second set of points from the
interface region for the EDXS analysis,
a ball from Series I at
π‘
= 1000 h. The
corresponding data is listed in Table
4.2.
Point Cu Sn
1 97.4 2.6
2 94.4 5.6
3 0.0 100.0
4 1.0 99.0
5 9.5 90.5
6 11.8 88.2
Tab. 4.2:
Data obtained from the
EDXS analysis for the points in Fig. 4.6,
values are given in at%.
a reaction proceeds according to the following scheme:
Cu6Sn5β3 Sn β2 Cu3Sn.(4.4)
As a result of reaction
(4.4)
and according to the phase diagram in Fig. 4.1, a
new Cu
3
Sn IMC should then be formed at the Cu β Cu
6
Sn
5
interface. However,
in Series I this compound was not observed, possibly due to an extremely small
thickness of the Cu3Sn layer, which could not be detected.
The set of examined points for a specimen from Series II after 1120 hours
heating is shown in Fig. 4.7 and the results of the chemical composition analysis
are listed in Tab. 4.3.
It can be seen from Fig. 4.7 that, unlike Series I, there are no βislandsβ in
the IMC layer and the boundaries of the layers show less roughness (comb-like
structure). During heat treatment a composition of two IMC configurations
was temporarily observed. However, at the later stages of the experiment (after
360 h) the Cu
3
Sn phase vanished. The reason for this is unknown. In any
case this phase is not so easy to detect. To quote from Ross, Vuorinen, and
Paulasto-KrΓΆckel [2016]: βHowever, it is well known that there is typically
a thin layer of Cu
3
Sn present between Cu and Cu
6
Sn
5
, which, immediately
after reflow, is generally only measurable by transmission electron microscopy
(TEM), although it can be resolved also by optical microscopy. Based on
literature data, we estimate the thickness of the Cu
3
Sn layer to be about 0.1
πm.β
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 77
.1
.2
3.
.4
25πm
Fig. 4.7:
Set of points from the interξΏ
face region for the EDXS analysis a
ball from Series II at
π‘
= 1120 h. The
corresponding data is listed in Table
4.3.
Point Cu Sn
1 0.0 100.0
2 31.5 68.5
3 54.7 45.3
4 97.0 3.0
Tab. 4.3:
Data obtained from the
EDXS analysis for the points in Fig. 4.7,
values are given in at%.
4.3.2 Thickness data for the IMC layers
In order to determine the characteristics of the growth kinetics of the inter-
metallic compounds a set of experiments was conducted using Series I and II
BGAs when exposed to a constant temperature. For each series measurements
were carried out as follows.
The heat treatment was performed at a constant temperature of 150
Β°
C.
Measurements of the thickness were done at the same location of 3 balls, namely
near the micro-indenter mark. This procedure allows to track the thickness of
the layers at certain points and to average between the balls obtained within
each set of microprocessor BGAs.
As it was mentioned above, in the Series I balls only the Cu
6
Sn
5
phase was
detected. It is known from the literature that the formation of the Cu
6
Sn
5
phase is largely determined by the diffusion processes of Cu, and the diffusion
rate determines the thickness of the IMC layer (Schaefer, Fournelle, and Liang
[1998]; Lee and Mohamad [2013]). Fig. 4.8 shows typical micrographs of a
sample after 120 hours storage at a constant temperature of 150 Β°C.
By comparison of Figs. 4.4 (initial state) and 4.8 it follows that the thickness
of the Cu
6
Sn
5
layer significantly increased. The comb-structure of the IMC
phase is preserved. Note that in one set of samples the thickness of Cu
6
Sn
5
varies along the interface from 3.1 to 9.1
π
m (in comparison, the initial layer
thickness was 1.9-3.3
π
m, i.e., it increased 2 - 3 times). A large dispersion of
thickness, as well as a formation of IMC agglomerations, can be explained by
uneven grain growth rate (Yu and Wang [2008]). Fig. 4.8 shows that there are
Section 4.3. Experimental results
78
(a) 100πm(b)
Cu6Sn5
(96-98)Sn - (4-2)Cu
10πm
Fig. 4.8:
Cross-section of a solder bump from Series I after
π‘
= 120 h heat treatment,
general view (a) and enlarged region indicated by the red rectangle (b).
local border areas, βislandsβ, in which the initial composition of the solder is
preserved.
Intermetallic layers after 120 and 240 hours heat treatment are shown in
Fig. 4.9 ((a) and (b), respectively). The approximate shape of the interfaces
between the solder, IMC and the substrate is outlined by red color. One
should note that the chemical composition of the IMC remains unchanged and
the appearance of new crystalline phases was not detected. The intermetallic
compound grows in both directions, toward the solder and toward the substrate.
However, the growth rate of the latter is much less than the one on the solder
side. From Fig. 4.9 it is also seen that the area of the βislandβ regions decreases
with time, which indicates that the Cu
6
Sn
5
grains grow in all directions and
not only perpendicular to the copper interface.
(a) (b)
10πm 10πm
Fig. 4.9:
IMC profiles in the solder bump from Series I after heat treatment:
π‘
= 120 h
(a) and π‘= 240 h (b). Spheroids of Cu6Sn5are marked with green.
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 79
The dark βislandβ in the solder area (marked by green color in Fig 4.9) are
spheroids of Cu
6
Sn
5
that formed by spalling from the IMC-solder interface
(Liu et al. [1996]; Lee and Mohamad [2013]).
Both Sn and Cu
6
Sn
5
islands are independent spheroids. It is most likely, that
if more material is polished from the surface of the specimen, these spheroids
will vanish. Therefore, since these areas do not really belong to the interface,
they are not considered, i.e., neither added nor subtracted during further
evaluation of the IMC phase thickness.
With further heat treatment a consistent increase in the thickness of the IMC
layer is observed. Fig. 4.10 shows the interface profiles at treatment times 360,
680 and 1000 h. It can be seen that the formation of a locally thick intermetallic
layer occurs. When comparing Fig. 4.10 (a) and (c), it follows that the areas
of the solder βislandsβ slightly increased upon heat treatment. This could be
explainable by surface effects after to the initial cut so that material might
have moved out of plane and was then removed during polishing.
(a) (b) (c)
10πm 10πm 10πm
Fig. 4.10:
IMC profiles in a solder bump from Series I after heat treatment:
π‘
= 360 h
(a), π‘= 680 h (b) and π‘= 1000 h (c).
Microphotographs for the Series II specimen are shown in Fig. 4.11 for times
from 0 h up to 120 h. One should note that for this series the thickness of the
IMC layer was less planar than for Series I. The initial thickness of the Cu
6
Sn
5
phase varied from 1.7 to 4.0 πm.
(a) 20πm(b)
Sn
Cu6Sn5
Cu3Sn
Cu
10πm
Fig. 4.11:
Cross-section of a solder ball from Series II in the interface region before
treatment at π‘= 0, (a), and after treatment, π‘= 120 h, (b).
Section 4.3. Experimental results
80
As shown in Fig. 4.11 (b), two intermetallic phases formed after treatment:
Cu
6
Sn
5
and Cu
3
Sn. This is different from the result obtained for balls of
Series I. As it is known from the literature the formation of Cu
3
Sn occurs on
the interface between the copper substrate and the Cu
6
Sn
5
compound. The
following chemical reaction describes its formation (Hu and Ke [2014]):
Cu6Sn5+ 9 Cu β5 Cu3Sn.(4.5)
Despite storing in a vacuum oven, oxides and other compounds formed on the
surface of the specimens cross-section. This made it hard to evaluate the exact
thickness of the Cu
6
Sn
5
layer for specimens of Series II at time
π‘
= 480 h, since
part of the layer was covered with other materials (βoxidation dirtβ), see Fig.
4.12.
5πm
Fig. 4.12:
Cross-section of the interξΏ
face region in a ball from Series II afξΏ
ter 480 h heat treatment. The chemiξΏ
cal composition of the marked area obξΏ
tained by the EDXS analysis is given
in Tab. 4.4.
Element Average, at %
C20.57
O4.62
Si 0.33
Cl 0.47
Cu 70.62
Sn 3.39
Tab. 4.4:
Results of the EDXS analysis
of the marked area in Fig. 4.12.
The EDXS analysis shows the presence of 4
.
62% oxygen in the region where
only Cu is expected (see points 1 and 2 in Fig. 4.6, or point 4 in Fig. 4.7). To
remove the dirt cover an additional polishing was carried out. In the following
results, for this particular measurement, the thickness of the only visible IMC
layer is presented. However, this measurement is excluded during the analysis
for validation of the theoretical model. The thickness of the polished layer
can be estimated from Fig. 4.13, which shows the IMC interface profiles of
Series II microchips at the same spot at time
π‘
= 800 h (a) and
π‘
= 1120 h (b).
Therefore the thickness of the removed layer is not greater than 10 πm.
In order to get statistical data about the thickness of the IMC all micrographs
were postprocessed with Python scripts. The positions of one interface were
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 81
(a) (b)
20πm 20πm
Fig. 4.13:
Cross-section of the IMC profile in a solder bump from Series II after heat
treatment: π‘= 800 h (a) and π‘= 1120 h after additional polishing (b).
subtracted from the other along the specimen. The mean value was chosen as
the current thickness and the standard deviation was calculated in order to
obtain an error estimate.
Results for Cu
6
Sn
5
in Series I and II specimens are shown in Fig. 4.14. Due
to the uneven growth of the IMC layer and its comb-like shape, the standard
deviation from the mean value increases in time.
Time, π‘, [h]
Cu6Sn5thickness, β, [πm]
Fig. 4.14:
Growth of the mean thickness of Cu
6
Sn
5
phase for the specimen from Series
I and II. Experimental points are connected by dashed lines for visual clarity only.
The IMC type Cu
3
Sn was observed only in Series II specimens and even
then only for a limited time period. The change of its thickness is shown in
Section 4.3. Experimental results
82
Fig. 4.15. Thus, during annealing for more than 480 hours, the Cu
3
Sn phase
dissolved, but the mechanism of this process is not understood.
Time, π‘, [h]
Cu3Sn thickness, β, [πm]
Fig. 4.15:
Growth of the mean thickness of Cu
3
Sn for the specimens from Series I
and II. Experimental points are connected by dashed lines for visual clarity only.
4.4 Theoretical model
One of the aims of this work is to validate and quantify the model of reaction
front kinetics described in Section 2.1 on the basis of the experimental results
for the IMC growth.
4.4.1 Analytical solution of a model problem
As mentioned in the previous sections, only the Cu
6
Sn
5
intermetallic phase was
detected at all stages of the experiment. Therefore, for simplicity of the further
analytical analysis, only this IMC is considered. It was pointed out that the
interface between Cu
6
Sn
5
and Sn moves much faster than one between Cu
6
Sn
5
and Cu. Therefore, it is also assumed that only the interface between the
IMC and the Sn is moving. The maximal thickness of the IMC layer observed
in the experiment was about 15
π
m (Fig. 4.14) which is much less than the
diameter of the solder ball (approx. 500
π
m, see Fig. 4.4). The above gives a
reason to consider a simple boundary value problem of mechanochemistry with
a planar chemical reaction front propagating in an infinite layer. Similar plane
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 83
problems for moving interfaces have been previously studied analytically and
numerically, see, e.g., Freidin et al. [2016]; Morozov et al. [2018a]; Morozov
et al. [2018b] and also in the previous chapters of this manuscript.
Due to the experimental setup, the cross-sections of the specimens were
stress free surfaces during the whole experiment. The chemical reaction front
propagation was observed at these stress free surfaces. Strictly speaking, one
cannot use neither a plane strain nor a plane stress simplification in order
to find the stresses at the reaction front. Nevertheless, since one can find a
kinetic equation for the interface movement for both plane stress and plane
strain formulations in a closed form, the model problem is solved in this work
in both formulations, assuming that the real behavior may be somewhat βin-
between.β This gives the opportunity to fit the theoretical prediction with the
experimental results, as described in the next Section 4.4.2. Based on the fitted
data, an estimate is obtained for the diffusion coefficient and for the chemical
reaction kinetic constant.
Consider an elastic layer with a cross-section in the
π₯π¦
-coordinate plane and
the plane reaction front propagating in the
π¦
-direction from
π¦
= 0 to
π¦
=
π»
where
π»
is the layer thickness and
β
is the current reaction front position (Fig.
4.16). The storage temperature in the oven is
π
and the reference temperature
is π0,π=πβπ0.
Assume that, by boundary conditions, the upper side of the layer is traction-
free, then
πβ
π¦|π¦=π»= 0,(4.6)
and the lower side is fixed, then the displacement
u+|π¦=0
= 0. Assume that
the displacement is zero in π₯-direction: π’Β±
π₯= 0 and, thus,
πΒ±
π₯= 0.(4.7)
The displacement and traction continuity conditions at the reaction front
will give, in particular, the continuity of
π¦
-components of the displacement
vector and the stress tensor:
Jπ’π¦Kπ¦=β= 0,Jππ¦Kπ¦=β= 0.(4.8)
The plane strain or the plane stress conditions are
πΒ±
π§
= 0 or
πΒ±
π§
= 0,
respectively. The equilibrium equations and boundary and interface conditions
in both plane statements will be satisfied for zero non-diagonal components of
stress and strain tensors, that is why these components are not mentioned above.
The displacement continuity condition at the reaction front also demands
continuity of
ππ§
, which strictly speaking is not fulfilled in the plane stress
formulation, but the input of this formal incompatibility is neglected.
Section 4.4. Theoretical model
84
100πm10πm
Cu6Sn5
Sn
Cu
π₯
π¦
β
π»
0
-
β
+
β
Fig. 4.16: Model description.
By constitutive equations (2.7),
πΒ±
π₯βπΌΒ±πβπtr
Β±=1
πΈΒ±(οΈπΒ±
π₯βπΒ±(πΒ±
π¦+πΒ±
π§))οΈ,(4.9)
where
πΈΒ±
and
πΒ±
are the Youngβs moduli and Poissonβs ratios of the materials
π΅Β±
,
πtr
β
= 0
, πtr
+β‘πtr
. The formulas for
πΒ±
π¦
and
πΒ±
π§
follow from
(4.9)
by cyclic
permutations of π₯, π¦, and π§.
Substituting strains and stresses found from
(4.9)
by using the conditions
(4.6)
β
(4.8)
and by taking into account all constraints regarding plane strain
and plane stress formulations into
(2.14)
, one can obtain that the contribution
πin π΄ππ does not depend on the front position, and
π=β§
βͺ
βͺ
β¨
βͺ
βͺ
β©
πΎ+1
2πΈβπΌ2
βπ2β1
2πΈ+(πΌ+π+πtr)2,plane stress,
πΎ+πΈβ
1βπβ
πΌ2
βπ2βπΈ+
1βπ+
(πΌ+π+πtr)2,plane strain.(4.10)
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 85
Then, by (2.5), the reaction rate becomes equal to
ππ=π*(πβπ
π*),(4.11)
where the stoichiometric coefficient
π*
is renormalized to one with respect to
all other stoichiometric coefficients, as shown in Chapter 2.1, and
π
is defined
as
π
= exp (οΈβπβπβ
πβ
π
π
π )οΈ.(4.12)
Diffusing Cu atoms are supplied through the lower side
π¦
= 0 and the
formation of the intermetallic phase occurs at the interface between Cu
6
Sn
5
and Sn at π¦=β. The diffusion equation takes the form
d2π
dy2= 0, π¦ β[0, β](4.13)
with the boundary conditions
π·dπ
dyξξξξπ¦=β
+ππ= 0,
π·dπ
dyξξξξπ¦=0
+π(π*βπ(0)) = 0,
(4.14)
where ππis given by (2.5).
The solution of the diffusion problem finally gives the concentration at at
the reaction front as a function of the reaction front positon
π|π¦=β=π*
1 + π
(οΈπ*
π+π*β
π·)οΈ
1 + π*
π+π*β
π·
.(4.15)
Then the reaction rate takes a form of the dependence of the front position
ππ=π*π*(1 βπ
)
1 + π*
π+π*β
π·
.(4.16)
Substitution of the expression
(4.16)
of the reaction rate into the formula
(2.6)
for the interface velocity results in the explicit equation for the dependence of
the front position on time:
dβ
dπ‘=πβπβ
πβ
π*π*(1 βπ
)
1 + π*
π+π*β
π·
.(4.17)
Section 4.4. Theoretical model
86
Integration yields an explicit dependence of time on the front position:
π‘(β) =
1
2
π*
π·(β2ββ2
0) + (οΈ1 + π*
π)οΈ(βββ0)
πβπβ
πβ
π*π*(1 βπ
)
,(4.18)
where
β0
is the initial thickness of the IMC layer appeared just after the solder
bump attachment.
The dependence
(4.18)
be rewritten in the standard form of the parabolic
law
β(π‘) = βοΈπΆ1π‘+πΆ2+πΆ3,(4.19)
where
πΆ1= 2π·πβπβ
πβ
π*(1 βπ
),
πΆ2=(οΈβ0βπ·(οΈ1
π*
+1
π)οΈ)οΈ2
,
πΆ3=βπ·(οΈ1
π*
+1
π)οΈ.
(4.20)
Note that the parameters
π*
and
π
occur in
(4.20)
only in a combination of
the sum of the inverse values. Note also that the influence of mechanical
actions is represented by the parameter
π
, which is constant in the considered
case. External mechanical loading may lead to the dependence
π
on the front
position, which would affect the front behavior.
The parabolic law
(4.19)
will be used in the next section for fitting parameters
to experimental data.
4.4.2 Fitting the model parameters
The mechanical properties of Cu, Sn and corresponding IMCs can be found in
the literature. Youngβs modulus, Poissonβs ratio, and the coefficient of thermal
expansion are listed in Tab. 4.5 (Yang et al. [2008]; Jiang et al. [1997]). Molar
masses and densities are given in Tab. 4.6 (Sun and Yin [2009]).
Material πΈ, GPa π πΌ, Kβ1
Sn 50 0.36 22 Β·10β6
Cu6Sn5118 0.31 18 Β·10β6
Tab. 4.5:
Mechanical material properties used in the model, Yang et al. [2008]; Jiang
et al. [1997].
The solubility of Cu in the IMC is defined as the maximum achievable
concentration of copper in a Cu
6
Sn
5
lattice. According to the phase diagram in
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 87
Sn Cu Cu6Sn5Cu3Sn
π, g/cm37.280 8.960 8.280 9.140
π, g/mol 118.7 63.55 974.8 309.3
Tab. 4.6: Molar masses and densities used in the model, Sun and Yin [2009].
Fig. 4.1 two stable phases can appear in the IMC at
π
= 150
o
C: Cu
6
Sn
5
and
Cu
3
Sn. By the reaction
(4.5)
the phase Cu
3
Sn forms when one mole of Cu
6
Sn
5
reacts with nine moles of Cu. Therefore it is estimated that the maximum
amount of copper which can be dissolved in Cu
6
Sn
5
is nine moles of Cu per
one mole of the IMC. In addition, here it is assumed that diffusing atoms of
Cu do not change the volume of the IMC. Hence, with the aforementioned
assumptions and based on the general definition of molar concentration,
π*
=
9
πCu6Sn5/πCu6Sn5
= 76
Γ
10
β6
mol/mm
3
was computed for the reference
value.
As mentioned earlier, the chemical energy parameter
πΎ
in
(4.10)
is determined
by the Cu
6
Sn
5
formation energy. This energy is defined empirically. The
dependence of this energy on temperature in operating temperature range can
be approximated by (Huang et al. [2015])
ΞπΊ=β7747.65 β0.371π[J/mol],(4.21)
where
π
is taken in K. In order to calculate
πΎ
, the formation energy has to
be divided by the molar volume
πCu6Sn5
= 11
.
28
cm3/mol
(value from Sobiech
et al. [2011]), since πΎunit is an energy density per unit volume.
By
(2.11)
, for the reaction
(4.5)
between Sn and Cu the ratio of volumes one
can find
π½tr
= 1
.
44 (44% volume expansion) if
π
= 0 (solid skeleton approach)
and
π½tr
= 0
.
92 (8% volume shrinkage) if
π
=
β
1. The question about the value
of the transformation strain accompanying IMC formation remains open and
estimates of the relative volume change are in the range from
β
10%, (e.g.,
Mei, Sunwoo, and Morris [1992]; Lee and Lee [1998]; Bordere et al. [2018])
up to +44%, (e.g., Jadhav et al. [2010]; Chudnovsky [2017]). Further the
model parameters are fitted for
π
=
β
0
.
5, which corresponds to volume ratio
π½tr
= 1
.
18. In particular, such a value of deformation may be consistent with a
small strain approach used in the modeling.
Certain material and chemical parameters remain unknown, namely
π·
,
π*
,
π
. In order to get an estimate for these values one can try to fit the square root
curve, Eq.
(4.19)
, constraining the parameters such that physically acceptable
values for
π·
and
π*
will result. The fitting process was performed by using
the weighted least square method, where the weights of the points depended
on their estimated error.
It should be noted that the equations proposed in Section 2.1 for the chemical
Section 4.4. Theoretical model
88
reaction front kinetics represent the velocity of the interface in the initial,
undeformed configuration, while the micrographs track the position of the
interface in the current, deformed configuration. Therefore, before fitting, the
analytical curve and the experimental points were transferred to the initial
configuration, by multiplying by the corresponding transformation strain value.
The fitting curves for Series I and II specimens are shown in Fig. 4.17 and
Fig. 4.18, respectively.
Time, π‘, [h]
Cu6Sn5thickness, β, [πm]
Fig. 4.17: Fitted square root dependence (4.19) for Series I specimens.
As it was noted in Section 4.2, the specimens of Series II were covered with
oxides and other compounds at time
π‘
= 480 h so that a polishing procedure
was required. Therefore it is not surprising that the measured thickness of the
IMC layer at that data point departs from the general pattern, Fig. 4.18, and
contains an error that is hard to assess. Because of that the data point was
excluded from the fitting curve procedure.
Then the unknown material parameters
π·
,
π*
,
π
can be calculated from
(4.20)
. As it was noted, the kinetic parameter
π*
and the mass transfer
coefficient
π
appear only in combination, so they cannot be resolved uniquely.
Found values of π·and the combination (1/π*+ 1/π)are shown in Tab. 4.7.
In Yuan et al. [2015] diffusion coefficients of Cu in the IMCs were calculated
based on the measured composition profiles of the diffusion zones within the
temperature range of 130
o
C β 200
o
C. The authors showed that the diffusion
coefficient depends highly on temperature, stating that
π·
= 0
.
38
Γ
10
β17m2/s
at
π
= 130
o
C,
π·
= 9
.
5
Γ
10
β17m2/s
at
π
= 150
o
Cand
π·
= 60
Γ
10
β17m2/s
at
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 89
Time, π‘, [h]
Cu6Sn5thickness, β, [πm]
Fig. 4.18:
Fitted square root dependence
(4.19)
for Series II specimens. The point
marked with asterisk was excluded from the fitting procedure due to bad surface
condition of the specimen.
Parameter Formulation Series I Series II Average
π·, [m2/s] plane stress 2.1 1.5 1.8
Γ10β17 plane strain 6.4 4.8 5.6
1/π*+ 1/π, [s/m] plane stress 89 57 73
Γ10β7plane strain 28 18 23
Tab. 4.7: Estimated diffusion coefficient for two series of specimens.
π
= 170
o
C. Keeping this in mind, one can conclude that obtaining the result
of the same order of magnitude can be considered as a good agreement. A
comparison of the diffusion parameters, averaged from the two experimental
results, with other works is shown in Tab. 4.8.
[1] [2] [3] [4] This work
plane stress plane strain
π·[m2/s] 1.5 5.64 1.82 9.5 1.8 5.6
Γ10β17
Tab. 4.8:
Comparison of the obtained estimation of the diffusion coefficient with
literature data. [1] - Onishi and Fujibuchi [1975], [2] - Paul, Ghosh, and Boettinger
[2011], [3] - Kumar, Handwerker, and Dayanada [2011], [4] - Yuan et al. [2015].
Section 4.4. Theoretical model
90
Both values of diffusion coefficients obtained in this work (basing on plane
stress and plane strain formulations) correlate with the data from various
papers and can be used as a reference for more complicated models based on
the chemical affinity tensor concept. These models can be extended since large
strains are involved and non-linear anisotropic materials as well.
4.5 Conclusions
A high-temperature storage test was carried out for two groups of microchips
with eutectic SnAg solder ball grid arrays. The specimens were cut and polished
before the heat treatment. The growth of the intermetallic phase was examined
by using the same set of solder balls through the entire experiment. The
proposed experimental procedure allowed to determine IMC growth kinetics
analyzing the small set of specimens. With the use of microindenter marks a
relative movement of the IMC interfaces and change of layer thickness were
obtained. This in turn lead to the evaluation and quantification of the kinetics
of intermetallic growth.
The growth of the Cu
6
Sn
5
intermetallic compound was modeled analytically
based on the chemical affinity tensor concept and by taking a temperature
dependence into account. The experimental results were used to determine the
kinetic parameters in the chemical affinity tensor model for the first time. The
simplest model of infinitely wide layers of linear elastic solids was analyzed
and a theoretical prediction of the growth kinetics was obtained. The influence
of mechanical stresses on IMC growth was taken into account. By comparison
of the experimental data with the theoretical model the values of the diffusion
coefficient and of the chemical reaction constant were estimated. The obtained
diffusion coefficients correlate with the results from works of other researchers.
The kinetic parameters of the model can now be used as reference values for
more general cases with the kinetic equation based on the chemical affinity
tensor.
In work Morozov et al. [2018b] the influence of the temperature regimes,
in particular temperature cycling, is studied numerically. The change in
temperature has a triplicate effect on the IMC growths kinetics: (i) as a source
of the shear load, (ii) through thermal stresses in the solder bump, and (iii)
through the chemical energy. Numerically it was shown that under all these
conditions the IMC growth might go non-uniformly. However, the analysis was
performed only with approximately estimated diffusion constants and reaction
parameters, because the experiment from this chapter was carried out much
later. An essential outlook for current work would be a recalculation of the old
results with the new data given in this chapter.
Chapter 4. Experimental and theoretical studies of Cu-Sn intermetallic phase growth
Numerical and analytical studies of the chemical reaction front kinetics in solids 91
5 Conclusions
In this work coupled problems of mechanochemistry were considered, namely
the stress-affected kinetics of a chemical reaction front propagation in solids.
A coupling of mechanical stresses, diffusion, and chemical reactions was com-
prehensively investigated: analytically, numerically, and experimentally. The
main results of the thesis are as follows:
(i)
A numerical procedure was developed and verified for simulating the
chemical reaction front propagation in elastic solids;
(ii)
The analytical procedure for a linear stability analysis of an equilibrium
phase interface was extended to the case of a chemical reaction front.
The stability of the propagating reaction front was studied numerically
for cases of stable and unstable thermodynamic equilibrium interface
positions;
(iii)
The competition between the global kinetics of the interface propagation
and the local kinetics of interface perturbations was demonstrated in the
case of an unstable equilibrium interface position;
(iv)
A high-temperature storage test was carried out for microchips with
eutectic SnAg solder ball grid arrays. Based on the experimental results,
diffusion and reaction kinetics parameters for the case of intermetallic
growth in the context of the chemical affinity tensor concept were obtained
for the first time.
The problem of chemical interface propagation was modeled based on the
chemical affinity tensor concept. This concept allows studying the influence
of mechanical stresses on the kinetics reaction front propagation based on
fundamental thermodynamic principles. A number of problems were solved
analytically and numerically. It was demonstrated that the stresses (which
arise from external loading, due to the chemical transformation strains, etc.)
could accelerate, retard, or even block the reaction front.
One of the main aims of the present work was developing a stability analysis
procedure for chemical reaction fronts. This was motivated by the following
reasons. The supply of the diffusing reactant governs the chemical reaction.
Therefore, the reaction front may be forced to propagate towards the unstable
equilibrium position. The growing instabilities may lead to plasticity and failure.
92
For the stability study a linearized perturbed boundary value problem was
formulated and solved analytically for two cases, namely planar and cylindrical
interfaces.
Even for these simple geometries the solution of the moving chemical in-
terface problem and especially the stability analysis might be complicated.
Therefore, a numerical procedure was developed using FEM for simulations of
the chemical reaction front propagation. The implementation of the procedure
was accomplished with the aid of the commercial FE software Abaqus. The
used method was verified by the analytical solutions for various problems of
chemical reaction front propagation. The proposed numerical procedure was
verified by comparison with the analytical predictions and also compared and
cross-verified with CutFEM and IGA-based procedures. It was shown that
the used numerical procedure can adequately reveal physical stability and
instability. The verified numerical procedure allowed to simulate the chemi-
cal reaction front kinetics when it approached stable or unstable equilibrium.
Numerical simulations allowed to observe the competition between global and
local kinetics of the chemical interface.
Many material, diffusion, and chemical reaction parameters have to be
defined in order to model the reaction front propagation, even if the linear
elastic materials and the simplest Fickβs diffusion are considered. Some of
the parameters (e.g., diffusion coefficients, the reaction rate constant and the
solubility of the diffusing reactant in the solid phase of the reaction product)
cannot be easily estimated. That is why one of the aims was to validate and
quantify the model of reaction front kinetics on the basis of experimental results.
To do this, a high-temperature storage test was carried out for microchips
with eutectic SnAg solder ball grid arrays and the kinetics of intermetallic
compound growth was evaluated. These experimental results were used then
to determine the kinetic parameters of the model. Such an evaluation in the
context of the chemical affinity tensor concept was done for the first time and
naturally required a comparison with the data given in other sources. The
obtained diffusion coefficients correlate with the results from works of other
researchers. Therefore, one can conclude that the chemical affinity tensor
approach can be used to model the reaction front kinetics and that the kinetic
parameters of the model can now be used as reference values for more general
cases.
Thus a model of the propagation of chemical reaction fronts in deformable
solids was formulated and numerically implemented. A comprehensive analysis
of the kinetics and stability of such fronts was performed for various conditions
and various types of mechanical loading. A practical basis has been created
for expanding this model and numerical procedures for subsequent studies. An
outlook for future research motivated by the obtained results may include:
Chapter 5. Conclusions
Numerical and analytical studies of the chemical reaction front kinetics in solids 93
(i)
In this work the correlation of the analytical predictions with numerical
results is shown and the analytical model is quantitatively verified with
the use of experimental results. To complete the cross-validation circle,
the numerical procedure might be utilized to simulate existing engineering
problems with the parameters obtained from the experiment.
(ii)
The proposed numerical procedure may be extended to the general case
of large deformations. This development is needed, e.g., for modeling the
reaction of silicon lithiation, which is accompanied with 300% volumetric
expansion due to the chemical transformation or the silicon oxidation
reaction with 100% volumetric expansion.
(iii)
In this work only stationary diffusion equation was considered. However,
if the characteristic time of the reaction is much less than the relative
time of the diffusion, then the chemical reaction is diffusion controlled.
For this case a dynamic diffusion equation should be utilized.
(iv)
The processes of initial accumulation of diffusing reactant before starting
the reaction and separation of the reaction front from the boundary of
the body require an additional study.
Numerical and analytical studies of the chemical reaction front kinetics in solids XI
List of Figures
1.1
Cross section of a solder bump prior (a) and after (b) current
stressing, reprinted from Chao et al. [2007]........... 2
1.2
SEM morphology of 250 nm a-Si film on Cu after 1 (a), and
30 cycles (b), reprinted from Kasavajjula, Wang, and Appleby
[2007] ............................... 2
2.1
A schematic representation of a localized chemical reaction in
solids. ............................... 10
2.2 Unperturbed and perturbed interfaces Ξ0and Ξ. . . . . . . . 17
2.3
A schematic representation of the planar chemical reaction front.
24
2.4 πΏπ
-lines for the stability analysis of the planar chemical reaction
front for stable (π)and unstable (π)sets of elastic parameters. 28
2.5
The contour plots of: (
π
)β
πΎ
for given external load in material
parameter space; (
π
)β
πΏ1
for given external load in material
parameter space at corresponding
πΎ
from figure (a). Black solid
lines in (b) represent the πΏ1= 0 level. ............. 29
2.6 A hollow cylinder undergoing a phase transformation. . . . . 31
2.7 Dependence of the driving force on the interface radius. . . . 32
2.8
Dependence of thermodynamically equilibrium interface radius
on the external load for (
π
)β solid and (
π
)β hollow cylinders.
Dashed line in (
π
)corresponds to the inner radius of the hollow
cylinder............................... 33
2.9 πΏπ
-lines for the stability analysis of the cylindrical phase in-
terface for: (
π
)- stable set of elastic parameters for the solid
cylinder, (
π
)- unstable set of elastic parameters for the solid
cylinder, (
π
)- stable set of elastic parameters for the hollow
cylinder with internal radius
π
= 0
.
1
π
,(
π
)- stable set of elastic
parameters for the hollow cylinder with internal radius
π
= 0
.
5
π
.
................................... 36
2.10 πΏπ
-lines for the stability analysis of the cylindrical chemical
reaction front for: (π)β stable set of elastic parameters for the
hollow cylinder with internal radius
π
= 0
.
1
π
,(
π
)β unstable set
of elastic parameters for the solid cylinder. . . . . . . . . . . 37
List of Figures
XII
2.11
Phase transition zones for the stable set of material parameters
from Tab. 2.1. Blue and red dots correspond to the strain state
for planar interface configuration. Dashed lines correspond to
the condition (2.86), solid lines β to condition (2.88). . . . . . 40
2.12
Phase transition zones for the stable set of material parameters
from Table 2.3. Blue and red dots correspond to the strain state
for cylindrical interface configuration. Dashed lines correspond
to the condition (2.86), solid lines β to condition (2.88). . . . 40
3.1 The flowchart for the numerical algorithm. . . . . . . . . . . . 45
3.2
Geometry described by NURBS curves which is directly used
for the numerical simulation (a). Classical FE model for the
numerical simulation (b). . . . . . . . . . . . . . . . . . . . . 46
3.3
Definition of the normal vector to the piecewise linear interface
ontheFEmesh.......................... 47
3.4
Automated adjustment of the interface position during its prop-
agation. .............................. 48
3.5
(a) β set of reference points on the interface for the output.
(b) β results of numerical simulations for stable (red curve) and
unstable (blue curve) configurations for the phase transformation
front propagation problem. . . . . . . . . . . . . . . . . . . . 49
3.6
Typical shape of the phase interface after losing stability near
the equilibrium radius during numerical simulation. . . . . . . 49
3.7
Shape of initially circular interface after the loss of stability on
theFEmesh. ........................... 50
3.8
Results of numerical simulations for stable (red curve) and
unstable (blue curve) configurations for the chemical reaction
frontpropagation. ........................ 51
3.9
Kinetics of the initially perturbed interface with different pre-
defined modes of instability for the unstable set of material
parameters. Top row - βbackward reactionβ, bottom row - direct
reaction. ............................. 52
3.10
The von Mises stress distribution for stable (a) and unstable (a)
configurations. .......................... 52
3.11
Unstable configuration of a solid cylinder with predefined shape
of the interface. Stress distributions are analyzed along radial
linesOAandOB.......................... 53
3.12
Distributions of von Mises stress along corresponding directions
53
3.13 Distributions of the hoop stress along corresponding directions 54
3.14
Distributions of the radial stress along corresponding directions
54
List of Figures
Numerical and analytical studies of the chemical reaction front kinetics in solids XIII
3.15
Results of numerical simulations by CutFEM and by the remesh-
ing procedures. Propagation of a planar chemical reaction front
for stable configuration is illustrated by the evolution of
π¦
-
coordinate of three points of the front. . . . . . . . . . . . . . 57
3.16
Results of numerical simulations by CutFEM and by the remesh-
ing procedures. Kinetics of the initially perturbed interface for
the unstable configuration is illustrated by four snapshots of the
reaction front at different times. . . . . . . . . . . . . . . . . 58
3.17
Cropped examples of the finite-element meshes used in CutFEM
and remeshing procedures are shown for the snapshot at
π‘
=
17600 inFig.3.16. ........................ 59
3.18
Propagation of the initially planar chemical reaction front in a
body under shear loading. The kinetics is illustrated by four
snapshots of the reaction front at different times. . . . . . . . 60
3.19
Propagation of the initially planar chemical reaction front in
a body under shear loading. The kinetics is illustrated by the
evolution of π¦-coordinate of three points of the front. . . . . . 61
3.20
Considered loading cases: biaxial stretching (a) and piecewise-
linear prescribed displacement (b). . . . . . . . . . . . . . . . 62
3.21
Intermediate and final configurations of a closed-curve reaction
front in a 2D body under biaxial stretching. . . . . . . . . . . 63
3.22
Intermediate and final configurations of a closed-curve reaction
front in a 2D body under complex loading. . . . . . . . . . . 64
3.23
Cropped examples of the finite-element meshes for final config-
uration of a of a closed-curve reaction front in a 2D body under
complex loading. (a) β CutFEM results, (b) β remeshing results.
64
3.24 Model for the finite layer with circular pore or inclusion. . . . 65
3.25
Kinetics of the initially planar chemical interface in the vicinity
of a cylindrical pore (a), a rigid inclusion (b), and in a uniform
layer(c)............................... 65
3.26
The von Mises stress distribution in the vicinity of a pore for
the reaction front position shown in the right inset. . . . . . . 66
3.27
The von Mises stress distribution in the vicinity of a solid
inclusion for the reaction front position shown in the right inset.
66
4.1 Cu-Sn phase diagram, reprinted from FΓΌrtauer et al. [2013]. . 70
4.2 Photos of the microprocessor boards. . . . . . . . . . . . . . . 72
4.3 Cutting scheme applied to the microprocessor. . . . . . . . . 73
4.4
Cross-sections of a solder ball from Series I before heat treatment,
general view (a) and interface zone (b). . . . . . . . . . . . . 74
List of Figures
XIV
4.5
First set of points from the interface region for EDXS analysis,
a ball from Series I at
π‘
= 1000 h. The corresponding data is
listedinTab.4.1.......................... 75
4.6
Second set of points from the interface region for the EDXS
analysis, a ball from Series I at
π‘
= 1000 h. The corresponding
data is listed in Table 4.2. . . . . . . . . . . . . . . . . . . . . 76
4.7
Set of points from the interface region for the EDXS analysis
a ball from Series II at
π‘
= 1120 h. The corresponding data is
listedinTable4.3. ........................ 77
4.8
Cross-section of a solder bump from Series I after
π‘
= 120 h heat
treatment, general view (a) and enlarged region indicated by
the red rectangle (b). . . . . . . . . . . . . . . . . . . . . . . . 78
4.9
IMC profiles in the solder bump from Series I after heat treat-
ment:
π‘
= 120 h (a) and
π‘
= 240 h (b). Spheroids of Cu
6
Sn
5
are
markedwithgreen......................... 78
4.10
IMC profiles in a solder bump from Series I after heat treatment:
π‘= 360 h (a), π‘= 680 h (b) and π‘= 1000 h (c). . . . . . . . . 79
4.11
Cross-section of a solder ball from Series II in the interface region
before treatment at
π‘
= 0, (a), and after treatment,
π‘
= 120 h,
(b).................................. 79
4.12
Cross-section of the interface region in a ball from Series II after
480 h heat treatment. The chemical composition of the marked
area obtained by the EDXS analysis is given in Tab. 4.4. . . 80
4.13
Cross-section of the IMC profile in a solder bump from Series
II after heat treatment:
π‘
= 800 h (a) and
π‘
= 1120 h after
additional polishing (b). . . . . . . . . . . . . . . . . . . . . . 81
4.14
Growth of the mean thickness of Cu
6
Sn
5
phase for the specimen
from Series I and II. Experimental points are connected by
dashed lines for visual clarity only. . . . . . . . . . . . . . . . 81
4.15
Growth of the mean thickness of Cu
3
Sn for the specimens from
Series I and II. Experimental points are connected by dashed
lines for visual clarity only. . . . . . . . . . . . . . . . . . . . 82
4.16 Model description. . . . . . . . . . . . . . . . . . . . . . . . . 84
4.17 Fitted square root dependence (4.19) for Series I specimens. . 88
4.18
Fitted square root dependence
(4.19)
for Series II specimens.
The point marked with asterisk was excluded from the fitting
procedure due to bad surface condition of the specimen. . . . 89
List of Figures
Numerical and analytical studies of the chemical reaction front kinetics in solids XV
List of Tables
2.1
Stable and unstable sets of material properties and parameters
used in analytical study and numerical simulation for the planar
interfaceproblem. ........................ 30
2.2
Diffusion parameters used in analytical study and numerical
simulation for the planar interface problem. . . . . . . . . . . 30
2.3
Material properties and parameters used in analytical and nu-
merical simulation. Material β+β refers to the outer, material
βββ to the inner region, respectively. . . . . . . . . . . . . . . 37
4.1
Data obtained from the EDXS analysis for the points in Fig. 4.5,
values are given in at%. . . . . . . . . . . . . . . . . . . . . . 75
4.2
Data obtained from the EDXS analysis for the points in Fig. 4.6,
values are given in at%. . . . . . . . . . . . . . . . . . . . . . 76
4.3
Data obtained from the EDXS analysis for the points in Fig. 4.7,
values are given in at%. . . . . . . . . . . . . . . . . . . . . . 77
4.4 Results of the EDXS analysis of the marked area in Fig. 4.12. 80
4.5
Mechanical material properties used in the model, Yang et al.
[2008]; Jiang et al. [1997]. .................... 86
4.6
Molar masses and densities used in the model, Sun and Yin
[2009]................................ 87
4.7 Estimated diffusion coefficient for two series of specimens. . . 89
4.8
Comparison of the obtained estimation of the diffusion coefficient
with literature data. [1] - Onishi and Fujibuchi [1975], [2] - Paul,
Ghosh, and Boettinger [2011], [3] - Kumar, Handwerker, and
Dayanada [2011], [4] - Yuan et al. [2015]. . . . . . . . . . . . 89
List of Tables
Numerical and analytical studies of the chemical reaction front kinetics in solids XVII
Bibliography
Abeyaratne, R. and Knowles, J. K. (2006). Evolution of Phase Transitions: A
Continuum Theory. Cambridge University Press.
Ahmad, Z. and Viswanathan, V. (2017). βStability of electrodeposition at solid-
solid interfaces and implications for metal anodesβ. In: Physical Review
Letters 119 (5), p. 056003.
doi
:
10.1103/PhysRevLett.119.056003
.
url
:
https : / / link . aps . org / doi / 10 . 1103 / PhysRevLett . 119 .
056003.
Antimonov, M. A., Cherkaev, A., and Freidin, A.B. (2016). βPhase transforma-
tions surfaces and exact energy lower boundsβ. In: International Journal
of Engineering Science 98. Special Issue Dedicated to Sergey Kanaun
Β΄
s
70th Birthday, pp. 153β182.
Arafat, Y. et al. (2020). βA Model for Intermetallic Growth in Thin Sn Joints
Between Cu Substrates: Application to Solder Microjointsβ. In: Journal
of Electronic Materials, pp. 1β16.
Baggetto, L., Danilov, D., and Notten, P. H. L. (2011). βHoneycomb-structured
silicon: remarkable morphological changes induced by electrochemical
(de)lithiationβ. In: Advanced Materials 23, pp. 1563β1566.
Ball, J. M. and James, R. D. (1987). βFine phase mixtures as minimizers of
energyβ. In: Archive for Rational Mechanics and Analysis 100, pp. 13β52.
Barvosa-Carter, W. and Aziz, M. J. (2004). βInterfacial roughening during
solid phase epitaxy: Interaction of dopant, stress, and anisotropy effectsβ.
In: Journal of Applied Physics 96.10, pp. 5462β5468.
Barvosa-Carter, W. et al. (1998). βKinetically Driven Growth Instability in
Stressed Solidsβ. In: Physical Review Letters 81.7, pp. 1445β1448.
Bordere, S. et al. (2018). βUnderstanding of Void Formation in Cu/Sn-Sn/Cu
System During Transient Liquid Phase Bonding Process Through Dif-
fusion Modelingβ. In: Metallurgical and Materials Transactions B 49,
pp. 3343β3356.
Bourderau, S., Brousse, T., and Schleich, D. M. (1999). βAmorphous silicon
as a possible anode material for Li-ion batteriesβ. In: Journal of Power
Sources 81β82, pp. 233β236.
XVIII
Bowen, R. M. (1967). βToward a thermodynamics and mechanics of mixturesβ.
In: Archive for Rational Mechanics and Analysis 24.5, pp. 370β403.
Bower, A. F. and Guduru, P. R. (2012). βA simple finite element model
of diffusion, finite deformation, plasticity and fracture in lithium ion
insertion electrode materialsβ. In: Modelling and Simulation in Materials
Science and Engineering 20.4, p. 045004.
Brassart, L. and Suo, Z. (2013). βReactive flow in solidsβ. In: Journal of the
Mechanics and Physics of Solids 61.1, pp. 61β77.
Burman, E. and Hansbo, P. (Apr. 2012). βFictitious domain finite element
methods using cut elements: II. A stabilized Nitsche methodβ. In: Applied
Numerical Mathematics 62.4, pp. 328β341.
doi
:
10.1016/j.apnum.
2011.01.008.
Burman, E. et al. (Jan. 2018). βShape optimization using the cut finite element
methodβ. In: Computer Methods In Applied Mechanics and Engineering
328, pp. 242β261. doi:10.1016/j.cma.2017.09.005.
β (June 2019a). βCut topology optimization for linear elasticity with cou-
pling to parametric nondesign domain regionsβ. In: Computer Methods
In Applied Mechanics and Engineering 350, pp. 462β479.
doi
:
10.1016/
j.cma.2019.03.016.
β
(2019b). βHybridized CutFEM for elliptic interface problemsβ. In: SIAM
Journal On Scientific Computing 41.5, A3354βA3380.
doi
:
10.1137/
18M1223836.
Buttner, C. C. and Zacharias, M. (2006). βRetarded oxidation of Si nanowiresβ.
In: Applied Physics Letters 89.
Callister Jr., W. D. and Rethwisch, D (2010). Material Science and engineering.
An introduction. 8th ed. Wiley.
Chang, S., Moon, J., and Cho, M. (2015). βStress-diffusion coupled multiscale
analysis of Si anode for Li-ion batteryβ. In: Journal of Mechanical Science
and Technology 29(11), pp. 4807β4816.
Chao, B. et al. (2007). βInvestigation of diffusion and electromigration pa-
rameters for CuβSn intermetallic compounds in Pb-free solders using
simulated annealingβ. In: Acta Materialia 55.8, pp. 2805β2814.
Chao, B. H.-L. et al. (2009). βRecent advances on kinetic analysis of electromi-
gration enhanced intermetallic growth and damage formation in Pb-free
solder jointsβ. In: Microelectronics Reliability 49.3, pp. 253β263.
Bibliography
Numerical and analytical studies of the chemical reaction front kinetics in solids XIX
Chenchiah, I. V. and Bhattacharya, K. (2008). βThe relaxation of two-well ener-
gies with possibility unequal moduliβ. In: Archive for Rational Mechanics
and Analysis 187.3, pp. 409β479.
Chudnovsky, B. H. (2017). Transmission, Distribution, and Renewable Energy
Generation Power Equipment. 2nd ed. CRC Press.
Cogan, S. F. et al. (1984). βDiffusion in the CuSn binary system: application
to Nb
3
Sn compositesβ. In: Journal of Materials Science 19, pp. 497β500.
Cui, Z., Gao, F., and Qu, J. (2012a). βA finite deformation stress-dependent
chemical potential and its applications to lithium ion batteriesβ. In:
Journal of the Mechanics and Physics of Solids 60, pp. 1280β1295.
β
(2012b). βInterface-reaction controlled diffusion in binary solids with
applications to lithiation of silicon in lithium-ion batteriesβ. In: Journal
of the Mechanics and Physics of Solids 61, pp. 293β310.
Dariavach, N. et al. (2006). βIntermetallic Growth Kinetics for Sn-Ag, Sn-Cu,
and Sn-Ag-Cu Lead-Free Solders on Cu, Ni, and Fe-42Ni Substratesβ.
In: Journal of Electronic Materials 35.7, pp. 1581β1592.
Duddu, R. et al. (Feb. 2011). βDiffusional evolution of precipitates in elastic
media using the extended finite element and the level set methodsβ. In:
Journal of Computational Physics 230.4, pp. 1249β1264.
doi
:
10.1016/
j.jcp.2010.11.002.
Dybkov, V. I. (2010). Reaction Diffusion and Solid State Chemical Kinetics:
Handbook. Trans Tech Publications Ltd.
Eremeev, V. A., Freidin, A. B., and Sharipova, L. L. (2003). βNonuniqueness
and Stability in Problems of Equilibrium of Elastic Two-Phase Bodiesβ.
In: Doklady Physics 48 (7), pp. 359β363.
β
(2007). βThe stability of the equilibrium of two phase elastic solidsβ. In:
Journal of Applied Mathematics and Mechanics 71, pp. 61β84.
Freidin, A. B. (1989). βCrazes and shear bands in glassy polymer as layers of a
new phaseβ. In: Mechanics of Composite Materials 1, pp. 1β7.
β
(1999). βSmall strains approach in the theory of strain induced phase
transformationsβ. In: Strenth and Fracture of materials (Ed. N.F. MoroξΏ
zov), Studies on Elasticity and Plasticity (St. Petersburg University) 18,
266β290, (in Russian).
β
(2007). βOn new phase inclusions in elastic solidsβ. In: ZAMM Journal
of Applied Mathematics and Mechanics 81.2, pp. 102β116.
Bibliography
XX
Freidin, A. B. (2009). βOn chemical reaction fronts in nonlinear elastic solidsβ.
In: Proceedings of XXXVII International Summer School-Conference
Advanced Problems in Mechanics. Ed. by Indeitsev, D.A. and Krivtsov,
A. M., pp. 231β237.
β
(2010). Fracture mechanics: Eshelby problem (in Russian). St.Petersburg
State Polytechnic University, St.Petersburg.
β
(2013). βChemical affinity tensor and stress-assist chemical reactions
front propagation in solidsβ. In: Proceedings of the ASME 2013 InternaξΏ
tional Mechanical Engineering Congress and Exposition. Vol. 9. American
Society of Mechanical Engineers, V009T10A102.
β
(2015). βOn the chemical affinity tensor for chemical reactions in de-
formable materialsβ. In: Mechanics of Solids 50.3, pp. 260β285.
Freidin, A. B. and Chiskis, A.M. (1994a). βPhase transition zones in nonlin-
ear elastic isotropic materials. Part 1: Basic relationsβ. In: Izv. RAN,
Mekhanika Tverdogo Tela (Mechanics of Solids) 29.4, pp. 91β109.
β
(1994b). βPhase transition zones in nonlinear elastic isotropic materials.
Part 2: Incompressible materials with a potential depending on one
of deformation invariantsβ. In: Izv. RAN, Mekhanika Tverdogo Tela
(Mechanics of Solids) 29.5, pp. 46β58.
Freidin, A. B. and Sharipova, L. L. (2018). βForbidden Strains and Stresses
in Mechanochemistry of Chemical Reaction Frontsβ. In: Generalized
Models and Non-classical Approaches in Complex Materials 1. Advanced
Structured Materials. Ed. by Altenbach, H. et al. Vol. 89. Springer Cham.
Chap. 17, pp. 335β348.
β
(2019). βTwo-phase equilibrium microstructures against optimal com-
posite microstructuresβ. In: Archive of Applied Mechanics 89, pp. 561β
580.
Freidin, A. B. and Vilchevskaya, E. N. (2020). βChemical Affinity Tensor in
Coupled Problems of Mechanochemistryβ. In: Encyclopedia of Continuum
Mechanics. Ed. by Altenbach, H. and Γchsner, A. Berlin, Heidelberg:
Springer Berlin Heidelberg.
Freidin, A. B., Vilchevskaya, E. N., and Korolev, I.K. (2014). βStress-assist
chemical reactions front propagation in deformable solidsβ. In: InternaξΏ
tional Journal of Engineering Science 83, pp. 57β75.
Freidin, A. B. et al. (2006). βSpherically symmetric two-phase deformations
and phase transition zonesβ. In: International Journal of Solids and
Structures 43(14), pp. 4484β4508.
Bibliography
Numerical and analytical studies of the chemical reaction front kinetics in solids XXI
Freidin, A. B. et al. (2015). βChemical reactions in spherically-symmetric
problems of mechanochemistryβ. In: Acta Mechanica 227(1), pp. 43β56.
Freidin, A. B. et al. (2016). βChemical affinity tensor and chemical reaction
front propagation: Theory and FE-simulationsβ. In: International journal
of Fracture 202, pp. 245β259.
Fried, E. (1993). βStability of a two-phase process in an elastic solidβ. In:
Journal of Elasticity 31.3, pp. 163β187.
Fu, Y. B. and Freidin, A. B. (2004). βCharacterization and stability of two-
phase piecewise-homogeneous deformationsβ. In: The Royal Society 46,
pp. 3065β3094.
FΓΌrtauer, S. et al. (2013). βThe CuβSn phase diagram, Part I: New experimental
resultsβ. In: Intermetallics 34, pp. 142 β147.
Gao, F. et al. (2006). βMicrostructure and mechanical properties evolution of
intermetallics between Cu and Sn-3.5Ag solder doped by Ni-Co additivesβ.
In: Journal of Electronic Materials 35, pp. 905β911.
Gao, H. et al. (2019). βEffect of nickel (Ni) on the growth rate of Cu
6
Sn
5
intermetallic compounds between SnβCuβBi solder and Cu substrateβ.
In: Journal of Materials Science: Materials in Electronics 30, pp. 2186β
2191.
Gibbs, J. (1948). The Collected Works of J.W. Gibbs, Vol. 1: Thermodynamics.
Yale University Press.
Glansdorff, P. and Prigogine, I. (1971). Thermodynamic theory of stability and
fluctuation. Wiley-interscience, New-York.
Grabovsky, Y. and Truskinovsky, L. (2011). βRoughening Instability of Broken
Extremalsβ. In: Archive for Rational Mechanics and Analysis 200 (1),
pp. 183β202.
β (2013). βMarginal Material Stabilityβ. In: Journal of Nonlinear Science
23, pp. 891β969.
Grinfeld, M. A. (1980). βOn conditions of thermodynamic equilibrium of phases
of a nonlinearly elastic materialβ. In: Soviet Mathematics. Doklady. 251,
pp. 824β827.
β
(1982). βStability of heterogeneous equilibrium in system containing solid
elastic phasesβ. In: Doklady Akademii Nauk SSSR 253.6, pp. 836β840.
β
(1990). Methods of continuum mechanics in the theory of phase transforξΏ
mations. in Russian. Nauka, Moscow.
Bibliography
XXII
Grinfeld, M. A. (1991). Thermodynamic methods in the theory of heterogeneous
systems. Longman Sc & Tech.
Gross, D., Mueller, R., and Kolling, S. (2002). βConfigurational forces - mor-
phology evolution and finite elementsβ. In: Mechanics Research CommuξΏ
nications 29, pp. 529β536.
Gurtin, M. E. (1983). βTwo-phase deformations of elastic solidsβ. In: Archive
for Rational Mechanics and Analysis 84.1, pp. 1β29.
β
(2000). Configurational forces as basic concepts of continuum physics.
Springer, New York, Berlin, Heidelberg.
Han, J. et al. (2017). βEffects of Grain Orientation on Cu
6
Sn
5
Growth Behavior
in Cu
6
Sn
5
-Reinforced Composite Solder Joints During Electromigrationβ.
In: Journal of Electronic Materials 47.2, pp. 1705β1712.
Hansbo, P. (2005). βNitscheβs method for interface problems in computa-tional
mechanicsβ. In: GAMM-Mitteilungen 28.2, pp. 183β206.
Hansbo, P., Larson, M. G., and Larsson, K. (2017). βCut finite element methods
for linear elasticity problemsβ. In: Geometrically Unfitted Finite Element
Methods and Applications. Ed. by Bordas, S. P. A. et al. Springer Inter-
national Publishing, pp. 25β63.
Heidemeyer, H. et al. (2000). βSelf-limiting and pattern dependent oxidation
of silicon dots fabricated on silicon-on-insulator materialβ. In: Journal of
Applied Physics 87.9, pp. 4580β4585.
Hu, X. and Ke, Z. (2014). βGrowth behavior of interfacial CuβSn intermetallic
compounds of Cu/Sn reaction couples during dip soldering and agingβ. In:
Journal of Materials Science: Materials in Electronics 25.2, pp. 936β945.
Huang, L. et al. (2015). βThermodynamic understanding of Sn whisker growth
on the Cu surface in Cu(top)-Sn(bottom) bilayer system upon room
temperature agingβ. In: Journal of Applied Physics 117.21, p. 215308.
Hughes, T. J. R., Cottrell, J. A., and Bazilevs, Y. (2005). βIsogeometric analysis:
CAD, finite elements, NURBS, exact geometry and mesh refinementβ. In:
Computer Methods in Applied Mechanics and Engineering 194, pp. 4135β
4195.
HΓΌter, C. et al. (2017). βElectrode-electrolyte interface stability in solid state
electrolyte systems: Influence of coating thickness under varying residual
stressesβ. In: AIMS Materials Science 4 (4), pp. 867β877.
Bibliography
Numerical and analytical studies of the chemical reaction front kinetics in solids XXIII
Jadhav, N. et al. (2010). βUnderstanding the Correlation Between Intermetal-
lic Growth, Stress Evolution, and Sn Whisker Nucleationβ. In: IEEE
Transactions on Electronics Packaging Manufacturing 33.3, pp. 183β192.
James, R. D. (1981). βFinite deformation by mechanical twinningβ. In: Archive
for Rational Mechanics and Analysis 77.2, pp. 143β176.
Jia, Z. and Li, T. (2015). βStress-modulated driving force for lithiation reaction
in hollow nano-anodesβ. In: Journal of Power Sources 275, pp. 866β876.
Jiang, Nan et al. (1997). βThermal Expansion of Several Sn-Based Intermetallic
Compoundsβ. In: Scripta Materialia 37.12, pp. 1851 β1854.
Juntunen, M. and Stenberg, R. (2009). βNitscheβs Method for General Boundary
Conditionsβ. In: Mathematics of Computation 78.267, pp. 1353β1374.
Kao, D. et al. (1988). βTwo-dimensional thermal oxidation of silicon. II. Mod-
eling stress effects in wet oxidesβ. In: IEEE Transactions on Electron
Devices 35.1, pp. 25β37.
Kasavajjula, U., Wang, C., and Appleby, J. A. (2007). βNano- and bulk-silicon-
based insertion anodes for lithium-ion secondary cellsβ. In: Journal of
Power Sources 163, pp. 1003β1039.
Kim, H. K. and Tu, K. N. (1996). βKinetic analysis of the soldering reaction
between eutectic SnPb alloy and Cu accompanied by ripeningβ. In:
Physical Review B 53 (23), pp. 16027β16034.
Kim, K. S., Huh, S. H., and Suganuma, K. (2003). βEffects of intermetallic
compounds on properties of SnβAgβAu lead-free soldered jointsβ. In:
Journal of Alloys and Compounds 352.1, pp. 226β236.
Knyazeva, A. G. (2003). βCross Effects in Solid Media with Diffusionβ. In:
Journal of Applied Mechanics and Technical Physics 44.3, pp. 373β384.
Kubanov, L. and Freidin, A. B. (1988). βSolid phase seeds in a deformable
materialβ. In: Journal of Applied Mathematics and Mechanics 52, pp. 382β
389.
Kumar, S., Handwerker, C. A., and Dayanada, M. A. (2011). βIntrinsic and
Interdiffusion in CuβSn Systemβ. In: Journal of Phase Equilibria and
Diffusion 32, pp. 309β319.
Kunin, I. (1983). Elastic media with microstructure. Springer, Berlin.
Lee, B.-Z. and Lee, D. N. (1998). βSpontaneous growth mechanism of tin
whiskersβ. In: Acta Materialia 46, pp. 3701β3714.
Bibliography
XXIV
Lee, L. M. and Mohamad, A. A. (2013). βInterfacial Reaction of SnβAg-Cu
Lead-Free Solder Alloy on Cu: A Reviewβ. In: Advances in Materials
Science and Engineering 2013, pp. 1β11.
Levitas, V. I. and Hamed, A. (2013). βAnisotropic Compositional Expansion
and Chemical Potential for Amorphous Lithiated Silicon under Stress
Tensorβ. In: Scientific Reports 3:1615.
Liashenko, O., Gusak, A., and Hodaj, F. (2014). βPhase growth competition in
solid/liquid reactions between copper or Cu3Sn compound and liquid tin-
based solderβ. In: Journal of Materials Science: Materials in Electronics
25.10, pp. 4664β4672.
Lin, E. J. et al. (2017). βEffect of Cu solubility on electromigration in Sn(Cu)
micro jointβ. In: Journal of Applied Physics 122.9, p. 095702.
Liu, A. A. et al. (1996). βSpalling of Cu6Sn5 spheroids in the soldering reaction
of eutectic SnPb on Cr/Cu/Au thin filmsβ. In: Journal of Applied Physics
80.5, pp. 2774β2780.
Liu, X. H. et al. (2012a). βIn situ atomic-scale imaging of electrochemical
lithiation in siliconβ. In: Nature Nanotechnology 7, pp. 749β756.
Liu, X. H. et al. (2012b). βSize-Dependent Fracture of Silicon Nanoparticles
During Lithiationβ. In: ACSNANO 6.2, pp. 1522β1531.
Liu, X. H. et al. (2013). βSelf-Limiting Lithiation in Silicon Nanowiresβ. In:
ACSNANO 7.2, pp. 1495β1503.
Loeffel, K. and Anand, L. (2011). βA chemo-thermo-mechanically coupled the-
ory for elasticβviscoplastic deformation, diffusion, and volumetric swelling
due to a chemical reactionβ. In: International Journal of Plasticity 27.9,
pp. 1409 β1431.
Marcus R. B.and Sheng, T. T. (1982). βThe Oxidation of Shaped Silicon
Surfacesβ. In: Journal of The Electrochemical Society 129.6, pp. 1278β
1282.
Maugin, G. A. (2010). Configurational forces. Thermomechanics, physics,
mathematics, and numerics. Chapman & Hall/CRC.
McDowell, Matthew T. et al. (2013). βIn Situ TEM of Two-Phase Lithiation of
Amorphous Silicon Nanospheresβ. In: Nano letters 13, pp. 758β764.
Mei, Z., Sunwoo, A. J., and Morris, J. W. (1992). βAnalysis of low-temperature
intermetallic growth in copper-tin diffusion couplesβ. In: Materials SciξΏ
ence and Engineering: A 23, pp. 857β864.
Bibliography
Numerical and analytical studies of the chemical reaction front kinetics in solids XXV
Mihalyi, A., Jaccodine, R. J., and Delph, T. J. (1999). βStress effects in the
oxidation of planar silicon substratesβ. In: Applied Physics Letters 74.14,
pp. 1981β1983.
Min-Suk, S., Chan-Jin, P., and Hyuk-Sang, K. (2008). βGrowth kinetics of
CuβSn intermetallic compounds at the interface of a Cu substrate and
42Snβ58Bi electrodeposits, and the influence of the intermetallic com-
pounds on the shear resistance of solder jointsβ. In: Materials Chemistry
and Physics 110.1, pp. 95β99.
Moghadam, M. M. and Voorhees, P. W. (2016). βThin film phase transformation
kinetics: From theory to experimentβ. In: Scripta Materialia 124, pp. 164
β168.
Monroe, C. and Newman, J. (2004). βThe Effect of Interfacial Deformation on
Electrodeposition Kineticsβ. In: Journal of The Electrochemical Society
151.6, A880βA886.
Morganti, S. et al. (2015). βPatient-specific isogeometric structural analysis of
aortic valve closureβ. In: Computer Methods in Applied Mechanics and
Engineering 284, pp. 508β520.
Morozov, A. V., Freidin, A. B., and MΓΌller, W. H. (2019). βStability of chemical
reaction fronts in the vicinity of a blocking stateβ. In: PNRPU Mechanics
Bulletin 2019.3, pp. 58β64.
Morozov, A. V. et al. (2020). βExperimental and Theoretical Studies of Cu-
Sn Intermetallic Phase Growth During High-Temperature Storage of
Eutectic SnAg Interconnectsβ. In: Journal of Electronic Materials.
Morozov, N. F. and Freidin, A. B. (1998). βZones of Phase Transitions and
Phase Transformations in Elastic Bodies under Various Stress Statesβ.
In: Proceedings of the Steklov Institute of Mathematics 223, pp. 219β232.
Mueller, R. and Gross, D. (1998). β3D simulation of equilibrium morphologies
of precipitatesβ. In: Computational Materials Science 11, pp. 35β44.
β
(1999). β3D inhomogeneous, misfitting second phase particles equilibrium
shapes and morphological developmentβ. In: Computational Materials
Science 16, pp. 53β60.
Mueller, R., Gross, D., and Lupascu, D. C. (2006). βDriving forces on do-
main walls in ferroelectric materials and interaction with defectsβ. In:
Computational Materials Science 35, pp. 42β52.
Muhlstein, C., Brown, S., and Ritchie, R. (2002). βHigh-cycle fatigue and
durability of polycrystalline silicon thin films in ambient airβ. In: Sensors
and Actuators A94, pp. 177β188.
Bibliography
XXVI
Muhlstein, C. and Ritchie, R. (2003). βHigh-cycle fatigue of micronscale poly-
crystalline silicon films: fracture mechanics analyses of the role of the
silica/silicon interfaceβ. In: International Journal of Fracture 119/120,
pp. 449β474.
Muhlstein, C., Stach, E., and Ritchie, R. (2001). βA reaction-layer mechanism
for the delayed failure of micron-scale polycrystalline silicon structural
films subjected to high-cycle fatigue loadingβ. In: Acta Materialia 50,
pp. 3579β3595.
MΓΌller, W. H., Vilchevskaya, E. N., and Freidin, A. B. (2015). βStructural
changes in micro-materials: Phenomenology, theory, applications, and
simulationsβ. In: Lecture Notes of TICMI 16, pp. 3β72.
Natsiavas, P. P. et al. (2016). βEffect of prestress on the stability of elec-
trodeβelectrolyte interfaces during charging in lithium batteriesβ. In:
Journal of the Mechanics and Physics of Solids 95, pp. 92 β111.
Onishi, M. and Fujibuchi, M. (1975). βReaction-Diffusion in the CuβSn Systemβ.
In: Transactions of the Japan Institute of Metals 16, pp. 539β547.
Ortiz, M., Repetto, E. A., and Si, H. (1999). βA continuum model of kinetic
roughening and coarsening in thin filmsβ. In: Journal of the Mechanics
and Physics of Solids 47.4, pp. 697β730.
Osmolovsky, V. G. (2000). The variational problem of phase transitions in
continuum mechanics. in Russian. Publ. of SPb university, St.Petersburg.
Paul, A. (2004). βThe Kirkendall effect in solid state diffusionβ. English. PhD
thesis. Department of Chemical Engineering and Chemistry.
isbn
: 90-
386-2646-0.
Paul, A., Ghosh, C., and Boettinger, W. J. (2011). βDiffusion Parameters and
Growth Mechanism of Phases in the CuβSn Systemβ. In: Metallurgical
and Materials Transactions A 42A, pp. 952β963.
Phan, A.-V. et al. (2001). βOn Transient Layers as New Phase Domains in
Composite Materialsβ. In: Modelling and simulation in materials science
and engineering 9.1, pp. 309β325.
Poluektov, M. and Figiel, Ε. (2019). βA numerical method for finite-strain
mechanochemistry with localised chemical reactions treated using a
Nitsche approachβ. In: Computational Mechanics 63.5, pp. 885β911.
Prigogine, I. and Defay, R. (1954). Chemical thermodynamics. Longmans,
Green.
Bibliography
Numerical and analytical studies of the chemical reaction front kinetics in solids XXVII
Ross, G., Vuorinen, V., and Paulasto-KrΓΆckel, M. (2016). βVoid formation and
its impact on CuSn intermetallic compound formationβ. In: Journal of
Alloys and Compounds 677, pp. 127 β138.
Rusanov, A. I. (2005). βSurface thermodynamics revisitedβ. In: Surface Science
Reports 58.5-8, pp. 111β239.
β
(2006). Thermodynamic foundations of mechanochemistry. in Russian.
Nauka, St. Petersburg.
Sang, X., Du, K., and Ye, H. (2009). βAn ordered structure of Cu
3
Sn in Cu-Sn
alloy investigated by transmission electron microscopyβ. In: Journal of
Alloys and Compounds 469 (1), pp. 129β136.
Schaefer, M., Fournelle, R. A., and Liang, J. (1998). βTheory for intermetal-
lic phase growth between Cu and liquid Sn-Pb solder based on grain
boundary diffusion controlβ. In: Journal of Electronic Materials 27.11,
pp. 1167β1176.
Schneider, D. et al. (Aug. 2017). βOn the stress calculation within phase-
field approaches: A model for finite deformationsβ. In: Computational
Mechanics 60.2, pp. 203β217.
Schneider, D. et al. (Mar. 2018). βSmall strain multiphase-field model ac-
counting for configurational forces and mechanical jump conditionsβ. In:
Computational Mechanics 61.3, pp. 277β295.
Silhavy, M. (1997). The Mechanics and Thermodynamics of Continuous Media.
Springer-Verlag Berlin, Heidelberg.
Sobiech, M. et al. (2011). βPhase formation at the Sn/Cu interface during room
temperature aging: Microstructural evolution, whiskering, and interface
thermodynamicsβ. In: Journal of Materials Research 26.12, pp. 1482β
1493.
Suganuma, K. (2003). Lead-Free Soldering in Electronics: Science, Technology,
and Environmental Impact. New York: CRC Press, Marcel Dekker.
Sun, F. and Yin, Z. (2009). βThe interfacial Cu-Sn intermetallic compounds
(IMCs) growth behavior of Cu/Sn/Cu sandwich structure via induc-
tion heating methodβ. In: Journal of Materials Science: Materials in
Electronics 30, pp. 18878β18884.
Sunwoo, A. J., Morris, J. W., and Lucey, G. K. (1992). βThe growth of Cu-Sn
intermetallics at a pretinned copper-solder interfaceβ. In: Metallurgical
Transactions A 23 (4), pp. 1323β1332.
Bibliography
XXVIII
Svendsen, B., Shanthraj, P., and Raabe, D. (Mar. 2018). βFinite-deformation
phase-field chemomechanics for multiphase, multicomponent solidsβ. In:
Journal of the Mechanics and Physics of Solids 112, pp. 619β636.
Truesdell, C. (1969). Rational Thermodynamics. McGraw-Hill, London.
Vilchevskaya, E. N., Filippov, R. A., and Freidin, A. B. (2013). βOn Transient
Layers as New Phase Domains in Composite Materialsβ. In: Mechanics
of Solids 48.1, pp. 92β118.
Vilchevskaya, E. N. and Freidin, A. B. (2007). βOn Phase Transitions in a
Domain of Material Inhomogeneity. I. Phase Transitions of an Inclusions
in a Homogeneous External Fieldβ. In: Mechanics of Solids 42.5, pp. 823β
840.
β
(2013). βOn kinetics of chemical reaction fronts in elastic solidsβ. In:
Surface Effects in Solid Mechanics. Ed. by Altenbach, H. and Morozov,
N. F. Springer Berlin Heidelberg, pp. 181β194.
Wang, Y. et al. (2020). βEffect of Sn Grain c-Axis on Cu Atomic Motion in
Cu Reinforced Composite Solder Joints Under Electromigrationβ. In:
Journal of Electronic Materials 49, pp. 2159β2163.
Weinberg, K., Werner, M., and Anders, D. (Feb. 2018). βA chemo-mechanical
model of diffusion in reactive systemsβ. In: Entropy 20.2, p. 140.
Xin, M. et al. (2003). βDevelopment of CuβSn intermetallic compound at
Pb-free solder/Cu joint interfaceβ. In: Materials Letters 57.22, pp. 3361β
3365.
Yang, F. (2010). βEffect of local solid reaction on diffusion-induced stressβ. In:
Journal of Applied Physics 107.
Yang, P.-F. et al. (2008). βNanoindentation identifications of mechanical prop-
erties of Cu
6
Sn
5
, Cu
3
Sn, and Ni
3
Sn
4
intermetallic compounds derived
by diffusion couplesβ. In: Materials Science and Engineering: A 485,
pp. 305β310.
Yu, D. Q. and Wang, L. (2008). βThe growth and roughness evolution of
intermetallic compounds of SnβAgβCu/Cu interface during soldering
reactionβ. In: Journal of Alloys and Compounds 458.1, pp. 542β547.
Yu, D. Q. et al. (2005). βIntermetallic compounds growth between Snβ3.5Ag
lead-free solder and Cu substrate by dipping methodβ. In: Journal of
Alloys and Compounds 392.1, pp. 192β199.
Bibliography
Numerical and analytical studies of the chemical reaction front kinetics in solids XXIX
Yuan, Y. et al. (2015). βInvestigation of diffusion behavior in CuβSn solid state
diffusion couplesβ. In: Journal of Alloys and Compounds 661, pp. 282β
293.
Zeeshan, A. and Venkatasubramanian, V. (2017). βKinetically Driven Growth
Instability in Stressed Solidsβ. In: Physical Review Letters 119.5,
p. 056003.
Zenga, Z. et al. (2016). βIn situ measurement of lithiation-induced stress in
silicon nanoparticles using micro-Raman spectroscopyβ. In: Nano Energy
22, pp. 105β110.
Zhao, K. et al. (2012a). βConcurrent Reaction and Plasticity During Initial
Lithiation of Crystalline Silicon in Lithium-Ion Batteriesβ. In: ECS
Meeting Abstracts.
Zhao, K. et al. (2012b). βFracture and debonding in lithium-ion batteries with
electrodes of hollow core shell nanostructuresβ. In: Journal of Power
Sources 218, pp. 6β14.
Zhao, X. J., Bordas, S. P. A., and Qu, J. M. (Dec. 2013). βA hybrid smoothed
extended finite element/level set method for modeling equilibrium shapes
of nano-inhomogeneitiesβ. In: Computational Mechanics 52.6, pp. 1417β
1428.
β
(Aug. 2015). βEquilibrium morphology of misfit particles in elastically
stressed solids under chemo-mechanical equilibrium conditionsβ. In: JourξΏ
nal of the Mechanics and Physics of Solids 81, pp. 1β21.
Zhao, X. J. et al. (June 2013). βEffects of elastic strain energy and interfacial
stress on the equilibrium morphology of misfit particles in heterogeneous
solidsβ. In: Journal of the Mechanics and Physics of Solids 61.6, pp. 1433β
1445.
Bibliography
Numerical and analytical studies of the chemical reaction front kinetics in solids i
A Analytical solutions for planar and
cylindrical reaction fronts kinetics
One should note that analytical solutions based on the chemical affinity tensor
concept for various simple geometries were obtained in many works, e.g. for
planar interface in Freidin, Vilchevskaya, and Korolev [2014]; Morozov et al.
[2018b], for cylindrical problems in Vilchevskaya and Freidin [2013]; Morozov et
al. [2018a], and for spherically symmetric problems in Freidin [2015]; Freidin et
al. [2015]. Here the solutions for planar and cylindrical problems are presented
in a form used as unperturbed solutions in the stability analyses given in
Sections 2.3 β 2.5. Also, these solutions were used as a reference to validate
the numerical results in Chapter 3.
A.1 Planar reaction front kinetics
As an example, a chemical reaction in the infinite layer of thickness
π»
is
considered, as illustrated in Fig. 2.3. The diffusing constituent is supplied
through the lower boundary. This gives the following boundary conditions for
equations (2.19) and (2.21):
u=0,at π¦= 0,
u=π’0eπ¦,at π¦=π»,
uΒ·eπ₯= 0,π:eπ₯eπ¦= 0,at π₯= 0 and π₯=πΏ,
π·nΒ· βπ+πΌ(πβπ*)=0,at π¦= 0,
nΒ· βπ= 0,at π₯= 0 and π₯=πΏ.
(A.1)
Solid constituents are assumed to be linear elastic and isotropic:
πβ=Cβ:πβ,π+=C+: (π+βπtr),(A.2)
CΒ±=πΒ±II + 2πΒ±4I,(A.3)
where πΒ±and πΒ±are LamΓ© parameters.
Due to the symmetry of the problem solutions for displacements and con-
ii
centration depend only on π¦coordinate. Hence, they can be found as
uΒ±=(οΈπ΄Β±π¦+π΅Β±)οΈeπ¦,
π=π΄ππ¦+π΅π.
Introducing the parameter
π
=
β/π»
, the constants can be found from the
boundary conditions (A.1) as:
π΄+=π’0/π»(πβ+ 2πβ) + π(π++π+)(1 βπ)
π(πβ+ 2πβ) + (1 βπ)(π++ 2π+),
π΄β=π’0/π»(π++ 2π+)βπ(π++π+)π
π(πβ+ 2πβ) + (1 βπ)(π++ 2π+),
π΅+= 0, π΅β=π»(οΈπ’0
π»βπ΄β)οΈ.
(A.4)
One should note that
πΒ±
=
π΄Β±eπ¦eπ¦
and stresses can be obtained from
(A.2)
.
These stresses and strains are substituted into the expression for π΄ππ
π΄π=πβπβ
πβ
ππ+π*π
π ln π
π*
,(A.5)
where
ππ
is defined by Eq.
(2.23)
and represents the mechanical contribution
to the chemical affinity tensor at the interface position π.
The reaction rate can be rewritten as
π=π*π(οΈ1βexp (οΈβπ΄ππ
π
π )οΈ)οΈ=
=π*π(οΈ1βexp (οΈβπβπβ
πβπ
π ππβln π
π*)οΈ)οΈ.
(A.6)
Now if the following notation is introduced,
ππ= exp (οΈβπβπβ
πβπ
π ππ)οΈ.(A.7)
a splitting of the exponent leads to
π=π*π(οΈ1βππ
π*
π)οΈ=π*(πβπππ*).(A.8)
With this expression for the reaction rate, boundary conditions for the diffusion
Chapter A. Analytical solutions for planar and cylindrical reaction fronts kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids iii
problem reads
π·dπ
dyξξξξπ¦=β
+π*(π(β)βπππ*)=0,
βπ·dπ
dyξξξξπ¦=0
βπ(π*βπ(0)) = 0.
(A.9)
Then, the constants for the solution can be found uniquely
π΄π=βπ*π*(1 βππ)
π·(οΈ1 + π*
π)οΈ+π*β
, π΅π=π*+π·
ππ΄π,(A.10)
and the reaction rate as a function of the current thickness has the form
π(β) = π*π*(1 βππ)π·
π·(οΈ1 + π*
π)οΈ+π*β
.
One should note that
ππ
is also a function of current thickness. Then according
to Eq. (2.6) the reaction front velocity reads:
ππ¦=dβ
dπ‘=πβπβ
πβ
π=πβπβ
πβ
π*π*(1 βππ)π·
π·(οΈ1 + π*
π)οΈ+π*β
.(A.11)
This equation can be integrated numerically with explicit Euler scheme with
given initial thickness β0and time increment πΏπ‘:
βπ+1 =βπ+ππ
π¦πΏπ‘, (A.12)
where
ππ
π¦
is the front velocity calculated at
ππ‘β
iteration, and
π
= 1
,
2
Β· Β· Β·
is the
iteration number.
A.2 Cylindrical reaction front kinetics
The plane strain problem for the hollow and solid cylinder undergoing chemical
reaction (Figure 2.6) is considered. The internal and external radii are
π
and
π
, respectively, the current position of the reaction front is a cylindrical surface
of the radius
π
. Diffusing reactant is supplied through the outer surface of the
cylinder. The transformation strain is assumed to be plane:
πtr =πtr(eπeπ+eπeπ),(A.13)
Section A.2. Cylindrical reaction front kinetics
iv
where
ππ
and
ππ
are the unit vectors of the cylindrical coordinates. Solid
constituents are assumed to be linear elastic and isotropic as in the previous
example, (A.2).
In order to find the stresses and strains at the reaction front we use the
equilibrium equation, which in the case of axial symmetry takes the form
dππ
dπ+ππβππ
π= 0 (A.14)
with boundary conditions for the hollow cylinder
ππ(π)=0, ππ(π)=0,Jπ’(π)K= 0,Jππ(π)K= 0 (A.15)
or for the solid cylinder
π’π|πβ0= 0, ππ(π)=0,Jπ’(π)K= 0,Jππ(π)K= 0.(A.16)
The radial displacement is given by the LamΓ© solution:
π’Β±(π) = π΄Β±π+π΅Β±
π,(A.17)
and the strains can be found through
πΒ±
π(π) = dπ’Β±(π)
dπ=π΄Β±βπ΅Β±
π2,πΒ±
π(π) = π’Β±(π)
π=π΄Β±+π΅Β±
π2.(A.18)
Four unknown constants
π΄Β±
and
π΅Β±
can be found from the boundary
conditions (A.15) or (A.16), and e.g. for the case of solid cylinder:
π΄+=
π’0
π(πβ+πβ+π+) + πtrπ2(π++π+)
π2(π++ 2π+) + (1 βπ2)(πβ+πβ+π+),
π΄β=1
π2(οΈπ’0
π+ (π2β1)π΄+)οΈ,
π΅+=π2(οΈπ’0
πβπ΄+)οΈ, π΅β= 0
(A.19)
where π=π/π.
The stresses and strains at the reaction can then be found from
(A.2)
and
(A.18)
, respectively. These stresses and strains are substituted into the
expression for π΄ππ
π΄π=πβπβ
πβ
ππ+π*π
π ln π
π*
,(A.20)
Chapter A. Analytical solutions for planar and cylindrical reaction fronts kinetics
Numerical and analytical studies of the chemical reaction front kinetics in solids v
where
ππ
is defined by Eq.
(2.23)
and represents the mechanical contribution
to the chemical affinity tensor at the interface position π.
As in the previous example in A.1, the reaction rate takes form
(A.8)
with
same notation for the
ππ
(which depends on the current interface radius) being
introduced.
The corresponding stationary diffusion equation in the case of axial symmetry
becomes: d
dπ(οΈ1
π
dπ
dπ)οΈ= 0, π β[π, π].(A.21)
Boundary and interface conditions according to the notation of Figure 2.6 can
be written as follows (as described in Section 2.1):
π·dπ
drξξξξπ=π
βπΌ(π*βπ(π)) = 0, π· dπ
drξξξξπ=π
βπ*(π(π)βπ*ππ)=0.(A.22)
In a cylindrical coordinate system the solution of the Laplace equation can
be found as follows:
π(π) = πΆ1ln (οΈπ
π)οΈ+πΆ2.(A.23)
The two unknown constants
πΆ1
and
πΆ2
can be found uniquely from the
boundary conditions (A.22):
πΆ1=π*(ππβ1)
ln π
πβπ·(οΈ1
π*π+1
πΌπ)οΈ, πΆ2=π*βπ·
πΌππΆ1.(A.24)
From (A.23), (A.24) it follows that
π(π) = π*
ln π
πβπ·(οΈ1
π*π+1
πΌπ)οΈ{οΈ(οΈln π
πβπ·
πΌπ)οΈππβπ·
π*π}οΈ(A.25)
Since n=βeπthe normal component of the reaction front velocity is
ππ=βdπ
dπ‘=πβπβ
πβ
π*(π(π)βπ*ππ),(A.26)
and the reaction front kinetics β the dependence
π
(
π‘
)β can be obtained by
integration of Eq.
(A.26)
where the dependencies
π
(
π
)and
ππ
are already
found.
Section A.2. Cylindrical reaction front kinetics
