- Email: [email protected]
Mathematical model of copper corrosion
Applied Mathematical Modelling 38 (2014) 4804–4816
Contents lists available at ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
Mathematical model of copper corrosion F. Clarelli a,⇑, B. De Filippo b, R. Natalini b a Istituto per le Applicazioni del Calcolo ‘‘M. Picone’’, Consiglio Nazionale delle Ricerche, Polo Scientifico – Edificio F CNR, Via Madonna del Piano 10, I-50019 Sesto Fiorentino (FI), Italy b Istituto per le Applicazioni del Calcolo ‘‘M. Picone’’, Consiglio Nazionale delle Ricerche, via dei Taurini 19, I-00185 Roma, Italy
a r t i c l e
i n f o
Article history: Received 29 November 2012 Received in revised form 19 December 2013 Accepted 25 March 2014 Available online 2 April 2014 Keywords: Free boundary model Parabolic problems Finite difference methods Corrosion Copper Brochantite
a b s t r a c t A new partial differential model for monitoring and detecting copper corrosion products (mainly brochantite and cuprite) is proposed to provide predictive tools suitable for describing the evolution of damage induced on bronze specimens by sulfur dioxide (SO2) pollution. This model is characterized by the movement of a double free boundary. Numerical simulations show a nice agreement with experimental result. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction Deterioration of copper and bronze artifacts is one of the main concerns for people working in cultural heritage [1]. More specifically, bronze, a copper-tin alloy, has been widely employed for daily-use and artistic purposes from the Bronze Age up to present. Conservation studies are based on a knowledge of the environmental conditions to which copper and copper alloys may be exposed and include all the information on the material technologies and the nature of the corrosion films or patina, which cover their surfaces. In particular a significant effort has been devoted to study the corrosion due to environmental conditions, such as temperature, moisture and the concentration of pollutants [2–4]. Although in recent years air pollution in European urban areas has decreased considerably, there still remain concentrations of pollutants such as sulfur dioxide (SO2) from combustion of fossil fuels, being one of the most important factors in the deterioration of bronze. Indeed SO2, mixed with water vapor, reacts to produce sulfuric acid (H2SO4), which causes corrosion phenomena on copper surfaces and produces several corrosion products as basic copper sulphates, such as antlerite, posnjakite, brochantite [5]. The latter is the final product of several reaction steps, which can be approximated by two main chemical reactions: cuprite formation, which occurs after a few weeks of exposure to atmospheric conditions, and brochantite formation, which is the final reaction step (for more details on chemical background see Section 2.1) [6]. The complexity of corrosion processes creates the necessity for a quantitative model approach to develop predictive tools, which simultaneously provide both quantitative information as well as simulations of the various processes involved. These methods, similar to those introduced in [7], are useful for the monitoring and detection of surface alterations even before they are visible, making it possible to determine optimal intervention strategies. In this paper we introduce a new partial ⇑ Corresponding author. Tel.: +39 055 522 5813. E-mail addresses: [email protected] (F. Clarelli), b.defi[email protected] (B. De Filippo), [email protected] (R. Natalini). http://dx.doi.org/10.1016/j.apm.2014.03.040 0307-904X/Ó 2014 Elsevier Inc. All rights reserved.
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
4805
differential model, which is used to describe the evolution of damage induced on a bronze specimen by atmospheric pollution. It is based on fluid dynamical and chemical relations and it is characterized by a double free boundary: one between copper and cuprite, the other between cuprite and brochantite. A similar approach was first proposed in the seminal paper [8], and then explored in many other works, see for instance [9–11]. In our paper we apply this approach to metal corrosion, and we validate it by calibrating the model according to the experimental results in [12]. The paper is organized as follows: in Section 2.1 we analyze the main chemical corrosion phenomena and in Section 2.2 we review the main mathematical models already proposed in literature. The Section 3 is entirely focused on the description of the model’s equations and the numerical schemes used, meanwhile in Section 4 we describe the experimental setting and the calibration of the model. Finally, in Section 5, we present the main results related to the simulations produced by our model and in Section 6 the related conclusions. 2. Modeling backgrounds 2.1. Chemical backgrounds Most of the oxidation processes occurring on bronze metal artifacts, when exposed to environmental conditions, are electrochemical and involve interactions between the metal surface, the adsorbed moisture and various atmospheric gases (SO2, CO2, NOx , hydrocarbons) [13]. Electrochemical corrosion processes in electrolyte and condensed moisture layers have been the subject of extensive studies, based on numerous and different approaches [14–16]. When exposed to the atmosphere, copper and its alloys form a thin layer of corrosion, from brownish-green up to greenish-blue colors, which is designated as patina. In the case of copper in a low pollutant levels atmosphere, a native cuprite (copper (I) oxide or Cu2O) film of approximately a few nanometers thick, protects the metal surface from further oxidation. The general reaction is well described in literature, where, in aerated solution, copper can dissolve electrochemically forming copper (I) oxide formation, due to the reaction of copper with oxygen. It is represented by the following schematic reaction [17]:
1 2Cu þ O2 ! Cu2 O: 2
ð1Þ
However, in an aggressive environment, like urban atmosphere, the protective nature of this oxide layer is altered and there is the formation of a non-protective, multi-component, tarnish layer. When the copper is exposed to humidity and sulfur dioxide, three kinds of basic copper sulfate hydroxide are mainly produced: Brochantite Cu4SO4(OH)6, which is a wellknown patina constituent, or other similar products like Antlerite or Posnjakite. In the following we focus our attention to the formation of Brochantite, which is the main observed product. Gaseous sulfur dioxide and sulfate particles are deposited on the electrolyte on the cuprite. Their deposition reduces the pH of the adsorbed water, and this promotes the dissolution of cuprous ion (Cu+) and its oxidation, thus forming cupric ions (Cu2+). In detail, copper (I) ions in solution disproportionate to give copper (II) ions and a precipitate of copper
2CuþðaqÞ ! Cu2þ ðaqÞ þ CuðsÞ : When the cupric and sulfate ion concentrations in the electrolyte are high enough to form brochantite, this phase starts to precipitate on the cuprite [18,19]. In [20] it is indicated that, in the initial oxidation process, cuprite formation is followed by posnjakite, as a precursor phase to brochantite, see also [21], but we are going to neglect this intermediate transformations, due to the elevated speed of the reaction with respect to the time scale of the mathematical model which we are going to present. 2.2. Existing mathematical models Graedel and his collaborators [22,23] have studied atmospheric copper sulfidation at AT&T Bell Laboratories in both experimental investigation and physical or mechanistic model development; furthermore a systematic investigation of copper sulfidation kinetics has been performed. In the paper [23] Tidblad and Graedel developed a model able to describe the SO2 copper corrosion. Their work is based on aqueous chemistry and without considering spatial dimensions. In the paper [24], they introduce a model which describes the corrosion of copper exposed to moist air with a low SO2 concentration. Here, they consider four stages in the development of a corrosion patina: the metal, a non-protective oxide film which has a high ionic transport property, an outer layer of corrosion products which can permit the penetration of water and gases and an external absorbed water layer. More recently, Larson proposed in [25] a different model that describes the atmospheric sulfidation of copper by proposing a physical copper-sulfidation model that includes four distinct phases: the substrate metal, a cohesive cuprous sulfide (Cu2S) product layer, a thin aqueous film adsorbed on the sulphide and the ambient gas. Larson postulated that transport through the sulfide layer occurs via diffusion and electromigration of copper vacancies and electron holes. Later, in the 1990s, copper sulfation by SO2 was investigated by Payer et al., who focused their attention on the early stage of corrosion in moist air (75% relative humidity (RH) at 25o C) with a sulfur dioxide concentration of 0:5%. The techniques employed (SEM, AES and TEM) allowed for an analysis and characterization of the oxide film (composed of
4806
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
cuprous oxide, copper sulfate and sulfide) on copper surfaces and for the mechanism of evolution of the corrosion chemistry to be described [2]. In the last decade the effect of sulfur dioxide in acid rain on copper and bronze has been investigated. Robbiola and his co-workers studied the effect of its cyclic action on bronze alloys, underlining the different types of patina formed in ’’sheltered’’ and ’’unsheltered’’ areas of bronze monuments [26,16]. Subsequently, Larson refined the copper-sulfidation model by focusing on the transport of charged lattice defects in a growing Cu2 product layer between the ambient gas and the substrate metal. As previously, this transport is postulated to occur via both diffusion and electromigration. 3. The model In this section, we aim to introduce a mathematical model able to describe the corrosion effects on a copper layer, which is subject to deposition of SO2. The present model is based on the mathematical approach used in [7]. We assume to have a copper sample on which is formed a non protective oxide layer (Cu2O), and, over this layer, a corrosion product (brochantite) grows. Over the brochantite layer is assumed to be the atmospheric air with SO2. An example of these three layers can be seen in Fig. 1. The reaction producing cuprite Cu2O, can be approximated by Eq. (1). Namely, two moles of copper combined with one-half mole of oxygen produce cuprite. Laboratory tests, with high concentration of SO2 and high relative humidity, near to 100%, show that brochantite is the primary product of the reaction, thus it can be considered as the final state reached by the whole reaction. The overall simplified reaction of brochantite formation can be approximate by the following reaction (2)
3 2Cu2 O þ SO2 þ 3H2 O þ O2 ! Cu4 ðOHÞ6 SO4 ; 2
ð2Þ
where two moles of cuprite combined with three moles of water and three-halves of oxygen produce one mole of brochantite. In the following, we assume that these two reactions are instantaneous; thus, we obtain a sharp free boundary between cuprite and the unreacted copper, due to the reaction (1), and a second free boundary between brochantite and cuprite, due to the reaction (2). Assuming these two reactions as instantaneous, the effective time of reaction is implicitly included in the diffusivity coefficients. 3.1. Swelling We define the boundary between copper and cuprite by aðtÞ, the boundary between cuprite and brochantite by bðtÞ and the boundary between brochantite and external air by cðtÞ. Summarizing, the geometry of our problem is one-dimensional, and we have 4 regions (e.g. see Fig. 2): 1. 2. 3. 4.
Copper (inner region). Cuprite Cu2O, between aðtÞ and bðtÞ. Brochantite Cu4 (OH)6 SO4, between bðtÞ and cðtÞ. Laboratory controlled atmosphere on the external side of cðtÞ.
Fig. 1. Example of cuprite and brochantite deposition on a copper sample.
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
4807
Fig. 2. Example of cuprite and brochantite theoretical growth (in time) on a copper sample.
When we speak about three dimensional quantities (e.g. volumes), we consider just one dimension of them, assuming a standard cross sectional area. It is assumed that, in the climatic chamber at a 100% of RH, a thin film of water (with SO2 dissolved) is formed over the external boundary, and it plays an important role in both reactions. The reaction (1) implies a production of Cu2O on the boundary between copper and the oxide layer. The transformation of copper into Cu2O proceeds if sufficient oxygen is available, and it is accompanied by a volume change (swelling rate). The swelling rate can be calculated easily, because the molar ratio in reaction between Cu and Cu2O is 2 : 1. Thus two moles of copper change into one mole of Cu2O, and a different volume of the new matter formed is obtained. The derivative with respect to the time of the swelling s1 of reaction (1) is
s_ 1 ¼ xp a_ ;
ð3Þ
where
xp ¼
lc 2lp
1;
ð4Þ
represents the variation of volume ratio, lc and lp being the molar density (moles/cm3) of Cu and Cu2O respectively. In the present case xp > 0, since it represents an expansion of the material. It is assumed that lc and that lp are constant (i.e. they are homogeneous materials). Under these assumptions, if s1 ð0Þ ¼ að0Þ ¼ 0, then it is possible to conclude that
s1 ðtÞ ¼ xp aðtÞ:
ð5Þ
On the boundary bðtÞ, between cuprite and brochantite, SO2 reacts with Cu2O in presence of water and oxygen, following the reaction (2). Here, 2 mol of cuprite change in 1 mol of brochantite, wasting 1 mol of SO2, 3 mol of H2O and 3=2 moles of oxygen. Thus, we have a cuprite production from the reaction (1) and a cuprite consumption from the reaction (2). Thus, the physical boundary between cuprite and brochantite bðtÞ is given by
bðtÞ ¼ bðtÞ þ s1 ðtÞ:
ð6Þ
Here, bðtÞ is the consumption of cuprite (assumed to be positive in our case) due to the reaction (2), and the derivative with respect to the time of the swelling rate of brochantite s2 is given by
_ s_ 2 ðtÞ ¼ xb bðtÞ;
ð7Þ
where
xb ¼
lp 2lb
1;
ð8Þ
and lb is the molar density of brochantite assumed constant. The variation of the overall external boundary c_ is given by
_ c_ ¼ xp a_ ðtÞ xb bðtÞ:
ð9Þ
3.2. Equation of the model It is well known that relative humidity plays a key role in regulating the speed of oxydation and sulfation. It has been observed that when RH exceeds some threshold, then SO2 reacts completely such as the oxidation of copper happens with
4808
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
full speed. We can interpret this phenomenon as follows. According to Eq. (2), when a molecule of SO2 comes in contact with cuprite, it reacts if three molecules of H2O are available at the same point (we suppose that there is always enough O2). To be more precise, this is true only if vapor condenses on the unreacted specimen surface forming a liquid film, and this happens with high RH values. The brochantite formation (Eq. (2)) has been assumed to develop on the cuprite layer, due to the SO2, H2O and O2, which move through the brochantite layer and react with Cu2O. This reaction implies a wasting of Cu2O, with the formation of a new layer of brochantite. Summarizing, the region of brochantite is given by cðtÞ 6 x 6 bðtÞ and the region of Cu2O is bðtÞ 6 x 6 aðtÞ. We define the concentration of SO2 in the pores of brochantite by S, the water concentration by W and the oxygen concentration by O throughout brochantite, and by G throughout cuprite. The flow of SO2 relative to air is governed by Fick’s law. Thus, in the frame of reference where copper is at rest, the SO2 flux has the following expression
@S @S J s ¼ nb Ds Sxp a_ Sxb b_ ¼ nb Ds þ Sc_ ; @x @x
ð10Þ
where the first term refers to SO2 diffusion in the brochantite layer, nb is the brochantite porosity, Ds is the diffusivity. The Sxp a_ and Sxb b_ terms refer to the swelling caused by cuprite and brochantite formation respectively. Hence the mass balance of SO2 in the brochantite layer cðtÞ 6 x 6 bðtÞ is
@S @2S @S Ds 2 þ c_ ¼ 0; @t @x @x
ð11Þ
The value of S at the external boundary cðtÞ is the environment SO2 concentration Sa ðtÞ, which is a known function of time
SðcðtÞÞ ¼ Sa ðtÞ;
ð12Þ
We assumed that SO2 reacts totally with Cu2O at the front bðtÞ, thus we have
SðbðtÞÞ ¼ 0:
ð13Þ
Now, we need a further condition (first free boundary). Since the SO2 flux of moles at the boundary bðtÞ is proportional to Cu2O moles consumption, we have
nb
Ds @S 1 qp _ b; ¼ M s @x 2 M p
ð14Þ
where M s ; M p are the molar weight of SO2 and Cu2O respectively, qp is the mass density of cuprite. Similarly, the water flux is
@W @W J w ¼ nb Dw W xp a_ W xb b_ ¼ nb Dw þ W c_ ; @x @x
ð15Þ
where Dw is the water diffusivity. Thus, the mass balance in the brochantite layer is
@W @2W @W Dw þ c_ ¼ 0: @t @x2 @x
ð16Þ
The value of W at the external boundary cðtÞ is the experimental water concentration W a ðtÞ
WðcðtÞÞ ¼ W a ðtÞ;
ð17Þ
Since some moles of water are wasted by the reaction (2), at the front bðtÞ we have
Jw 3 qp _ W _ b þ nb b: ¼ Mw Mw 2 Mp
ð18Þ
Finally, the oxygen flux is
@O J o ¼ nb Do þ Oc_ ; @x
ð19Þ
where Do is the oxygen diffusivity. Thus, the mass balance in the brochantite layer is
@O @2O @O Do 2 þ c_ ¼ 0: @t @x @x
ð20Þ
The value of O at the external boundary cðtÞ is the experimental oxygen concentration Oa ðtÞ
OðcðtÞÞ ¼ Oa ðtÞ: Since the oxygen is also wasted by the reaction (2), at the front bðtÞ we have
ð21Þ
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
Jo 3 qp _ O _ b þ nb b: ¼ Mo Mo 4 Mp
4809
ð22Þ
3.3. Cuprite layer equations We assumed that the reaction between oxygen and copper occurs on the copper-cuprous oxide boundary aðtÞ. Here, we indicate oxygen by G, just to avoid confusion with the previous region. We make a mass balance of oxygen concentration G, and we assume that all oxygen moles arriving on the inner boundary aðtÞ react. This assumption is not really true, in fact the water plays a key role in the speed of reaction, but in our experiments we have a RH near to 100%, and this fact justify our assumption; also, the diffusivity Dg includes implicitly the finite time of reaction. The flux of oxygen is
@G J g ¼ np Dg Gxp a_ @x
ð23Þ
and the mass balance equation is
@G @2G @G Dg 2 xp a_ ¼ 0: @t @x @x
ð24Þ
The value of G at the boundary bðtÞ is given by the value of oxygen on the boundary b, given by Eq. (22).
GðbðtÞÞ ¼ OðbðtÞÞ;
ð25Þ
Since oxygen reacts totally with copper at the boundary aðtÞ, we have
GðaðtÞÞ ¼ 0;
ð26Þ
Now, we need a further condition to close the system (second free boundary). Since the oxygen react totally at the boundary
aðtÞ, the copper moles wasted are given by np
Dg @G 1 qc ¼ a_ : M g @x 4 M c
ð27Þ
3.4. Rescaling The geometry of the problem is given by two regions, one is given by x 2 ½bðtÞ; aðtÞ which describes the layer of Cu2O, the other is given by x 2 ½cðtÞ; bðtÞ i.e. the brochantite layer. Also we have c_ ¼ xp a_ þ xb b_ . 3.4.1. New variables In the first region we adopt ðx; tÞ ! ðy; sÞ, in the second one ðx; tÞ ! ðz; sÞ. We find out
y¼
z¼
xb
;
@ x ¼ 1=ða bÞ@ y ;
@ xx ¼ 1=ða bÞ2 @ yy ;
xc ; bc
@ x ¼ 1=ðb cÞ@ z ;
@ xx ¼ 1=ðb cÞ2 @ zz ;
ab
@ t ¼ f ðtÞ@ y þ
@ t ¼ qðtÞ@ z þ
1 @s; tr
1 @s; tr
ð28aÞ
ð28bÞ
where
f ðy; sÞ ¼
1 yð@ s b @ s aÞ @ s b ; tr ab
ð29aÞ
qðz; sÞ ¼
1 zð@ s c @ s bÞ @ s c : tr bc
ð29bÞ
Now, we rescale the following parameters and variables in non-dimensional form:
^ ¼ W ; W Wr ^ g ¼ t r Dg ; D k2
^¼ G; G Gr
^S ¼ S ; Sr
^ w ¼ t r Dw ; D k2
a b b c a^ ¼ ; b^ ¼ ; b^ ¼ ; c^ ¼ ; k
^ s ¼ t r Ds ; D k2
k
^ðz; sÞ ¼ t r q; q
k
k
^f ðy; sÞ ¼ t r f
ð30aÞ
ð30bÞ
From now on, we assume that dotted variables are derivatives with respect to the dimensionless time s ðV_ ¼ dV=dsÞ. ^_ the non-dimensional system is: ^ s Mp Sr , Cw ¼ 3 1 Mw qp , Xg ¼ 4np D ^ g Mc Gr . Also, c ^_ ¼ ðxp a ^_ þ xb bÞ, Assuming Xs ¼ 2nb D 2 n Mp W r Ms q Mg q p
b
c
4810
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
3.4.2. Outer region For c 6 x 6 b ! 0 6 z 6 1
0 1 ^s ^_ D @ ^S @ 2 ^S @ c @ ^S ^A ; þq ¼ 2 2 @s @z ^c ^ ^c b ^ @z b
ð31Þ
^Sð0; sÞ ¼ ^Sa :
ð32Þ
^Sð1; sÞ ¼ 0;
ð33Þ
@ ^S ^_ ð1; sÞ ¼ b; ^c ^ @z b
Xs
ð34Þ
0 1 ^ ^w ^ ^ ^_ D @W @2W c @W ^A þq ¼ @ ; 2 2 @s @z @z ^c ^ ^c b ^ b
ð35Þ
^ ^ a ðsÞ: Wð0; sÞ ¼ W
ð36Þ
^ ^ D @W _ ^ 3 qp M w b; w ð1; sÞ ¼ c_ b_ W 2nb W r M p ^c ^ @z b
ð37Þ
0 1 ^ ^ ^ ^o ^_ D @O @2O c @O ^A þq ¼ @ ; 2 2 @s @z ^c ^ ^c b ^ @z b
ð38Þ
^ sÞ ¼ O ^ a ðsÞ: Oð0;
ð39Þ
^ ^ @O D _ ^ 3 qp M o b; o ð1; sÞ ¼ c_ b_ O 4nb Or M p @z ^ ^ bc
ð40Þ
3.4.3. Inner region For b 6 x 6 a ! 0 6 y 6 1,
! ^ ^ ^ ^g @G @2G a^_ @G D ^ ¼ þ xp f ðy; sÞ ; 2 @y2 ^ @ s ða @y ^ ^ b a ^ bÞ
ð41Þ
(we have to highlight that a_ is a derivative with respect to the s)
^ sÞ ¼ 0; Gð1;
ð42Þ
^ sÞ ¼ O ^ ; Gð0; b
^ Xw @ G ^ @y ^b a
^_ : ¼a
ð43Þ
ð44Þ
From now on we use the non-dimensional system, so we indicate all non-dimensional terms without hat. 3.5. Numerical solutions To solve our model, we have to set up an appropriate numerical scheme which is able to describe the process in a short range of time (about 40 h), but also simulations of 1 year taking under control the numerical stability. For these reasons we use finite differences schemes with implicit-explicit terms. It can be found in [27,28].
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
4811
3.5.1. Initial conditions Our numerical procedure requires að0Þ – 0 and bð0Þ – að0Þ to avoid singularities in the internal region confining with copper. By physical assumption, we know that að0Þ > 0 because it represents the first copper consumption. In order to initialize properly the external equations, where SO2 is present, we need to set up bð0Þ > 0 and cð0Þ – bð0Þ. 3.5.2. Numerical scheme Indicating explicit term as HðUÞ and implicit term by GðUÞ. We present our system in the external region (Eqs. 31, 35 and 38), which supplemented by boundary and transmission conditions (32)–(34), (36), (37), (39) and (40).
U t ¼ HðUÞ þ GðUÞ;
ð45Þ
where
0
S
1
B C U ¼ @ W A; O the explicit term HðUÞ is
1 _ ðbc cÞ þ q @S @z B C B C c_ @W C B HðUÞ ¼ B ðbcÞ þ q @z C; @ A _ ðbc cÞ þ q @O @z 0
and the implicit term GðUÞ is
0
Ds @2 S ðbcÞ2 @z2
B B Dw GðUÞ ¼ B ðb 2 @ cÞ
@2 W @z2
Do @2 O ðbcÞ2 @z2
1 C C C; A
since GðUÞ is a stiff term which will be integrated implicitly to avoid excessively small time steps. A more general scheme is IMEX-DIRK Runge–Kutta, which is given by, for t ¼ nDt,
uðiÞ ¼ un þ Dt
i1 i X X ~ik H uðkÞ þ Dt aik G uðkÞ ; a k¼1
unþ1 ¼ un þ Dt
i ¼ 1; . . . m
ð46Þ
k¼1
m X
~ i H uðiÞ þ Dt x
i¼1
m X
xi G uðiÞ :
ð47Þ
i¼1
e ¼ ða ~ik Þ, where a ~ik ¼ 0 for j P i and A ¼ ðaik Þ are m m matrices, x ~ i and xi are m-dimensional vectors, Here, the matrices A such that the resulting scheme is explicit in H and implicit in G. The DIRK formulation requires aik ¼ 0 for j > i [27]. Following the IMEX formalism we shall use the notation Name (s; m; p) to identify a scheme, where s is the number of stages of the implicit scheme, m is the number of explicit stages and p is the combined order of the scheme. In our case ~ 1; x ~ 2 Þ ¼ ð0; 1Þ and ðx1 ; x2 Þ ¼ ð0; 1Þ. we choose to use an Implicit-Explicit Midpoint (1; 2; 2): s ¼ 1; m ¼ 2 and p ¼ 2, where ðx By the same method, and at the same time-step of previous integration, we integrated also Eqs. (41)–(44). (See Table 1). 4. Assessment of the model 4.1. Experiments The calibration of the model has been made with reference to a precise experimental campaign, which is described in [12], where the main characteristics of the experimental setting are reported. Experiments were carried out in a cyclic corrosion cabinet (Erchisen Mod. 519/AUTO) and the Cu–12Sn cast bronze specimens were chosen as a representative alloy used in the past (ancient Greek type). For the corrosion test, all bronze specimens were exposed to an atmosphere containing about 200 ppm of SO2 at 40 C and 100% RH for 8 h (wet cycle), subsequently they were exposed to room conditions for 16 h (dry cycle). Each wet and dry cycle was repeated 20 times. The measures of patina thickness (see Fig. 3), useful for the model calibration, were performed using a SEM-EDS instrument equipped with an Image Analyzer (IA). Each measurement (after 8, 24 and 40 h) has been obtained by an average of 20 measures. The early stage of exposure of bronze samples mainly produced copper hydroxyl-sulfate (brochantite and chalcanthite), cuprous oxide, and tin sulfide (ottemannite). The tarnish product layer also contained a trace of tin oxide. As determined by XRD analysis, the brochantite formation was predominant with the increase of the number of corrosion cycles. Otherwise,
4812
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816 Table 1 Parameters table. Param.
Value
Dimension
Indications
qp
6:00
g/cm3
Cu2O mass density
Mp
143:09 8:94 63:55 3:97 452:3 1:46 64:07
Cu2O molar mass Copper mass density Copper molar mass Brochantite mass density Brochantite molar mass SO2 mass density SO2 molar mass Reference value
1 106
g/mol g/cm3 g/mol g/cm3 g/mol g/cm3 g/mol g/cm3
Gr
5:86 105
g/cm3
Reference value
Or
5:86 105
g/cm3
Reference value
Wr
2:3 105
g/cm3
Reference value
Oa
2:931 105
g/cm3
Reference value
tr
3:1536 107 15 0:3 0:002
sec
Time reference value
cm nondim. nondim.
Length reference value Porosity of brochantite Porosity of cuprite
qc Mc
qb Mb
qs Ms Sr
k nb np
the chalcanthite formation was detected during the early stage of the corrosion processes and the intensity of its characteristic peaks decreases as a function of time. The same consideration can be made for the ottemannite formation, detected with SEM-EDS in localized micro-crack areas, which are rapidly covered by basic copper sulfates [12,29]. 4.2. Calibration The laboratory corrosion tests have been used to calibrate the model. We suppose that the diffusivity coefficients are the parameters to calibrate by experimental tests. To do that, we used the thickness of corrosion products measured at three fixed times. Each thickness value has been obtained by a sample average and standard deviation, computed at each time, based on twenty random measurements, see Table 2. We used the least square method to find the best parameters. We obtained (in cm2/sec): Dg ¼ 9:9 109 , Ds ¼ 3:96 105 , Do ¼ 9:9 106 and Dw ¼ 3:96 105 . The evolution in time of the corrosion products thickness, using the best parameters, is in Fig. 4. We can see the difference between the experimental points and the best simulation in Fig. 5. Remark. here, we test the numerical scheme to reproduce the conservation of mass implicit in the chemical reactions described in Eqs. (1) and (2). Simulations, after 40 h, give us the following values: c ¼ 9:505 104 (cm), a ¼ 3:1693 104 (cm), b ¼ 7:9916 104 (cm). Knowing the molar volume of copper V c ¼ qc =M c , cuprite V p ¼ qp =M p , and brochantite V b ¼ qb =M b , we find that the number of copper moles wasted are a=V c ¼ 4:4588 105 , the number of cuprite moles formed are hp =V p ¼ 2:2294 105 , i.e. two copper moles are wasted to form one cuprite mole, as we expect from Eq. (1). Similarly, we obtain the number of cuprite moles wasted by reaction (2) is b=V p ¼ 2:2173 105 , while the number of brochantite moles formed by the same reaction is about 1:1086 105 , i.e. for each brochantite mole formed two cuprite moles are wasted, as expected from (2).
Fig. 3. Example of patina thickness measurement with the SEM employment. The green lines measure the thickness of corrosion products, mainly brochantite. From [12]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4813
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816 Table 2 Patina thickness measures. Time of measure (h)
Averaged value (cm)
8
1:7331 104
9:2672 104
1:8473 104
5:4418 10
24 40
Standard deviation
4
4
2:4102 104
13:2522 10
Fig. 4. Evolution of cðtÞ (black line), bðtÞ (blue line) and aðtÞ (red line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
−3
1.6
x 10
Experimental points Simulation
1.4
Thickness of corrosion products (cm)
1.2
1
0.8
0.6
0.4
0.2
0 −5
0
5
10
15
20
25
30
35
40
45
Time (hours) Fig. 5. Simulation (continuous line) and experimental points (in red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4814
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816 −11
4
x 10
3.5
3
SO mass density [g/cm ]
3 2.5
2
2 1.5 1 0.5 0 Jan
Feb Mar
Apr
May
Jun
Jul
Aug
Sep
Oct
Nov
Dec
Fig. 6. SO2 concentration detected at Villa Ada (Rome, Italy) during the 2006.
35 30 25
o
Celsius
20 15 10 5 0 −5 −10 Jan
Feb Mar
Apr
May
Jun
Jul Aug Time
Sep
Oct
Nov
Dec
Fig. 7. Temperature detected at Villa Ada (Rome, Italy) during the 2006.
−5
4
x 10
3.5
Saturated vapor density
SVD (g/cm3)
3 2.5 2 1.5 1 0.5 0 Jan
Feb
Mar
Apr
May
Jun
Jul
Aug Time
Sep
Oct
Nov
Dec
Fig. 8. Saturated vapor density obtained by temperature and RH detected at Villa Ada (Rome, Italy) in 2006.
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
4815
Fig. 9. Simulation of corrosion with environmental data detected at Villa Ada (Rome, Italy) in 12 months (2006).
5. Application In this section we use the calibrated model with environmental data, detected in Rome at Villa Ada during the year 2006. This application is just an example of how the model can be used. Since the model has been calibrated for high RH and SO2 concentration, the calibration should be improved by experiments with lower values of RH and SO2 concentration. The SO2 concentration detected is in Fig. 6 and the temperature is in Fig. 7. Assuming good calibration with our environmental data, we need the saturated vapor density (SVD in [g/cm3]) as function of RH and temperature T ½ C detected at Villa Ada. This is used to obtain an exact quantity of water vapor in the air from relative humidity, the density of water in the air is given by RH SVD ¼ ActualVaporDensity. The behavior of water vapor density is a non-linear function, but an approximate calculation of saturated vapor density can be made from an empirical fit of the vapor density curve in Eq. (48) (between 0–45°C is a good approximation).
SVDðTÞ ¼ 5:018 þ 0:32321T þ 8:1847 103 T 2 þ 3:1243 104 T 3 :
ð48Þ
The figure of saturated vapor density obtained is in Fig. 8. Thus, using SO2, temperature and humidity data detected at Villa Ada in Rome during the 2006, we obtained the behavior of corrosion products illustrated in Fig. 9. Results of this simulation represent a qualitative example of corrosion formation on copper under SO2 attack. They have to be improved by other laboratory experiments, because our model is calibrated using a very high SO2 concentration (never present in the atmospheric environment in this concentration), and under high RH (close to 100%). To have a more accurate calibration, it would be necessary set up mainly experiments with different SO2 concentration and different concentration of humidity in the air, so to calibrate the model under these variations. Our conditions produce a thick cuprite layer in 1-year simulations with respect to the laboratory experiments, because of the air SO2 concentration is lower than that present in the experimental room, and slowing the growth of brochantite layer. 6. Conclusions In conclusion a new mathematical model was developed to describe and simulate the evolution of brochantite formation on Cultural Heritage artifacts exposed to sulfur dioxide atmosphere. The aim was to create a new approach to forecasting corrosion behavior without the necessity of an extensive use of laboratory testing using chemical-physical technologies, while taking into account the main chemical reactions. Although the model was kept simple, just describing the main reaction and transport processes involved, the mathematical simulations and the related model calibration are in agreement with the laboratory experiments. References [1] D.A. Scott, Copper and Bronze in Art: Corrosion, Colorants, Conservation, The Getty Conservation Institute, Los Angeles, 2002, ISBN 0-89236-638-9. [2] S.K. Chawla, J.H. Payer, The early stage of atmospheric corrosion of copper by sulfur dioxide, J. Electrochem. Soc. 137 (1) (1990). [3] K. Nassau, A.E. Miller, T.E. Graedel, The reaction of simulated rain with copper, copper patina, and some copper compound, Corros. Sci. 27 (1987) 703– 719. [4] K.P. Fitzgerald, J. Nairn, G. Skennerton, A. Atrens, Atmospheric corrosion of copper and the colour, structure and composition of natural patina on copper, Corros. Sci. 48 (2006) 2480–2509. [5] T.E. Graedel, K. Nassau, J.P. Franey, Copper Patinas Formed in the atmosphere-I. Introduction, Corros. Sci. 27 (7) (1987) 639–657. [6] K.P. Fitzgerald, J. Nairn, A. Atrens, The chemistry of copper patination, Corros. Sci. 40 (12) (1998) 2029–2050. [7] F. Clarelli, A. Fasano, R. Natalini, Mathematics and monument conservation: free boundary models of marble sulfation, SIAM J. Appl. Math. 69 (1) (2008) 149–168.
4816
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
[8] I. Chadam, A. Peirce, P. Ortoleva, Stability of reactive flows in porous media: coupled porosity and viscosity changes, SIAM J. Appl. Math. 51 (3) (1991) 684–692. [9] J.D. Evans, J.R. King, The Stefan problem with nonlinear kinetic undercooling, Q. J. Mech. Appl. Math. 56 (1) (2003) 139–161. [10] T. Fatima, N. Arab, E.P. Zemskov, A. Muntean, Homogenization of a reaction-diffusion system modeling sulfate corrosion of concrete in locally periodic perforated domains, J. Eng. Math. 69 (2–3) (2011) 261–276. [11] A. Muntean, M. Boehm, J. Kropp, Moving carbonation fronts in concrete: a moving-sharp-interface approach, Chemical Eng. Sci. 66 (2011) 538–547. [12] B. De Filippo, L. Campanella, A. Brotzu, S. Natali, D. Ferro, Characterization of bronze corrosion products on exposition to sulphur dioxide, Adv. Mater. Res. 138 (2010) 21–28. [13] J.H. Payer, Corrosion processes in the development of thin tarnish films, in: Proceedings of the Thirty Sixth IEEE Holm Conference on Electrical Contacts Meeting Jointly with the Fifteenth International Conference on Electrical Contacts, IEEE, Piscataway, NJ, 1990, pp. 203–211. [14] M. Watanabe, Y. Higashi, T. Ichino, Surface observation and depth profiling analysis studies of corrosion products on copper exposed outdoors, J. Electrochem. Soc. 150 (2003) B37–B44. [15] I. Sandberg, I. Odnevall Wallinder, C. Leygraf, N. Le Bozec, Corrosion-induced copper runoff from naturally and pre-patinated copper in a marine environment, Corros. Sci. 48 (2006) 4316–4338. [16] C. Chiavari, E. Bernardi, C. Martini, F. Passarini, F. Ospitali, L. Robbiola, The atmospheric corrosion of quaternary bronzes: the action of stagnant rain water, Corros. Sci. 52 (2010) 3002–3010. [17] I.D. Mac Load, Bronze disease: an electrochemical explanation, ICCM Bull. 7 (1) (1981) 16. [18] T. Aastrup, M. Wadsak, C. Leygraf, M. Schreiner, In situ studies of the initial atmospheric corrosion of copper influence of humidity, sulfur dioxide, ozone, and nitrogen dioxide, J. Electrochem. Soc. 147 (7) (2000) 2543–2551. [19] A. Kratschmer, I. Odnevall Wallander, C. Leygraf, The evolution of outdoor copper patina, Corros. Sci. 44 (2002) 425–450. [20] I. Odnevall, C. Leygraf, Atmospheric corrosion of copper in a rural atmosphere, J. Electrochem. Soc. 142 (1995) 3682–3689. [21] M. Watanabe, M. Tomita, T. Ichino, Characterization of corrosion products formed on copper in urban, rural/coastal, and hot spring areas, J. Electrochem. Soc. 148 (12) (2001) B522–B528. [22] T.E. Graedel, Gildes model studies of aqueous chemistry I. Formulation and potential applications of the multi-regime model, Corros. Sci. 38 (12) (1996) 2153. [23] J. Tidblad, T.E. Graedel, Gildes model studies of aqueous chemistry. III. Initial So2-induced atmospheric corrosion of copper, Corros. Sci. 38 (12) (1996) 2201. [24] J.H. Payer, G. Ball, B.I. Rickett, H.S. Kim, Role of transport properties in corrosion product growth, Mater. Sci. Eng. A198 (1995) 91. [25] R.S. Larson, A physical and mathematical model for the atmospheric sulfidation of copper by hydrogen sulfide, J. Electrochem. Soc. 149 (2) (2002) B40– B46. [26] C. Chiavari, A. Colledan, A. Frignani, G. Brunoro, Corrosion evaluation of traditional and new bronzed for artistic casting, Mater. Chem. Phys. 95 (2006) 252. [27] M. Briani, R. Natalini, G. Russo, Implicit-explicit numerical schemes for jump-diffusion processes, Calcolo 44 (2007) 33–57. [28] S.J. Ruuth, Implicit explicit methods for reaction-diffusion problems in pattern formation, J. Math. Bio. 34 (1995) 148–176. [29] B. De Filippo, A novel approach for the study of the characterization of the corrosion patina on bronze samples exposed to sulphur dioxide corrosion atmosphere, Ph.D Thesis, Università of Rome La Sapienza, 2011.
Contents lists available at ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
Mathematical model of copper corrosion F. Clarelli a,⇑, B. De Filippo b, R. Natalini b a Istituto per le Applicazioni del Calcolo ‘‘M. Picone’’, Consiglio Nazionale delle Ricerche, Polo Scientifico – Edificio F CNR, Via Madonna del Piano 10, I-50019 Sesto Fiorentino (FI), Italy b Istituto per le Applicazioni del Calcolo ‘‘M. Picone’’, Consiglio Nazionale delle Ricerche, via dei Taurini 19, I-00185 Roma, Italy
a r t i c l e
i n f o
Article history: Received 29 November 2012 Received in revised form 19 December 2013 Accepted 25 March 2014 Available online 2 April 2014 Keywords: Free boundary model Parabolic problems Finite difference methods Corrosion Copper Brochantite
a b s t r a c t A new partial differential model for monitoring and detecting copper corrosion products (mainly brochantite and cuprite) is proposed to provide predictive tools suitable for describing the evolution of damage induced on bronze specimens by sulfur dioxide (SO2) pollution. This model is characterized by the movement of a double free boundary. Numerical simulations show a nice agreement with experimental result. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction Deterioration of copper and bronze artifacts is one of the main concerns for people working in cultural heritage [1]. More specifically, bronze, a copper-tin alloy, has been widely employed for daily-use and artistic purposes from the Bronze Age up to present. Conservation studies are based on a knowledge of the environmental conditions to which copper and copper alloys may be exposed and include all the information on the material technologies and the nature of the corrosion films or patina, which cover their surfaces. In particular a significant effort has been devoted to study the corrosion due to environmental conditions, such as temperature, moisture and the concentration of pollutants [2–4]. Although in recent years air pollution in European urban areas has decreased considerably, there still remain concentrations of pollutants such as sulfur dioxide (SO2) from combustion of fossil fuels, being one of the most important factors in the deterioration of bronze. Indeed SO2, mixed with water vapor, reacts to produce sulfuric acid (H2SO4), which causes corrosion phenomena on copper surfaces and produces several corrosion products as basic copper sulphates, such as antlerite, posnjakite, brochantite [5]. The latter is the final product of several reaction steps, which can be approximated by two main chemical reactions: cuprite formation, which occurs after a few weeks of exposure to atmospheric conditions, and brochantite formation, which is the final reaction step (for more details on chemical background see Section 2.1) [6]. The complexity of corrosion processes creates the necessity for a quantitative model approach to develop predictive tools, which simultaneously provide both quantitative information as well as simulations of the various processes involved. These methods, similar to those introduced in [7], are useful for the monitoring and detection of surface alterations even before they are visible, making it possible to determine optimal intervention strategies. In this paper we introduce a new partial ⇑ Corresponding author. Tel.: +39 055 522 5813. E-mail addresses: [email protected] (F. Clarelli), b.defi[email protected] (B. De Filippo), [email protected] (R. Natalini). http://dx.doi.org/10.1016/j.apm.2014.03.040 0307-904X/Ó 2014 Elsevier Inc. All rights reserved.
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
4805
differential model, which is used to describe the evolution of damage induced on a bronze specimen by atmospheric pollution. It is based on fluid dynamical and chemical relations and it is characterized by a double free boundary: one between copper and cuprite, the other between cuprite and brochantite. A similar approach was first proposed in the seminal paper [8], and then explored in many other works, see for instance [9–11]. In our paper we apply this approach to metal corrosion, and we validate it by calibrating the model according to the experimental results in [12]. The paper is organized as follows: in Section 2.1 we analyze the main chemical corrosion phenomena and in Section 2.2 we review the main mathematical models already proposed in literature. The Section 3 is entirely focused on the description of the model’s equations and the numerical schemes used, meanwhile in Section 4 we describe the experimental setting and the calibration of the model. Finally, in Section 5, we present the main results related to the simulations produced by our model and in Section 6 the related conclusions. 2. Modeling backgrounds 2.1. Chemical backgrounds Most of the oxidation processes occurring on bronze metal artifacts, when exposed to environmental conditions, are electrochemical and involve interactions between the metal surface, the adsorbed moisture and various atmospheric gases (SO2, CO2, NOx , hydrocarbons) [13]. Electrochemical corrosion processes in electrolyte and condensed moisture layers have been the subject of extensive studies, based on numerous and different approaches [14–16]. When exposed to the atmosphere, copper and its alloys form a thin layer of corrosion, from brownish-green up to greenish-blue colors, which is designated as patina. In the case of copper in a low pollutant levels atmosphere, a native cuprite (copper (I) oxide or Cu2O) film of approximately a few nanometers thick, protects the metal surface from further oxidation. The general reaction is well described in literature, where, in aerated solution, copper can dissolve electrochemically forming copper (I) oxide formation, due to the reaction of copper with oxygen. It is represented by the following schematic reaction [17]:
1 2Cu þ O2 ! Cu2 O: 2
ð1Þ
However, in an aggressive environment, like urban atmosphere, the protective nature of this oxide layer is altered and there is the formation of a non-protective, multi-component, tarnish layer. When the copper is exposed to humidity and sulfur dioxide, three kinds of basic copper sulfate hydroxide are mainly produced: Brochantite Cu4SO4(OH)6, which is a wellknown patina constituent, or other similar products like Antlerite or Posnjakite. In the following we focus our attention to the formation of Brochantite, which is the main observed product. Gaseous sulfur dioxide and sulfate particles are deposited on the electrolyte on the cuprite. Their deposition reduces the pH of the adsorbed water, and this promotes the dissolution of cuprous ion (Cu+) and its oxidation, thus forming cupric ions (Cu2+). In detail, copper (I) ions in solution disproportionate to give copper (II) ions and a precipitate of copper
2CuþðaqÞ ! Cu2þ ðaqÞ þ CuðsÞ : When the cupric and sulfate ion concentrations in the electrolyte are high enough to form brochantite, this phase starts to precipitate on the cuprite [18,19]. In [20] it is indicated that, in the initial oxidation process, cuprite formation is followed by posnjakite, as a precursor phase to brochantite, see also [21], but we are going to neglect this intermediate transformations, due to the elevated speed of the reaction with respect to the time scale of the mathematical model which we are going to present. 2.2. Existing mathematical models Graedel and his collaborators [22,23] have studied atmospheric copper sulfidation at AT&T Bell Laboratories in both experimental investigation and physical or mechanistic model development; furthermore a systematic investigation of copper sulfidation kinetics has been performed. In the paper [23] Tidblad and Graedel developed a model able to describe the SO2 copper corrosion. Their work is based on aqueous chemistry and without considering spatial dimensions. In the paper [24], they introduce a model which describes the corrosion of copper exposed to moist air with a low SO2 concentration. Here, they consider four stages in the development of a corrosion patina: the metal, a non-protective oxide film which has a high ionic transport property, an outer layer of corrosion products which can permit the penetration of water and gases and an external absorbed water layer. More recently, Larson proposed in [25] a different model that describes the atmospheric sulfidation of copper by proposing a physical copper-sulfidation model that includes four distinct phases: the substrate metal, a cohesive cuprous sulfide (Cu2S) product layer, a thin aqueous film adsorbed on the sulphide and the ambient gas. Larson postulated that transport through the sulfide layer occurs via diffusion and electromigration of copper vacancies and electron holes. Later, in the 1990s, copper sulfation by SO2 was investigated by Payer et al., who focused their attention on the early stage of corrosion in moist air (75% relative humidity (RH) at 25o C) with a sulfur dioxide concentration of 0:5%. The techniques employed (SEM, AES and TEM) allowed for an analysis and characterization of the oxide film (composed of
4806
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
cuprous oxide, copper sulfate and sulfide) on copper surfaces and for the mechanism of evolution of the corrosion chemistry to be described [2]. In the last decade the effect of sulfur dioxide in acid rain on copper and bronze has been investigated. Robbiola and his co-workers studied the effect of its cyclic action on bronze alloys, underlining the different types of patina formed in ’’sheltered’’ and ’’unsheltered’’ areas of bronze monuments [26,16]. Subsequently, Larson refined the copper-sulfidation model by focusing on the transport of charged lattice defects in a growing Cu2 product layer between the ambient gas and the substrate metal. As previously, this transport is postulated to occur via both diffusion and electromigration. 3. The model In this section, we aim to introduce a mathematical model able to describe the corrosion effects on a copper layer, which is subject to deposition of SO2. The present model is based on the mathematical approach used in [7]. We assume to have a copper sample on which is formed a non protective oxide layer (Cu2O), and, over this layer, a corrosion product (brochantite) grows. Over the brochantite layer is assumed to be the atmospheric air with SO2. An example of these three layers can be seen in Fig. 1. The reaction producing cuprite Cu2O, can be approximated by Eq. (1). Namely, two moles of copper combined with one-half mole of oxygen produce cuprite. Laboratory tests, with high concentration of SO2 and high relative humidity, near to 100%, show that brochantite is the primary product of the reaction, thus it can be considered as the final state reached by the whole reaction. The overall simplified reaction of brochantite formation can be approximate by the following reaction (2)
3 2Cu2 O þ SO2 þ 3H2 O þ O2 ! Cu4 ðOHÞ6 SO4 ; 2
ð2Þ
where two moles of cuprite combined with three moles of water and three-halves of oxygen produce one mole of brochantite. In the following, we assume that these two reactions are instantaneous; thus, we obtain a sharp free boundary between cuprite and the unreacted copper, due to the reaction (1), and a second free boundary between brochantite and cuprite, due to the reaction (2). Assuming these two reactions as instantaneous, the effective time of reaction is implicitly included in the diffusivity coefficients. 3.1. Swelling We define the boundary between copper and cuprite by aðtÞ, the boundary between cuprite and brochantite by bðtÞ and the boundary between brochantite and external air by cðtÞ. Summarizing, the geometry of our problem is one-dimensional, and we have 4 regions (e.g. see Fig. 2): 1. 2. 3. 4.
Copper (inner region). Cuprite Cu2O, between aðtÞ and bðtÞ. Brochantite Cu4 (OH)6 SO4, between bðtÞ and cðtÞ. Laboratory controlled atmosphere on the external side of cðtÞ.
Fig. 1. Example of cuprite and brochantite deposition on a copper sample.
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
4807
Fig. 2. Example of cuprite and brochantite theoretical growth (in time) on a copper sample.
When we speak about three dimensional quantities (e.g. volumes), we consider just one dimension of them, assuming a standard cross sectional area. It is assumed that, in the climatic chamber at a 100% of RH, a thin film of water (with SO2 dissolved) is formed over the external boundary, and it plays an important role in both reactions. The reaction (1) implies a production of Cu2O on the boundary between copper and the oxide layer. The transformation of copper into Cu2O proceeds if sufficient oxygen is available, and it is accompanied by a volume change (swelling rate). The swelling rate can be calculated easily, because the molar ratio in reaction between Cu and Cu2O is 2 : 1. Thus two moles of copper change into one mole of Cu2O, and a different volume of the new matter formed is obtained. The derivative with respect to the time of the swelling s1 of reaction (1) is
s_ 1 ¼ xp a_ ;
ð3Þ
where
xp ¼
lc 2lp
1;
ð4Þ
represents the variation of volume ratio, lc and lp being the molar density (moles/cm3) of Cu and Cu2O respectively. In the present case xp > 0, since it represents an expansion of the material. It is assumed that lc and that lp are constant (i.e. they are homogeneous materials). Under these assumptions, if s1 ð0Þ ¼ að0Þ ¼ 0, then it is possible to conclude that
s1 ðtÞ ¼ xp aðtÞ:
ð5Þ
On the boundary bðtÞ, between cuprite and brochantite, SO2 reacts with Cu2O in presence of water and oxygen, following the reaction (2). Here, 2 mol of cuprite change in 1 mol of brochantite, wasting 1 mol of SO2, 3 mol of H2O and 3=2 moles of oxygen. Thus, we have a cuprite production from the reaction (1) and a cuprite consumption from the reaction (2). Thus, the physical boundary between cuprite and brochantite bðtÞ is given by
bðtÞ ¼ bðtÞ þ s1 ðtÞ:
ð6Þ
Here, bðtÞ is the consumption of cuprite (assumed to be positive in our case) due to the reaction (2), and the derivative with respect to the time of the swelling rate of brochantite s2 is given by
_ s_ 2 ðtÞ ¼ xb bðtÞ;
ð7Þ
where
xb ¼
lp 2lb
1;
ð8Þ
and lb is the molar density of brochantite assumed constant. The variation of the overall external boundary c_ is given by
_ c_ ¼ xp a_ ðtÞ xb bðtÞ:
ð9Þ
3.2. Equation of the model It is well known that relative humidity plays a key role in regulating the speed of oxydation and sulfation. It has been observed that when RH exceeds some threshold, then SO2 reacts completely such as the oxidation of copper happens with
4808
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
full speed. We can interpret this phenomenon as follows. According to Eq. (2), when a molecule of SO2 comes in contact with cuprite, it reacts if three molecules of H2O are available at the same point (we suppose that there is always enough O2). To be more precise, this is true only if vapor condenses on the unreacted specimen surface forming a liquid film, and this happens with high RH values. The brochantite formation (Eq. (2)) has been assumed to develop on the cuprite layer, due to the SO2, H2O and O2, which move through the brochantite layer and react with Cu2O. This reaction implies a wasting of Cu2O, with the formation of a new layer of brochantite. Summarizing, the region of brochantite is given by cðtÞ 6 x 6 bðtÞ and the region of Cu2O is bðtÞ 6 x 6 aðtÞ. We define the concentration of SO2 in the pores of brochantite by S, the water concentration by W and the oxygen concentration by O throughout brochantite, and by G throughout cuprite. The flow of SO2 relative to air is governed by Fick’s law. Thus, in the frame of reference where copper is at rest, the SO2 flux has the following expression
@S @S J s ¼ nb Ds Sxp a_ Sxb b_ ¼ nb Ds þ Sc_ ; @x @x
ð10Þ
where the first term refers to SO2 diffusion in the brochantite layer, nb is the brochantite porosity, Ds is the diffusivity. The Sxp a_ and Sxb b_ terms refer to the swelling caused by cuprite and brochantite formation respectively. Hence the mass balance of SO2 in the brochantite layer cðtÞ 6 x 6 bðtÞ is
@S @2S @S Ds 2 þ c_ ¼ 0; @t @x @x
ð11Þ
The value of S at the external boundary cðtÞ is the environment SO2 concentration Sa ðtÞ, which is a known function of time
SðcðtÞÞ ¼ Sa ðtÞ;
ð12Þ
We assumed that SO2 reacts totally with Cu2O at the front bðtÞ, thus we have
SðbðtÞÞ ¼ 0:
ð13Þ
Now, we need a further condition (first free boundary). Since the SO2 flux of moles at the boundary bðtÞ is proportional to Cu2O moles consumption, we have
nb
Ds @S 1 qp _ b; ¼ M s @x 2 M p
ð14Þ
where M s ; M p are the molar weight of SO2 and Cu2O respectively, qp is the mass density of cuprite. Similarly, the water flux is
@W @W J w ¼ nb Dw W xp a_ W xb b_ ¼ nb Dw þ W c_ ; @x @x
ð15Þ
where Dw is the water diffusivity. Thus, the mass balance in the brochantite layer is
@W @2W @W Dw þ c_ ¼ 0: @t @x2 @x
ð16Þ
The value of W at the external boundary cðtÞ is the experimental water concentration W a ðtÞ
WðcðtÞÞ ¼ W a ðtÞ;
ð17Þ
Since some moles of water are wasted by the reaction (2), at the front bðtÞ we have
Jw 3 qp _ W _ b þ nb b: ¼ Mw Mw 2 Mp
ð18Þ
Finally, the oxygen flux is
@O J o ¼ nb Do þ Oc_ ; @x
ð19Þ
where Do is the oxygen diffusivity. Thus, the mass balance in the brochantite layer is
@O @2O @O Do 2 þ c_ ¼ 0: @t @x @x
ð20Þ
The value of O at the external boundary cðtÞ is the experimental oxygen concentration Oa ðtÞ
OðcðtÞÞ ¼ Oa ðtÞ: Since the oxygen is also wasted by the reaction (2), at the front bðtÞ we have
ð21Þ
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
Jo 3 qp _ O _ b þ nb b: ¼ Mo Mo 4 Mp
4809
ð22Þ
3.3. Cuprite layer equations We assumed that the reaction between oxygen and copper occurs on the copper-cuprous oxide boundary aðtÞ. Here, we indicate oxygen by G, just to avoid confusion with the previous region. We make a mass balance of oxygen concentration G, and we assume that all oxygen moles arriving on the inner boundary aðtÞ react. This assumption is not really true, in fact the water plays a key role in the speed of reaction, but in our experiments we have a RH near to 100%, and this fact justify our assumption; also, the diffusivity Dg includes implicitly the finite time of reaction. The flux of oxygen is
@G J g ¼ np Dg Gxp a_ @x
ð23Þ
and the mass balance equation is
@G @2G @G Dg 2 xp a_ ¼ 0: @t @x @x
ð24Þ
The value of G at the boundary bðtÞ is given by the value of oxygen on the boundary b, given by Eq. (22).
GðbðtÞÞ ¼ OðbðtÞÞ;
ð25Þ
Since oxygen reacts totally with copper at the boundary aðtÞ, we have
GðaðtÞÞ ¼ 0;
ð26Þ
Now, we need a further condition to close the system (second free boundary). Since the oxygen react totally at the boundary
aðtÞ, the copper moles wasted are given by np
Dg @G 1 qc ¼ a_ : M g @x 4 M c
ð27Þ
3.4. Rescaling The geometry of the problem is given by two regions, one is given by x 2 ½bðtÞ; aðtÞ which describes the layer of Cu2O, the other is given by x 2 ½cðtÞ; bðtÞ i.e. the brochantite layer. Also we have c_ ¼ xp a_ þ xb b_ . 3.4.1. New variables In the first region we adopt ðx; tÞ ! ðy; sÞ, in the second one ðx; tÞ ! ðz; sÞ. We find out
y¼
z¼
xb
;
@ x ¼ 1=ða bÞ@ y ;
@ xx ¼ 1=ða bÞ2 @ yy ;
xc ; bc
@ x ¼ 1=ðb cÞ@ z ;
@ xx ¼ 1=ðb cÞ2 @ zz ;
ab
@ t ¼ f ðtÞ@ y þ
@ t ¼ qðtÞ@ z þ
1 @s; tr
1 @s; tr
ð28aÞ
ð28bÞ
where
f ðy; sÞ ¼
1 yð@ s b @ s aÞ @ s b ; tr ab
ð29aÞ
qðz; sÞ ¼
1 zð@ s c @ s bÞ @ s c : tr bc
ð29bÞ
Now, we rescale the following parameters and variables in non-dimensional form:
^ ¼ W ; W Wr ^ g ¼ t r Dg ; D k2
^¼ G; G Gr
^S ¼ S ; Sr
^ w ¼ t r Dw ; D k2
a b b c a^ ¼ ; b^ ¼ ; b^ ¼ ; c^ ¼ ; k
^ s ¼ t r Ds ; D k2
k
^ðz; sÞ ¼ t r q; q
k
k
^f ðy; sÞ ¼ t r f
ð30aÞ
ð30bÞ
From now on, we assume that dotted variables are derivatives with respect to the dimensionless time s ðV_ ¼ dV=dsÞ. ^_ the non-dimensional system is: ^ s Mp Sr , Cw ¼ 3 1 Mw qp , Xg ¼ 4np D ^ g Mc Gr . Also, c ^_ ¼ ðxp a ^_ þ xb bÞ, Assuming Xs ¼ 2nb D 2 n Mp W r Ms q Mg q p
b
c
4810
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
3.4.2. Outer region For c 6 x 6 b ! 0 6 z 6 1
0 1 ^s ^_ D @ ^S @ 2 ^S @ c @ ^S ^A ; þq ¼ 2 2 @s @z ^c ^ ^c b ^ @z b
ð31Þ
^Sð0; sÞ ¼ ^Sa :
ð32Þ
^Sð1; sÞ ¼ 0;
ð33Þ
@ ^S ^_ ð1; sÞ ¼ b; ^c ^ @z b
Xs
ð34Þ
0 1 ^ ^w ^ ^ ^_ D @W @2W c @W ^A þq ¼ @ ; 2 2 @s @z @z ^c ^ ^c b ^ b
ð35Þ
^ ^ a ðsÞ: Wð0; sÞ ¼ W
ð36Þ
^ ^ D @W _ ^ 3 qp M w b; w ð1; sÞ ¼ c_ b_ W 2nb W r M p ^c ^ @z b
ð37Þ
0 1 ^ ^ ^ ^o ^_ D @O @2O c @O ^A þq ¼ @ ; 2 2 @s @z ^c ^ ^c b ^ @z b
ð38Þ
^ sÞ ¼ O ^ a ðsÞ: Oð0;
ð39Þ
^ ^ @O D _ ^ 3 qp M o b; o ð1; sÞ ¼ c_ b_ O 4nb Or M p @z ^ ^ bc
ð40Þ
3.4.3. Inner region For b 6 x 6 a ! 0 6 y 6 1,
! ^ ^ ^ ^g @G @2G a^_ @G D ^ ¼ þ xp f ðy; sÞ ; 2 @y2 ^ @ s ða @y ^ ^ b a ^ bÞ
ð41Þ
(we have to highlight that a_ is a derivative with respect to the s)
^ sÞ ¼ 0; Gð1;
ð42Þ
^ sÞ ¼ O ^ ; Gð0; b
^ Xw @ G ^ @y ^b a
^_ : ¼a
ð43Þ
ð44Þ
From now on we use the non-dimensional system, so we indicate all non-dimensional terms without hat. 3.5. Numerical solutions To solve our model, we have to set up an appropriate numerical scheme which is able to describe the process in a short range of time (about 40 h), but also simulations of 1 year taking under control the numerical stability. For these reasons we use finite differences schemes with implicit-explicit terms. It can be found in [27,28].
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
4811
3.5.1. Initial conditions Our numerical procedure requires að0Þ – 0 and bð0Þ – að0Þ to avoid singularities in the internal region confining with copper. By physical assumption, we know that að0Þ > 0 because it represents the first copper consumption. In order to initialize properly the external equations, where SO2 is present, we need to set up bð0Þ > 0 and cð0Þ – bð0Þ. 3.5.2. Numerical scheme Indicating explicit term as HðUÞ and implicit term by GðUÞ. We present our system in the external region (Eqs. 31, 35 and 38), which supplemented by boundary and transmission conditions (32)–(34), (36), (37), (39) and (40).
U t ¼ HðUÞ þ GðUÞ;
ð45Þ
where
0
S
1
B C U ¼ @ W A; O the explicit term HðUÞ is
1 _ ðbc cÞ þ q @S @z B C B C c_ @W C B HðUÞ ¼ B ðbcÞ þ q @z C; @ A _ ðbc cÞ þ q @O @z 0
and the implicit term GðUÞ is
0
Ds @2 S ðbcÞ2 @z2
B B Dw GðUÞ ¼ B ðb 2 @ cÞ
@2 W @z2
Do @2 O ðbcÞ2 @z2
1 C C C; A
since GðUÞ is a stiff term which will be integrated implicitly to avoid excessively small time steps. A more general scheme is IMEX-DIRK Runge–Kutta, which is given by, for t ¼ nDt,
uðiÞ ¼ un þ Dt
i1 i X X ~ik H uðkÞ þ Dt aik G uðkÞ ; a k¼1
unþ1 ¼ un þ Dt
i ¼ 1; . . . m
ð46Þ
k¼1
m X
~ i H uðiÞ þ Dt x
i¼1
m X
xi G uðiÞ :
ð47Þ
i¼1
e ¼ ða ~ik Þ, where a ~ik ¼ 0 for j P i and A ¼ ðaik Þ are m m matrices, x ~ i and xi are m-dimensional vectors, Here, the matrices A such that the resulting scheme is explicit in H and implicit in G. The DIRK formulation requires aik ¼ 0 for j > i [27]. Following the IMEX formalism we shall use the notation Name (s; m; p) to identify a scheme, where s is the number of stages of the implicit scheme, m is the number of explicit stages and p is the combined order of the scheme. In our case ~ 1; x ~ 2 Þ ¼ ð0; 1Þ and ðx1 ; x2 Þ ¼ ð0; 1Þ. we choose to use an Implicit-Explicit Midpoint (1; 2; 2): s ¼ 1; m ¼ 2 and p ¼ 2, where ðx By the same method, and at the same time-step of previous integration, we integrated also Eqs. (41)–(44). (See Table 1). 4. Assessment of the model 4.1. Experiments The calibration of the model has been made with reference to a precise experimental campaign, which is described in [12], where the main characteristics of the experimental setting are reported. Experiments were carried out in a cyclic corrosion cabinet (Erchisen Mod. 519/AUTO) and the Cu–12Sn cast bronze specimens were chosen as a representative alloy used in the past (ancient Greek type). For the corrosion test, all bronze specimens were exposed to an atmosphere containing about 200 ppm of SO2 at 40 C and 100% RH for 8 h (wet cycle), subsequently they were exposed to room conditions for 16 h (dry cycle). Each wet and dry cycle was repeated 20 times. The measures of patina thickness (see Fig. 3), useful for the model calibration, were performed using a SEM-EDS instrument equipped with an Image Analyzer (IA). Each measurement (after 8, 24 and 40 h) has been obtained by an average of 20 measures. The early stage of exposure of bronze samples mainly produced copper hydroxyl-sulfate (brochantite and chalcanthite), cuprous oxide, and tin sulfide (ottemannite). The tarnish product layer also contained a trace of tin oxide. As determined by XRD analysis, the brochantite formation was predominant with the increase of the number of corrosion cycles. Otherwise,
4812
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816 Table 1 Parameters table. Param.
Value
Dimension
Indications
qp
6:00
g/cm3
Cu2O mass density
Mp
143:09 8:94 63:55 3:97 452:3 1:46 64:07
Cu2O molar mass Copper mass density Copper molar mass Brochantite mass density Brochantite molar mass SO2 mass density SO2 molar mass Reference value
1 106
g/mol g/cm3 g/mol g/cm3 g/mol g/cm3 g/mol g/cm3
Gr
5:86 105
g/cm3
Reference value
Or
5:86 105
g/cm3
Reference value
Wr
2:3 105
g/cm3
Reference value
Oa
2:931 105
g/cm3
Reference value
tr
3:1536 107 15 0:3 0:002
sec
Time reference value
cm nondim. nondim.
Length reference value Porosity of brochantite Porosity of cuprite
qc Mc
qb Mb
qs Ms Sr
k nb np
the chalcanthite formation was detected during the early stage of the corrosion processes and the intensity of its characteristic peaks decreases as a function of time. The same consideration can be made for the ottemannite formation, detected with SEM-EDS in localized micro-crack areas, which are rapidly covered by basic copper sulfates [12,29]. 4.2. Calibration The laboratory corrosion tests have been used to calibrate the model. We suppose that the diffusivity coefficients are the parameters to calibrate by experimental tests. To do that, we used the thickness of corrosion products measured at three fixed times. Each thickness value has been obtained by a sample average and standard deviation, computed at each time, based on twenty random measurements, see Table 2. We used the least square method to find the best parameters. We obtained (in cm2/sec): Dg ¼ 9:9 109 , Ds ¼ 3:96 105 , Do ¼ 9:9 106 and Dw ¼ 3:96 105 . The evolution in time of the corrosion products thickness, using the best parameters, is in Fig. 4. We can see the difference between the experimental points and the best simulation in Fig. 5. Remark. here, we test the numerical scheme to reproduce the conservation of mass implicit in the chemical reactions described in Eqs. (1) and (2). Simulations, after 40 h, give us the following values: c ¼ 9:505 104 (cm), a ¼ 3:1693 104 (cm), b ¼ 7:9916 104 (cm). Knowing the molar volume of copper V c ¼ qc =M c , cuprite V p ¼ qp =M p , and brochantite V b ¼ qb =M b , we find that the number of copper moles wasted are a=V c ¼ 4:4588 105 , the number of cuprite moles formed are hp =V p ¼ 2:2294 105 , i.e. two copper moles are wasted to form one cuprite mole, as we expect from Eq. (1). Similarly, we obtain the number of cuprite moles wasted by reaction (2) is b=V p ¼ 2:2173 105 , while the number of brochantite moles formed by the same reaction is about 1:1086 105 , i.e. for each brochantite mole formed two cuprite moles are wasted, as expected from (2).
Fig. 3. Example of patina thickness measurement with the SEM employment. The green lines measure the thickness of corrosion products, mainly brochantite. From [12]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4813
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816 Table 2 Patina thickness measures. Time of measure (h)
Averaged value (cm)
8
1:7331 104
9:2672 104
1:8473 104
5:4418 10
24 40
Standard deviation
4
4
2:4102 104
13:2522 10
Fig. 4. Evolution of cðtÞ (black line), bðtÞ (blue line) and aðtÞ (red line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
−3
1.6
x 10
Experimental points Simulation
1.4
Thickness of corrosion products (cm)
1.2
1
0.8
0.6
0.4
0.2
0 −5
0
5
10
15
20
25
30
35
40
45
Time (hours) Fig. 5. Simulation (continuous line) and experimental points (in red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4814
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816 −11
4
x 10
3.5
3
SO mass density [g/cm ]
3 2.5
2
2 1.5 1 0.5 0 Jan
Feb Mar
Apr
May
Jun
Jul
Aug
Sep
Oct
Nov
Dec
Fig. 6. SO2 concentration detected at Villa Ada (Rome, Italy) during the 2006.
35 30 25
o
Celsius
20 15 10 5 0 −5 −10 Jan
Feb Mar
Apr
May
Jun
Jul Aug Time
Sep
Oct
Nov
Dec
Fig. 7. Temperature detected at Villa Ada (Rome, Italy) during the 2006.
−5
4
x 10
3.5
Saturated vapor density
SVD (g/cm3)
3 2.5 2 1.5 1 0.5 0 Jan
Feb
Mar
Apr
May
Jun
Jul
Aug Time
Sep
Oct
Nov
Dec
Fig. 8. Saturated vapor density obtained by temperature and RH detected at Villa Ada (Rome, Italy) in 2006.
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
4815
Fig. 9. Simulation of corrosion with environmental data detected at Villa Ada (Rome, Italy) in 12 months (2006).
5. Application In this section we use the calibrated model with environmental data, detected in Rome at Villa Ada during the year 2006. This application is just an example of how the model can be used. Since the model has been calibrated for high RH and SO2 concentration, the calibration should be improved by experiments with lower values of RH and SO2 concentration. The SO2 concentration detected is in Fig. 6 and the temperature is in Fig. 7. Assuming good calibration with our environmental data, we need the saturated vapor density (SVD in [g/cm3]) as function of RH and temperature T ½ C detected at Villa Ada. This is used to obtain an exact quantity of water vapor in the air from relative humidity, the density of water in the air is given by RH SVD ¼ ActualVaporDensity. The behavior of water vapor density is a non-linear function, but an approximate calculation of saturated vapor density can be made from an empirical fit of the vapor density curve in Eq. (48) (between 0–45°C is a good approximation).
SVDðTÞ ¼ 5:018 þ 0:32321T þ 8:1847 103 T 2 þ 3:1243 104 T 3 :
ð48Þ
The figure of saturated vapor density obtained is in Fig. 8. Thus, using SO2, temperature and humidity data detected at Villa Ada in Rome during the 2006, we obtained the behavior of corrosion products illustrated in Fig. 9. Results of this simulation represent a qualitative example of corrosion formation on copper under SO2 attack. They have to be improved by other laboratory experiments, because our model is calibrated using a very high SO2 concentration (never present in the atmospheric environment in this concentration), and under high RH (close to 100%). To have a more accurate calibration, it would be necessary set up mainly experiments with different SO2 concentration and different concentration of humidity in the air, so to calibrate the model under these variations. Our conditions produce a thick cuprite layer in 1-year simulations with respect to the laboratory experiments, because of the air SO2 concentration is lower than that present in the experimental room, and slowing the growth of brochantite layer. 6. Conclusions In conclusion a new mathematical model was developed to describe and simulate the evolution of brochantite formation on Cultural Heritage artifacts exposed to sulfur dioxide atmosphere. The aim was to create a new approach to forecasting corrosion behavior without the necessity of an extensive use of laboratory testing using chemical-physical technologies, while taking into account the main chemical reactions. Although the model was kept simple, just describing the main reaction and transport processes involved, the mathematical simulations and the related model calibration are in agreement with the laboratory experiments. References [1] D.A. Scott, Copper and Bronze in Art: Corrosion, Colorants, Conservation, The Getty Conservation Institute, Los Angeles, 2002, ISBN 0-89236-638-9. [2] S.K. Chawla, J.H. Payer, The early stage of atmospheric corrosion of copper by sulfur dioxide, J. Electrochem. Soc. 137 (1) (1990). [3] K. Nassau, A.E. Miller, T.E. Graedel, The reaction of simulated rain with copper, copper patina, and some copper compound, Corros. Sci. 27 (1987) 703– 719. [4] K.P. Fitzgerald, J. Nairn, G. Skennerton, A. Atrens, Atmospheric corrosion of copper and the colour, structure and composition of natural patina on copper, Corros. Sci. 48 (2006) 2480–2509. [5] T.E. Graedel, K. Nassau, J.P. Franey, Copper Patinas Formed in the atmosphere-I. Introduction, Corros. Sci. 27 (7) (1987) 639–657. [6] K.P. Fitzgerald, J. Nairn, A. Atrens, The chemistry of copper patination, Corros. Sci. 40 (12) (1998) 2029–2050. [7] F. Clarelli, A. Fasano, R. Natalini, Mathematics and monument conservation: free boundary models of marble sulfation, SIAM J. Appl. Math. 69 (1) (2008) 149–168.
4816
F. Clarelli et al. / Applied Mathematical Modelling 38 (2014) 4804–4816
[8] I. Chadam, A. Peirce, P. Ortoleva, Stability of reactive flows in porous media: coupled porosity and viscosity changes, SIAM J. Appl. Math. 51 (3) (1991) 684–692. [9] J.D. Evans, J.R. King, The Stefan problem with nonlinear kinetic undercooling, Q. J. Mech. Appl. Math. 56 (1) (2003) 139–161. [10] T. Fatima, N. Arab, E.P. Zemskov, A. Muntean, Homogenization of a reaction-diffusion system modeling sulfate corrosion of concrete in locally periodic perforated domains, J. Eng. Math. 69 (2–3) (2011) 261–276. [11] A. Muntean, M. Boehm, J. Kropp, Moving carbonation fronts in concrete: a moving-sharp-interface approach, Chemical Eng. Sci. 66 (2011) 538–547. [12] B. De Filippo, L. Campanella, A. Brotzu, S. Natali, D. Ferro, Characterization of bronze corrosion products on exposition to sulphur dioxide, Adv. Mater. Res. 138 (2010) 21–28. [13] J.H. Payer, Corrosion processes in the development of thin tarnish films, in: Proceedings of the Thirty Sixth IEEE Holm Conference on Electrical Contacts Meeting Jointly with the Fifteenth International Conference on Electrical Contacts, IEEE, Piscataway, NJ, 1990, pp. 203–211. [14] M. Watanabe, Y. Higashi, T. Ichino, Surface observation and depth profiling analysis studies of corrosion products on copper exposed outdoors, J. Electrochem. Soc. 150 (2003) B37–B44. [15] I. Sandberg, I. Odnevall Wallinder, C. Leygraf, N. Le Bozec, Corrosion-induced copper runoff from naturally and pre-patinated copper in a marine environment, Corros. Sci. 48 (2006) 4316–4338. [16] C. Chiavari, E. Bernardi, C. Martini, F. Passarini, F. Ospitali, L. Robbiola, The atmospheric corrosion of quaternary bronzes: the action of stagnant rain water, Corros. Sci. 52 (2010) 3002–3010. [17] I.D. Mac Load, Bronze disease: an electrochemical explanation, ICCM Bull. 7 (1) (1981) 16. [18] T. Aastrup, M. Wadsak, C. Leygraf, M. Schreiner, In situ studies of the initial atmospheric corrosion of copper influence of humidity, sulfur dioxide, ozone, and nitrogen dioxide, J. Electrochem. Soc. 147 (7) (2000) 2543–2551. [19] A. Kratschmer, I. Odnevall Wallander, C. Leygraf, The evolution of outdoor copper patina, Corros. Sci. 44 (2002) 425–450. [20] I. Odnevall, C. Leygraf, Atmospheric corrosion of copper in a rural atmosphere, J. Electrochem. Soc. 142 (1995) 3682–3689. [21] M. Watanabe, M. Tomita, T. Ichino, Characterization of corrosion products formed on copper in urban, rural/coastal, and hot spring areas, J. Electrochem. Soc. 148 (12) (2001) B522–B528. [22] T.E. Graedel, Gildes model studies of aqueous chemistry I. Formulation and potential applications of the multi-regime model, Corros. Sci. 38 (12) (1996) 2153. [23] J. Tidblad, T.E. Graedel, Gildes model studies of aqueous chemistry. III. Initial So2-induced atmospheric corrosion of copper, Corros. Sci. 38 (12) (1996) 2201. [24] J.H. Payer, G. Ball, B.I. Rickett, H.S. Kim, Role of transport properties in corrosion product growth, Mater. Sci. Eng. A198 (1995) 91. [25] R.S. Larson, A physical and mathematical model for the atmospheric sulfidation of copper by hydrogen sulfide, J. Electrochem. Soc. 149 (2) (2002) B40– B46. [26] C. Chiavari, A. Colledan, A. Frignani, G. Brunoro, Corrosion evaluation of traditional and new bronzed for artistic casting, Mater. Chem. Phys. 95 (2006) 252. [27] M. Briani, R. Natalini, G. Russo, Implicit-explicit numerical schemes for jump-diffusion processes, Calcolo 44 (2007) 33–57. [28] S.J. Ruuth, Implicit explicit methods for reaction-diffusion problems in pattern formation, J. Math. Bio. 34 (1995) 148–176. [29] B. De Filippo, A novel approach for the study of the characterization of the corrosion patina on bronze samples exposed to sulphur dioxide corrosion atmosphere, Ph.D Thesis, Università of Rome La Sapienza, 2011.