A parametric study of crack propagation in paintings caused by temperature and relative humidity cycles based on irreversible cohesive zone model

: The current paper aims to use an irreversible cohesive zone model to investigate the effects of 9 temperature and relative humidity cycles on multilayer thin-film paintings crack pattern. The homogenous 10 one-dimensional paint layers composed of alkyd and acrylic gesso over a canvas foundation (support) with 11 known constant thicknesses are considered as the mechanical model of painting. Experimental data used for 12 mathematical modeling of canvas as a linear elastic material and paint as a viscoelastic material with the Prony 13 series. Fatigue damage parameters such as crack initiation time and maximum loads are calculated by an 14 irreversible cohesive zone model used to control the interface separation. With the increase of the painting 15 thickness and/or the initial crack length, the value of the maximum force increases. Moreover, by increasing 16 the relative humidity (RH) and the temperature difference at loading by one cycle per day, the values of 17 initiation time of delamination decrease. It is shown that the thickness of painting layers is the most important 18 parameter in crack initiation times and crack growth rate in historical paintings in museums and conservation


Introduction
A painting structure consists of support (wood panel and canvas), glue sizing, ground, and paint film (binding media and pigments).Different mechanical properties of the various layers (support, ground, and paint) through aging can lead to craquelure in paintings.Environmental changes such as humidity caused stress in the different paint layers which sometimes leads to strains larger than 1% (restrained layers) and produces a stress-rise in the painting or shrinkage of the glue size during dryness.Accumulation of such failures in cyclic load results in fatigue, plastic deformation (ductile failure), and brittle failure.
The interface between various materials, for example, the interface between a solid gravity dam and the bedrock, is constantly a powerless connection, advancing split initiation and prompting break even under administration loads [1].The irreversible cohesive zone model is an appropriate system to investigate and assess the potential crack at a bi-material interface [2].In light of the irreversible cohesive zone model, some interfacial crack parameters, let's say, crack and break durability were explored through exploratory and numerical studies [3][4].The exploratory examinations showed that the greatness of interfacial unpleasantness would affect the previously mentioned interfacial crack parameters, driving analysts to contemplate its impact by researching examples with smooth interfaces and false notching interfaces [5].The problem of stress analysis of a plate having an elliptical hole (

   
,0 1 2 first time, Griffith (based on thermodynamics, , ,tensile strain bending strain surface load ) proposed the energy-balance concept of fracture ( 0 dU dc  ). Figure 1 presents the schematic of the 2D plane stress problem of a plate having an elliptical hole.Historical paintings in museums are one of the materials that suffer various low-stress cycles which can cause the initiation of a crack or accelerate the crack growth rate [6].For semi-weak materials, the fracture propagation zone lies before the split-tip and pulls in huge concerns when considering the nonlinear reaction of a designing structure built with semi-brittle materials during the crack propagation [7].The impact of the fracture propagation zone on the crack parameters of cement, as a sort of semi-brittle material, has been widely examined over the most recent couple of decades [8][9][10].The size impact of the crack was observed to be connected with the fracture propagation zone properties [11,12], showing that the fracture propagation zone length specifically diminishes quickly when the split causes near the top surface of an example [13].Thusly, the locality crack was observed to be not consistent during the entire crack propagation and rather diminished with the decrease of the fracture propagation zone length [14].Consolidating the hypothetical and exploratory investigations, a bilinear model on neighborhood crack vitality conveyance was proposed to ascertain the genuine explicit break energy [15].Numerous non-linear models have been established to characterize the fatigue parameters such as size, shape, material, and, test method.J-integral is used by many methods.The cohesive crack model used fracture energy, strength in uni-axial, and elasticity modulus [16].As well, the crack band model also uses a width of micro-cracks [17].More fracture-based methods have used the benefit of critical stress intensity factor and critical crack tip opening displacement such as a two-parameter fracture model [18], while the size-effect model for infinitely large test specimens uses critical effective crack length extension (at peak load).Other parameters such as the critical effective extension of crack and critical stress intensity factor (effective crack model), unstable fracture toughness and initial cracking toughness (double-k fracture model ) [19], unstable fracture energy released and initiation fracture energy release (double-g fracture model) [20] are the base of other methods.The cohesive zone method is a classic and simpler method than the above-mentioned methods [16][17][18][19][20].In comparison to linear elastic fracture mechanics or crack tip open displacement methods, it can forecast the un-cracked configurations manners (such as those have blunt notches), larger non-linear zone, and start without initial crack.
The huge impact of a changing fracture propagation zone on solid break attributes and the whole crack propagation has drawn in logical and building networks.The significant examinations have been brought out through test investigations [21][22][23] and numerical simulations [22,23].Furthermore, as one three-dimensional impact on break investigation, a coupled break mode was found to exist in the split thick plate under shear or out-plane stacking, and the power of the coupled model was essentially affected by the thickness of the plate in the three-dimensional limited component analysis [24][25][26].Be that as it may, the examination on the development of the fracture propagation zone during the total crack propagation at a stone solid interface has been minimally detailed.As to shake concrete interfacial crack, it is beneficial to bring up that the determined break energy dependent on crack length without considering the fracture propagation zone is not as much as that dependent on nonlinear crack mechanics [28] by 83 %.Along these lines, it is huge to join the investigation of the fracture propagation zone development at the stone solid interface when investigating the crack system and evaluating the nonlinear reaction of a solid structure developed on bedrock [29].In the interim, the split engendering criteria in numerical strategies have been generally examined, which show the component of break development in semi-brittle materials like cement.Determination of the fracture energy of mortar and concrete using three-point bend tests on notched beams is performed by numerous specialists [30][31].Considering the complex stress distribution at the notch tip under mixed-mode loading, a strain energy density and crack zone model fracture criterion was used to predict the critical load for blunt U-and V-notched brittle specimens [32].
Artists' paintings are composed of polymeric layers.One of the main problems in the polymeric coating materials is to endure mechanical fractures over a continued loading.Painting on canvas made the use of binder for the pigmented paint layers [33].The pigment material provides the color, the binding medium a substance that guarantees that the colored material remains in the applied place.Common paint binders are Steam-pressed linseed oil, Acrylic Resins, and Alkyd resins.The impacts of temperature and changes in the vapor content of air on the mechanical properties of the painting layers in artwork have been examined in many works [34].In painting, which is made of various layers, delamination growth is likely to occur under mixed-mode loading.
Delamination between an alkyd configuration layer and acrylic prepared canvas because of cyclic changes in RH has recently been explored [35].Common environmental control specifications for galleries and museums are relative humidity at 50 or 55 ± 5% relative humidity (RH) and temperature in winter at 19 ± 1 °C and summer at 24 ± 1 °C [36].The work executes the irreversible firm zone model in a limited component examination to display the interface between alkyd paint and prepared the canvas, which results in an alteration to the footing detachment law to represent fatigue failure.Mecklenburg [37] demonstrated that the constituents show diverse dimensional and stress-strain reactions relying upon the ecological conditions.A straightforward order of splits in other works of literature was first efficiently connected by Keck [38].Recently, the type of separation in the interfacial interaction of modern paint layers has been distinguished by Young [39].
Craquelure and interfacial splits have likewise been distinguished, and poor bond characteristics featured when blended-media paints are utilized on canvas for example in a blend of acrylic and alkyd paint layers.Inverse analysis can be employed to optimize the coating design [40][41].
Creep [42] and fracture [43] are classic topics that initially manipulated by linear theories.The nonlinear models of pre-existent crack propagation are developed for the first time by Dugdale [44] and later by Barenblatt [45].
In the present study, a two-layer painting is simulated in 3D stress conditions with the finite element method.This simulation has been carried out in steady-state, isothermal, and single-phase and the effect of temperature variations, layer thickness, and initial crack length on the crack propagation time and the maximum load of the painting are investigated.As well, the distribution of the stress in the painting regions have been studied and evaluated.

Model details and validation
In this study, a single support canvas with a paint membrane is modeled in two dimensions and isothermal.
The schematic of the computational region, including the paint membrane, the canvas layer, and the interface, is illustrated in Fig. 2a.Forces are applied to the top edges at the cracked end and the center of the test specimen.
In this figure, the upper part of the geometry is the painting and the lower part is the canvas.In this figure, the crack of composite structures develops as delamination between plies.As the geometry of the problem is shown in Figure 2a, the layers cracked along a ply interface, and the test specimen is supported at the outermost bottom edges.Because of the symmetry, just half of the test specimen is considered and a symmetry boundary condition is applied.The experimental setup of loading is plotted in Figure 2b.Table 1 shows the geometry parameters of the experimental setup in [27].

 
).The properties of unidirectional laminates composite AS4/PEEK (APC2) which is a carbon fiber reinforced composite used for experimental setup [27].Correspondingly, the geometric, material properties of the laminate composite, and physical properties of the fracture modeling and validation are summarized in Table 1.This specification is derived from the experimental study by Camanho et al. [27].The orthotropic linear elastic properties assume that the longitudinal direction is alongside the global longitudinal direction.The experimental tests are performed by applying different loads in the middle and at the end of the test specimen.
The experimental results relate the load to the displacement of the point of application of the load in the lever (load-point displacement).The lever is not simulated here.In numerical modeling, the cohesive zone elements divided into two groups point and continuous cohesive zone elements.Here to calculate delamination onset and growth surface elements used for cohesive zone elements.The initial crack length is c l .Traction (obey bilinear traction-separation law) linearly increases (with a stiffness Kp ) in anticipation of the ultimate displacement jump (u 0 ) where opening crack reaches a failure initiation.After that the stiffness reductions by an increase of the damage (material softens irreversibly) till the failure at zero stiffness (u f ).In mode I separation displacement is normal to an interface while on mode II and III separation displacement is tangential.
The crack stress modes are:  Mode I (opening mode): the displacements of the crack surfaces are perpendicular to the plane of the crack.
 Mode II (sliding mode): the displacements of the crack surfaces are in the plane of the crack and perpendicular to the leading edge of the crack.
 Mode III (tearing mode): the crack surface displacements are in the plane of the crack and parallel to the leading edge of the crack.
For the current study the mixed-mode mode I, mode II, and mode III separation displacement are included then a combination of separation displacement is used as modeled by Kenane and Benzeggagh [38].Note that for mode I, damage only accumulates positive normal separation, while in shear loading, damage occurs for both shear displacement directions.To test the independence of results from the mesh, nine different meshes of a structured type have been produced and the results are compared with each other.Here meshes with several refinement degrees are used.
Table 2 shows the results of different meshes.From really coarse mesh (only 56 elements) to finer mesh (29928 elements).It is well known that the process zone size is a fraction of this characteristic length, therefore one should not provide calculations for mesh sizes that are larger than 1/10 or 1/20 of this characteristic length.
Although it is nonsense to use such coarse meshes (only 56 elements) regarding the cohesive zone characteristic length, its results could be used as an initial guess for finer cases.Considering the characteristic length of the cohesive zone model, the largest mesh size with 29928 elements is less than 1/200 of characteristic length of the cohesive zone model and all the presented calculations are meaningful accordingly.Percent of relative difference with benchmark solution is the average of percent of the solution with the solution of the 29928 elements case.To check the precision of different meshes, a bi-linear irreversible traction-separation curve is obtained for each mesh; the relative error shown in table 2 presents the average difference of traction-separation achieved from numerical simulation of a sample mesh and traction-separation from an experimental study by Camanho et al. [27].As shown in Table 2, the results are close to each other for the 12364 and 29928 meshes, and finally the mesh number has been selected 29928.Figure 3 reveals multi-grid mesh for two layers from the top view, right view, front view, an isometric view.The finer grid is used in the crack growth region.To increase the accuracy of simulation, given that the damage occurs in the midline of the domain, the mesh of this section is finer, and a multi-grid method is used for meshing (see Fig. 3).Results using irregular meshes are same as regular meshes and irregular meshes are presented in Figure 3 (c).
For validation, a numerical study has been used from the experimental work of Camanho et al, which is shown in Fig. 4. As can be seen, there is an acceptable agreement between the empirical work of [27] and numerical simulation presented here.Although figure 5 does not provide a comparison of data concerning crack propagation, it compares the linear elastic part of the curve for which no crack propagation occurs.As well the start stage of crack propagation is calculated perfectly.As shown material resistance to crack extension firstly is a linear and continuous balance between consumed energy and released energy is maintained during slow stable crack growth.Because of slow stable crack extension, finally, the rising shape is observed (it is flat for truly brittle material).

Governing equations
The stress-strain formula for homogeneous materials is as follows.
(1) 10 1 10 1 while 2nd Piola-Kirchhoff stress tensor is related to the Cauchy stress tensor through the geometric transformation of (4) The dynamics of the displacement of the solid structure is where a Rayleigh damping factor proportional to the stiffness is used for the beam and J is the determinant of F and the deformation gradient tensor is computed from and for the St. Venant-Kirchhoff material the Lagrange strain tensor E is calculated by ( 7) The initial condition of the system is the stationary condition.The presented dynamic formulation in Eq.
(5) would be applied to quasi-static simulations.The quasi-static solution is updated in each iteration of the dynamic solution.In a quasi-static scheme, the transient term is collect the residuals and let the unbalanced of the system to relax in longer times.
Ogden and van der Waals models are usually used in literature for uniaxial tension.Hagan et al. [29] announced that the mechanical reaction of latex paints under uniaxial stacking can be depicted utilizing the hyperelastic, van der Waals model, related to the time needy, viscoelastic Prony arrangement.They used pigments and coloring material of titanium white acrylic gesso and phthalo blue alkyd.Using the viscoelastic model, one can model the creep at constant stress [37], relaxation at a constant displacement, recovery without the stress, constant rate stress, and constant rate strain.
The boundary conditions are the supports fixed in two horizontal directions and the test are monotonic not cyclic.The time-dependent manner of the viscoelastic material is given by the Prony series as (see Table 3 for the constant of Prony Series):

  
where ( ) and the stress as a function of λ ( to uniaxial loading can be derived as follows: (9) The strain vitality capability of the van der Waals model [31] is given by: where α is the chain interaction parameter (0.5), λ m is the locking stretch (for Alkyd is 8 and for Gesso is 10), μ is the initial shear modulus (for Alkydis 75 MPa and for Gesso is 125 MPa) and I is the first stretch invariant which under uniaxial tension is given by: ( I    ); G is the shear modulus.Cyclic loading can be described using the stress amplitude, mean stress, and stress range, respectively.Based on Paris law the fatigue crack growth rate is defined The cohesive zone model assumes a linear relationship between cohesive stress response ( 16) where the energy release rate at maximum loading is larger than the G threshold or one percent of critical energy release rate which is 0.8.As well as the propagation of the delamination is defined by the rate of damage is defined as:

Results and discussion
General numerical solution procedure is to use a Newton-Raphson solution technique to solve the nonlinear system of equations.Other researchers [41][42] presented a Newton-Raphson solution technique to solve the nonlinear system of equations.Though, the Jacobian matrix is unsymmetrical and since is not suitable for FE.The overall procedure for static cohesive crack growth simulation is briefly shortlisted in the next steps: (1) solve the linear system of equations under the external load and calculate the external SIFs, then determine the crack growth; (2) let the damage parameter as a given constant parameter; (3) solve the nonlinear system of equilibrium Equations; (4) calculate the damage parameter under the current external loads for the present crack-tip opening; and (5) go to Step 2 and repeat the calculation for the next step.A quasi-static analysis approach is used for the current viscoelastic modeling.The accuracy of the analysis is controlled by the maximum difference between the creep strain at the beginning and the end of the increment given as Fatigue failures in painting take place in cyclic loading, after a definite time.It also shows evidence of through-thickness and interfacial cracks in paintings under mechanical stresses, as well as the delamination between oil and acrylic paints [32].Here a similar geometry of the experimental setup is used (see Figure 2b).Fig. 5 shows Von mises stress, as combination of other stress (sx, sy, txy ...).In Figure 5, the stress variation in the materials is displayed.Stress contours for various displacements (0.4, 1.8, 4.5, 6 mm) are illustrated in Figure 5.The stress contour with the unit of Pa (N/m 2 ) has the highest values at the crack edge and under the loading place in the middle.The maximum stress across the layer is 250 Pa at the lowest displacement and throughout the highest displacement is more than kPa.Although the stress has increased in both layers, the pressure drop in the painting is larger than the canvas.Since the stress in the painting layer is 2.5 times the stress in the canvas layer.As in the cohesive zone model, the strip plastic zone has cohesive traction equal to yield stress, it is clear to track it while crack propagates through the process in Figure 5.The effect of painting layer thickness and initial crack length is revealed in Figure 7.In Figure 7 a, the painting layer thickness is shown at 20 °C.Regarding the figure, the behavior of force-displacement is similar while the values are different.In contrast to low thickness cases, the increasing behavior of curves is observed after the initial displacement and the highest at the final displacement.In Figure 7 a, the thickness of the painting system is changed from 1 mm to 1 cm.By increase of painting system's thickness, the value of maximum force increases.In Figure 7 b, the painting layer thickness is shown at 20 °C.Regarding the figure, the behavior of force-displacement is similar while the values are different.In all cases, the decreasing behavior of curves is observed after the initial displacement and the highest at the initial of damage displacement.In Figure 7 b, the initial crack length of the painting system is changed from 2 cm to 4 cm.By increase of painting system's initial crack length, the value of maximum force increases.
The curves in Fig. 7a ad 7b are not differ in character.As shown In 7a, the curves tend to increase further after minor drop, but 7b shows no such increase.As seen the delta displacement, at which point is it measured.
The applied force F is similar to the applied force from Fig. 2 and calculated from the location of crack tip.After the initiation, the relationship of the crack extension versus time is found to be approximately linear with a constant extension rate.This is small; highlighting that damage propagation in works of art is likely to be a slow process.Non-destructive inspections could reveal damage and timely corrective action could be taken to allow conservation of fine-art paintings.
Figure 9 shows the effect of crack length on maximum components of stress in numerical modeling of delamination.As shown the maximum stress happens in normal stress and shear stress is not important in this case.The shear stress in this test is lower than normal stress by two orders of magnitude.As the applied load made the system to bend, the produced stress field is anticipated.As shown normal stress (in beam axial direction) is in order of von misses stress.As shown:  By increase of painting system's thickness, the value of maximum force increases.
 By increase of painting system's initial crack length, the value of maximum force increases.
 By increase of loading relative humidity percent at one cycle per day, the values of initiation time of delamination decrease.
 By increase of loading temperature difference at one cycle per day, the values of initiation time of delamination decrease.
The advantages of the present work is to calculate the crack length propagation in painting layers and limitation is constant coefficients.To continue this research the author suggests making accelerated fatigue delamination by humidity and temperature-controlled chamber.The material age, anisotropic behavior of the canvas, are also other parameters that are neglected in the current paper.Finally, the inverse analysis can be employed to optimize the coating design.
Comparison of various components of stress.

Figure 1 .
Figure 1.Schematic of 2D plane stress problem of a plate having the elliptical hole.

Figure 2 .
Figure 2. a) schematic of the computational domain and boundary conditions.b) Schematic of the experimental setup used for validation.

Figure 3 .
Figure 3. Multi-grid mesh for two layers (a) top view, right view, and front view (b) isometric view (c) irregular meshes.

Figure 4 .
Figure 4. Comparison between the numerical simulation results and the experimental study by Camanho et al [27].
is the fracture strain by the opening of the micro cracks, E is Young's elastic modulus,  is its Poisson's ratio, h c is the width of the fracture front, f  is the crack displacement, 0  is the strain at the end of strain-softening at which the micro-cracks coalesce into a continuous crack and z  vanishes.The 2nd Piola-Kirchhoff stress tensor for the material with the Poisson ratio of ν s and Young's modulus of E is defined by(3) δ ⩽ δ initial which the stiffness is changed with cyclic loading to (1-D c ) times.The onset of delamination growth is defined as:(19)

Figure 7 .Figure 8
Figure 7. a) The effect of painting layer thickness, b) The effect of initial crack length.

Fig. 9
Fig.9which shows the force magnitude versus crack extension for various thickness and initial length assesses the geometrical effects.Each data point is obtained by recording the force when the cohesive element at the current crack tip is fully damaged; the crack is then considered to be extended by one element length.

Figure 9 .
Figure 9.Comparison of various components of stress.

Figure 2 a
Figure 2

Figure 4 Comparison
Figure 4

Figure 6 Damage
Figure 6

Figure 7 a 8 a)
Figure 7
X 122.7 GPa Young's modulus, along fibers E Y ,E Z 10.1 GPa Young's modulus, across fibers ν YZ 0.45 -Poisson's ratio, along fibers ν XY =ν XZ 0.25 -Poisson's ratio, across fibers G YZ 3.7 GPa Shear modulus, along fibers G XY ,G XZ 5.5 GPa Shear modulus, across fibersThe main parameters in FEA modeling of the current system are stress (

Table 2 .
The relative difference with benchmark solution for nine different numbers of meshes.

Table 4 .
Material properties of the cohesive zone model interface.

Table 4
shows the material properties of the cohesive zone model interface.Parameters are defined in

Table 4 ,
which are not the same as those presented in Table1as Table1data is for benchmark Table4data is for painting.Both parameters are finally used in the CZM.Based on the analytical method vertical separation of the P denotes the loading force; E is Young's modulus; a is the crack length; B is the width of the beam; h is the thickness of beam; I is the second moment of area of the cantilever beam ( where