# Numerical modelling of mechanical degradation of canvas paintings under desiccation – Heritage Science

#### ByD. S.-H. Lee, N.-S. KIM, M. Scharff, A. V. Nielsen, M. Mecklenburg, L. Fuster-López, L. Bratasz and C. K. Andersen

Aug 12, 2022

The results from the earlier models showed the effective application of FEM in understanding crack formation, propagation, and later development in canvas paintings. Since computational models in general are mere abstractions of the detailed reality, they must be accepted as approximated solutions. Therefore, results from computational analysis should be verified experimentally.

For the extreme model sizes, the required computational facilities may not be readily available. Hence, many investigations often take the route of developing mathematically valid but highly abstract computational models at first, and gradually increase the complexity closer to reality. The current investigation is based on the same basis of understanding, and it starts by adopting the theoretically valid degradation models for oil paintings on canvas, and further attempts to increase the level of precision by introducing further geometrical and material complexities.

It can be observed from several paintings, which have been exposed to the comparable environmental conditions, that the crack patterns and sizes can vary significantly [21]. One may argue that the current conditions of paintings cannot be discussed without the thorough review of both environmental and physical history of specific objects. Such historical record is, however, rarely complete, for example, for a 200-year-old painting; and arguably, it only takes one incident such as a sudden drop in RH or a poor handling of the painting during transportation to cause significant cracks. Conservators may further argue that the different geometrical and material specifications between the paintings can also incur varying magnitudes of degradations. Indeed, the review of historical paintings can reveal that the range of material compositions and geometrical specifications such as the size of stretcher and the thickness of different layers can be widespread, which can be the reason behind the varying crack patterns. For these reasons, it is always difficult to provide a sharp and definite explanation for the local cracks on a specific painting based on the studied mechanisms, which were initially derived from parametrically generalized models. Consequently, this paper presents a more focused investigation with specific sets of materials in their different compositions.

### Oil paintings on canvas with different material compositions

This paper presents an investigation, which focuses on compositions of four specific material layers: Canvas layer, glue/size layer, ground layer and paint layer. They are considered the basic layers in oil paintings on canvas. In some paintings the glue layer can be absent, but the majority in the western painting tradition are composed with these layers. The specific materials occupying the layers can vary. However, earth pigments were typically used extensively in grounds in the 17th to eighteenth century, while lead white mixed with fillers, such as chalk or barium sulphate was widely used as ground material in the nineteenth century [22], as well as in the twentieth century to some degree [23]. Moreover, the common observation is, that paintings containing earth pigments in the ground exhibit more conservation issues and are in worse state of the preservation comparing to other types of grounds. The materials and the types of layer compositions investigated in this study are summarised in Table 1. They represent a simplified version of the tendency described above, where either an earth colour or a lead white based oil ground was used as the structurally significant ground layer applied on the (sized) canvas. Mechanically, these represent a strong and a weak ground. The presence of glue sizing cannot always be assumed and a case without glue is incorporated to investigate the significance of this layer. The glue layer is important since this layer is mechanically dominant due to its stiffness and strength in the lower RH region and can therefore affect crack initiation in a painting structure.

The material definitions used in the analysis were dependent on two environmental parameters: relative humidity and temperature. However, the research interest limits the material definitions to be functions of RH only, and the temperature is assumed to be constant at 25 °C. The material definitions of the applied materials, such as paint, ground and glue were derived based on the full-relaxation stress state which correlates better to the incurred stresses when the paintings are exposed to slow environmental variations, reflecting for example the effect of natural seasonal change in uncontrolled indoor climate.

With respect to RH, the paint samples took about 48 h to reach full equilibrium with the environment, but took 7–10 days to full stress relaxation. The corresponding material tests were carried out at the Smithsonian Institution, in which 20–21-year-old flake white safflower oil specimen was given the sufficient time to fully relax between strain increments. The results were compared with the results from rapid tests with only 30 s strain interval. The comparison showed that, the material behaved significantly different under full-relaxation condition; when the specimen were given the consistent time to relax between the strain intervals, the measured modulus and the ultimate stress were only 14% and 20% of when strained with only 30 s interval, respectively. Glues took much longer to reach fully relaxed state at the minimum of 160–200 days depending on the thickness [24], and the current study considered the glue at 7 to 10 days relaxation state (Fig. 1).

Glue is the most responsive material among the three materials in the desiccation model with the largest moisture expansion (, αG = 2.6E-4) and the elastic modulus. The effective elastic modulus can be calculated from the equations (Eq. 1, 2, 3) allowing for the reconstruction of the restrained tests at Smithsonian Institution.

$${E}_{G}=1378.8+15.1*left(75-{mu *RH}_{i}right)+(15.1*RH)$$

(1)

where RHi is initial relative humidity value, RH is the relative humidity at the point of interest, and µ is the offset coefficient for RHi, which can be calculated from below:

For RHi > 71,

$$mu =(-0.014*{RH}_{i})*{e}^{((3.5e-5*{{RH}_{i}}^{2} – 0.00499*{RH}_{i} + 0.183) * RH)}$$

(2)

For RHi < 71,

$$mu =left(3.06E-4*{{RH}_{i}}^{2}-0.0228*{RH}_{i}+1.0344right)*{e}^{((3.52E-5*{{RH}_{i}}^{2} – 5.0165e-3*{RH}_{i} + 0.1864) * RH)}$$

(3)

Figure 1 shows the comparison of test results and FEM results of desiccation induced stress in glue samples. The computational analysis was run in Abaqus (Dassault Systèmes, 2018), and the comparison shows that the FEM results are very accurate in predicting stresses in restrained materials under desiccation.

The material definitions introduced in this investigation are given in Table 2 which are based on fully relaxed stress conditions. The Poisson’s ratios of the materials were kept constant at 0.3. For canvas, elastic modulus of 35 MPa and moisture expansion coefficient of 1.0E-4 are used [20, 24].

#### Effect of strain rate on the material property

##### Characteristics of wooden stretcher

The structural constraint in a canvas painting is provided by the stretcher. The distributions of displacements in paintings are related to aspect dimension ratios of the stretchers. The typical dimensional ratios of the stretchers commonly used from 18th Century [31, 33] and onwards are used in current study. Pernety’s ratio of 1:1.2 is confirmed by a study of 38 Danish Golden age paintings [32] (Andersen 2013), thus the canvas size in the current study is set to 762 by 635 mm (30 in. × 25 in.).

Pine and Spruce is the commonly known types of stretchers, and the following material properties of Scots Pine were adopted for the current models (Table 3; the subscription number indicate the transverse strain direction and the stressed direction of the material, while 1, 2 and 3 represents directions in x, y and z-axis, respectively).

The moisture induced deformation of different timber species in conservation has been investigated in the past. However, the understanding about the response time varies [26, 27]. The amount of shrinkage in the wooden stretcher over time affects the induced stress in paint layers and stress in paint layer under desiccation is greater when the shrinkage of the stretcher is minimal. Therefore, in the consideration of worst-case scenario, the wood shrinkage is assumed negligible at the time of the other layers reached the equilibrium state.

### Computational models for global and local behaviours of oil canvas painting

#### Geometrical definitions

Models in two different scales were constructed; global scale and local scale. In the global scale models (global model), the full-size painting (762 mm by 635 mm) was digitally reconstructed in three-dimension, and the local scale models (local model) are two-dimensional cross-section models, which were constructed based on an image of an actual painting.

Two types of global model were tested initially. The first type has paint and ground layers with constant thickness (C-Type), while the thicknesses were varied for the second type (V-type) (Fig. 2). The constant thicknesses of C-Type model are the average of empirically measured values from 37 paint cross sections at the Royal Danish Academy, Institute of Conservation (see Valbjørn Nielsen et al., forthcoming), and the measured values of paint ranged between 0.012 mm to 0.5 mm, and 0.0069–0.4 mm for the ground. The computational model for the V-Type was constructed with randomly generated thickness within the given ranges for each layer at 200 grid points, which were spaced at equal distances. As the result of such geometrical refinement, the mesh was also refined with the approximate element size of 3 mm, which is similar to the maximum element size (5 mm) defined from the mesh independent study of constant thickness model.

Figure 3 shows the induced maximum principal stress (MPS) in painting models of the different thickness types. The painting models had material layer composition 1, which were exposed to RH change from 50 to 10%. In the figure, Type-V has up to 38% more stress than Type-C model. It is to note that, despite the thickness of C-Type is the average value of V-Type, the average of the undulating MPS values in V-Type is not the same, but greater than Type-C. Such result is often seen from analysis of thin structures; when geometrical imperfection is introduced to numerically perfect geometry, the stress level can increase by multiple factors. Similarly, the simulation result indicates that the impact of geometrical undulation (can be understood as introducing the realistic geometrical imperfection to the ‘ideal’ constant thickness model) on the stress development in the paint layer is actually greater than using ideally flat model with the corresponding thickness. Furthermore, the thickness in the model was varied only every 3 mm, and in reality, the thickness can vary as small as every 1 mm. Thus, it is expected that the stress variation can be greater in the real painting.

The immediate advantages of modelling in local scale in addition to the global scale is, firstly, that the analysis can be based on the models with the realistic geometrical definitions and secondly, that the interactions between different material layers can be closely observed for the possible crack formations. For the realistic geometrical definitions, the cross-section model was reconstructed based on the cross-section image of an actual painting (Fig. 4).

As discussed in the earlier section (Table 1), the glue layer is omitted in the study of composition type 2; the above figure is only the depiction of the models when all material layers are present. In the model with glue layer omitted, the void was replaced with ground layer.

#### Computation scheme for desiccation strain in the finite element analysis

The RH induced strains are calculated based on the assigned moisture expansion coefficients of materials (Table 2), and calculation used the following equation:

$${varepsilon }^{f}= {alpha }_{f}left({f}_{beta }right)left({f}_{n}-{f}_{n}^{0}right)-{alpha }_{f}left({f}_{beta }^{I}right)left({f}_{n}^{I}-{f}_{n}^{0}right)$$

(4)

where ({alpha }_{f}) is the material’s moisture expansion coefficient, which is a function of the field variable, ({f}_{beta } or {f}_{beta }^{I}), which in this case is the prescribed relative humidity at different steps. (left({f}_{n}-{f}_{n}^{0}right)) refers to the change in RH at the current RH step, ({f}_{n}) with respect to the reference RH value, at which the strain is assumed to be zero.

Therefore, the induced stress is eventually calculated from the equation below:

$$left[sigma right]=left[Dright](left[varepsilon right]-left[{varepsilon }^{f}right])$$

(5)

where [ɛ] is the total strain and [D] is the material stiffness matrix.

#### Boundary condition

The global models have been restrained in Z-axis translation at the four corners of the stretcher, while all the models have one more translation constraint in X and Y-axes at the centre point (Fig. 5). These nodal constraints were only applied to prevent the unwanted out-of-frame distortions of the entire model due to unbalanced forces, and do not have impact on the shrinkage behaviours of the elements.

The displacement boundary conditions are applied to the outer edges of the section model as shown in Fig. 6, where the horizontal displacements are applied. The different magnitudes of the displacements were obtained from simulations of the global models.

#### Numerical analysis of local crack formations

In the analysis of full oil canvas painting under RH changes, the specific moments and locations of crack formations have not been discussed in depth. For FEM, the simulations of crack propagation are time-consuming process since the crack paths must be followed and the mesh refinement around the crack-tips must be updated every step increment. Therefore, the current investigation adopts the framework of extended finite element method (XFEM) to approximate and understand the underlying mechanisms of crack initiation and propagation in the local scale. As the name implies, XFEM can be understood as extending the existing FEM through additions of the enrichment functions, which typically consist of discontinuous basis functions to standard polynomial functions, and also the near-tip asymptotic functions for the singularity around the crack tip.

However, the current crack propagation simulations were executed, in which traction-separation laws based cohesive segments damage modelling, which can follow crack developments in the element based on the calculated criteria, thus the near-tip asymptotic singularity functions (the right most term in Eq. 6) can be omitted [28]. The main reason behind the choice of cohesive segments method was, that the method better suits the intention of the study to examine crack initiation at any possible location, and the propagation of the initiated cracks are not dependent on predefined paths nor element boundaries [28]

$$uleft(xright)= sum_{i=1}^{n}{N}_{i}left(xright) {[u}_{i}+H(x){a}_{i}+sum_{alpha =1}^{4}{F}_{alpha }(x){b}_{i}^{alpha }]$$

(6)

where ({N}_{i}) is the shape function associated to node i, ({u}_{i}) is nodal displacement vector at node i, ({a}_{i}) is enriched degrees of freedom, and (H(x)) is the Heaviside function as given below:

$$Hleft(xright)=left{begin{array}{c}+1, when (x-{x}^{*})cdot nge 0\ -1, elseend{array}right.$$

(7)

where (x) is an integration point in the element, ({x}^{*}) is the point on the crack closest to (x), and (n) is the unit outward normal to the crack at ({x}^{*}). The initiation of crack occurs when the provided maximum principal stress is reached. For the propagation, the damage evolution model based on energy is used, where the energy value is calculated from the empirical data.

It is often promoted that XFEM does not require extensive mesh refinements as seen in FEM for fracture analysis. However, preliminary studies showed the effect of mesh size on the results from XFEM, especially for the location of first crack formation and the length of crack propagation. The mesh independency study confirmed, that the time of first crack formation could be comparable, while the exact locations and the size of the cracks can differ with the element sizes; mainly due to the fact that the traction–separation cohesive behaviour makes the cracks run across the full size of the elements when the failure criteria is met. The mesh independence study showed that the models with the maximum element size below 0.005 mm didn’t show further stress discrepancies at the measured locations.