Creep in oak material from the Vasa ship: verification of linear viscoelasticity and identification of stress thresholds

Creep deformation is a general problem for large wooden structures, and in particular for shipwrecks in museums. In this study, experimental creep data on the wooden cubic samples from the Vasa ship have been analysed to confirm the linearity of the viscoelastic response in the directions where creep was detectable (T and R directions). Isochronous stress–strain curves were derived for relevant uniaxial compressive stresses within reasonable time spans. These curves and the associated creep compliance values justify that it is reasonable to assume a linear viscoelastic behaviour within the tested ranges, given the high degree of general variability. Furthermore, the creep curves were fitted with a one-dimensional standard linear solid model, and although the rheological parameters show a fair amount of scatter, they are candidates as input parameters in a numerical model to predict creep deformations. The isochronous stress–strain relationships were used to define a creep threshold stress below which only negligible creep is expected. These thresholds ranges were 0.3–0.5 MPa in the R direction and 0.05–0.2 MPa in the T direction.


Introduction
The creep properties of conserved archaeological wood material in large artefacts have not been much investigated despite the importance for the dimensional stability and long-term preservation. For example, the Vasa ship in Stockholm, the Bremer Cog in Bremerhaven, Mary Rose and HMS victory in Portsmouth and the Oseberg viking ship in Oslo, are all carrying substantial loads with different types of support structures which are aimed to minimize stresses in the hull (Jones et al. 1986;Cederlund and Hocker 2007;Hoffmann 2009). However, long-term effects from creep strains have not been explicitly taken into consideration in these support constructions. In the case of the Vasa ship, increasing deformation of joints and load carrying members have been observed. This is to a considerable extent due to softening effect from polyethylene glycol (PEG) used as an impregnation agent to mitigate shrinkage on drying the waterlogged wood (Hocker et al. 2012;van Dijk et al. 2016). In addition, partial degradation of wood polysaccharides in the PEG-free interior of the beams is known to affect the mechanical properties (Bjurhager et al. 2012). Therefore, and as a result of planning a new support system for the Vasa, the need for more detailed information on the mechanical properties and in particular on the creep behaviour has been identified. Previous studies have shown that the main effects of PEG on the Vasa oak compared to recent oak are increased mass density, increased hygroscopicity of the dried material (for low molecular weight PEG) and material softening. This is because PEG acts as a plasticiser (Hoffmann 2013) and reservoir of water in the PEG-rich outer layers of the wood. Therefore, variable ambient conditions may have a high impact on moisture content and variation in mass (Vorobyev et al. 2019), with effects on creep behaviour. Nonetheless, museums tend to showcase their ships in controlled climate (e.g., Hocker 2010), or at least with smaller moisture and temperature variations than in regular buildings. As first approximation, the creep behaviour can be assumed to be an effect of load and time only, neglecting mechanosorptive effects due to a varying climate.
In addressing creep of PEG-treated waterlogged wood in large wooden structures, like shipwrecks, some key questions present themselves. Do these materials show linear creep within the stress levels anticipated in the structure? In a finite element model of the ship, which can potentially be used in the design of improved support structures, one would like to implement a relatively simple material model accounting for creep. This would imply a linear viscoelastic model accounting for creep deformation, where the creep strains depend linearly on the applied stresses within a given range. A second question to address with regard to a lower limit of this range can be formulated: Is there a stress threshold below which no creep is expected of any practical importance? If such a threshold exists, it could serve as target in the design of a support structure. If the support structure can be designed in such a way that the stresses in main load carrying components are relieved down to a level beyond this threshold, creep deformation would effectively not be an issue.
In this study, experimental creep data from a first set of samples from the Vasa have been analysed. The aim was to investigate the linearity of the isochronous stress-strain curves of the creep data as well as to find the viscoelastic parameters and a stress threshold of the samples under creep. The tests were performed at constant climate conditions under three different stress levels for each of the three orthogonal directions of the wooden samples. The standard linear solid (SLS) model was used to describe the material behaviour in one loading direction at a time. The viscoelastic material parameters are important for characterization of the Vasa oak as well as using them in a full-scale numerical model to simulate the actual stresses on the Vasa ship and understand its creep behavior. The approach is considered applicable also to other large wooden structures, specifically PEG-treated waterlogged wooden shipwrecks.

Experimental procedures
The experimental tests were performed in a room with controlled climate such that the relative humidity (RH) was 54 ± 1% and temperature was 20 ± 0.1 °C. Eight cubic samples of the Vasa oak of various origins in the ship prepared for testing in the three orthogonal directions, i.e., longitudinal (L), radial (R) and tangential (T), are presented in Table 1. The dimension of the cubic samples was 11 × 11 × 11 mm 3 . Since mechanical testing is destructive and the material is very precious, only small samples could be spared for mechanical characterization. These were taken primarily from drill cores left when ventilation holes were installed. Although the small specimens mean a larger material variability, one advantage was that the curvature of the annual rings was large compared to the specimen dimensions. The samples could therefore be regarded as orthotropic in a Cartesian sense rather than polar orthotropic, which facilitates a direct identification of the mechanical properties from uniaxial testing along the material axes (Shipsha and Berglund 2007). Two samples were loaded in R direction, four in T direction, and two in L direction. All tests were performed in compression, as this is more common than tension in large structures, where the loads come from the self-weight. The creep deformation was monitored at three stress levels in different orthogonal directions for different samples. According to earlier studies on Vasa wood the longitudinal strength is significantly decreased due to cellulose depolymerisation (Bjurhager et al. 2012). Therefore, more conservative stress levels were chosen in L direction compared to R and T direction. For each sample a camera recorded the deformation with an interval of every 2 s until 100 s, and then every thirty minutes until the end of each period at each stress level. The stress levels were maintained for about fifty days resulting in a total experimental time of 4 months, as shown in Fig. 1. For the sake of traceability, the sample IDs are the same as the unambiguous identification numbers used by the museum.
After recording a series of photographs from one side of each sample, image analysis from software GOM Aramis  1 Experimental set-up for the creep testing: L 1 is 70 cm and L 2 is 5 cm, where the beam L 1 + L 2 is one component stereo system 5 M (Braunschweig, Germany) was used to obtain the average strain of the samples at every time interval. The measurements of the side of the sample were expressed in two perpendicular strains ( y and x ). Since the purpose of the present study is only to investigate creep linearity and the presence of the creep threshold, we confine ourselves to a one-dimensional case, where only y was used, as it represents the strain in the loading direction. Later, the strain curves perpendicular to loading direction ( x ) and corresponding Poisson ratio (υ xy ) calculations are used to identify an orthogonal 3D creep model, which can be implemented for solid-shell elements in a finite element model.

Material model
The following procedure describes the method to obtain the viscoelastic material parameters (hereinafter referred to as material parameters, unless otherwise stated) in loading direction only.
In this study, the material model is chosen as a linear viscoelastic model called standard linear solid (SLS) (Ward and Sweeney 2004). This is one of the simplest creep models, which qualitatively shows the same behaviour as that observed in experiments. The trend of this model consists of an initial elastic deformation followed by decelerating creep deformation slowly approaching an asymptotic strain level for extended periods of testing time. The SLS model has two configurations, which depend on the connection of spring and damper. If spring and damper are in series, it yields a Maxwell model, while connecting a spring and damper in parallel yields a Kelvin-Voigt model of a material. Here, the Maxwell configuration is used. The Kelvin-Voigt configuration is adopted in the paper by Huč and Svensson (2018).
According to linear viscoelastic theory, the relaxation and creep functions are only a function of time, which means that behaviour of the material is independent of strain or stress. The SLS model is a three-parameter model (E r , E m , and τ) that contains a Maxwell arm in parallel with an elastic arm as shown in Fig. 2.
In Fig. 2, E r and E m are the spring parameters of the elastic arm and Maxwell arm, respectively and η is the viscous parameter of the dashpot in Maxwell arm.
The governing equation, obtained in time domain, where strain is a function of applied stress, Δ , and creep compliance modulus J (t) is given by The creep compliance modulus is simplified as follows where the glassy compliance (C g ) , which is the value of creep compliance when the material exhibits glassy behaviour at the beginning of the creep test, the rubbery compliance (C r ) the value of creep compliance at rubbery stage of the material and the characteristic time c are given in the following expressions: The creep model in Eq. (1) in combination with the curve fitting tool in Matlab was used to identify the two viscoelastic parameters of SLS model from the experimental creep data, namely as relative modulus α and relaxation time τ. The former, α, is the relative stiffness loss at infinite time, and the latter, τ, is the time related to the decay rate and defined as the ratio of viscosity (η) to the spring constant of Maxwell arm ( E m ). To calculate α and τ, first C g , C r and τ c were obtained through fitting Eq. (2) to experimental data. Note that Eq. (2) is a simplified version of Eq. (1). After that, E r and E m were obtained using Eqs. (3) and (4). Finally, the parameters α and τ are defined as follows: According to Boltzmann superposition principle (e.g., Ward and Sweeney 2004), as depicted in Fig. 3, at time t 1 , stress Δσ 1 is applied and induced strain is given by

Fig. 2 Standard linear solid model
According to linear viscoelasticity, the compliance J(t) is independent of stress, i.e. it is the same for all the stresses at a specific time. If stress increment Δσ 2 is applied at time t 2 , then the strain increase due to stress increment Δσ 2 can be given by Finally, the strain increase due to Δσ 3 can be written as follows

Results and discussion
In this section, before presenting the viscoelastic parameters for all of the samples, the viscoelastic parameter extraction for a sample curve, verification of linear viscoelasticity according to Boltzmann principle as well as stress threshold identification are demonstrated.

Viscoelastic parameters extraction
In Fig. 4, the strain curve and curve fitting of sample 65901-11 at the first stress level (Δσ = 0.4 MPa) in the T direction are shown. Some transients were observed during load increase followed by a period for the sample to settle in. Therefore, the first three hours of the experimental data were excluded from the curve fitting procedure. The corresponding results, in terms of α and τ for all eight samples are presented in the results section below.

Verification of linear viscoelasticity and identification of a creep threshold
It is well known that the theory of linear viscoelasticity is valid for the low stress levels, which is the condition of this study. As exemplified by sample 65901-11 (T direction), the linearity of the SLS model is demonstrated by the three stress increments Δσ i = 0.4, 0.2 and 0.2 MPa, for i = 1, 2 and 3, respectively. The creep response is shown in Fig. 5, where the effects of step-wise adding of the stress can be clearly observed.
To plot the isochronous stress-strain curves at each point in time, the strain values were taken from the creep curves that were also taken at specific moments in time at different stress level, as shown in Fig. 6. Each stress level can be obtained by adding the previous stress increments. For example, the second stress level for samples in T-direction is σ 2 = Δσ 1 + Δσ 2 = 0.4 + 0.2 MPa = 0.6 MPa. In addition, notice that for a different orthogonal direction, different stress increment was used to cover stress ranges considered relevant for the various directions (cf . Table 1). Usually, the measured strain has large fluctuations. Therefore, to extract the stress and strain at different time points the measured curve needs to be smoothed. Smoothing is a method of reducing the noise within a data set. Different methods such as moving average, Savitzky-Golay filters, and local regression with and without weights and robustness can be used. In this work, the RLOESS (locally weighted non-parametric regression fitting using a 2nd order polynomial) function in Matlab with the span of 0.5 was used to plot the strain as a function of time at different stress levels (Fig. 6). The span is the number of data points for calculating the smoothed value, specified as an integer or as a scalar value in the range (0,1) denoting a fraction of the total number of data points.
By extracting the stress and strain at different time points (t 1 , t 2 and t 3 in Fig. 6) and connecting the data points, it is possible to plot the stress-strain relationship, as shown in Fig. 7. Contrary to conventional stress-strain diagrams, it was chosen to present the data with stress along the abscissa and strain along the ordinate, since the creep compliance is then defined by the slope in the diagram. The data from Fig. 7 and the corresponding compliance moduli for each stress level are also given in Table 2. Due to limitations in time and the scarcity of material, samples sufficient to only 3 stress levels could be spared. Some material was saved by increasing the stress level in a step-wise manner for each sample. It is however then assumed that the previous loading steps do not affect the creep behaviour of the subsequent ones.
According to linear viscoelasticity, the compliance J(t) should be independent of stress, i.e. it is the same for all the stress levels at a particular point in time. Considering the large variation in properties in general for these materials, the relatively small variation in J in Table 2 confirms this hypothesis of linear viscoelasticity for practical purposes within the chosen range of stresses and time.
In addition to independence of creep compliance J with respect to stress level σ, presented in Table 2, the linearity can also be confirmed by the similarity of the isochronous stress-strain curves in Fig. 7.
Normally, a stress threshold in creep is defined by the stress level below which the creep rate, dε/dt, is zero (e.g., Zhang 2010). The limited creep data for the archaeological material does not allow characterization of this type of threshold. Instead, the isochronous stress-strain data were considered, which in the present case take a convex shape including the certain data point in the origin (no stress, no creep strain) in the ε-σ plot. If data points are extrapolated down to the intersection with the stress axis or x-axis, a creep threshold stress can be defined, as schematically shown in Fig. 7. Below this stress value, only negligible creep is expected. If a support structure could be designed to relieve stresses in the wooden artefact down to this level, creep could in principle be ignored.
No creep was detected in the samples loaded in longitudinal direction (Vorobyev et al. 2019). The main reason can be that the samples were loaded under the stress threshold in L direction. In addition, this could be understood from a composite mechanics perspective. In the L direction, the cellulose microfibrils carry most The stiff reinforcing cellulose predominantly oriented in the L direction shows an elastic behaviour, whereas the surrounding lignin-hemicellulose matrix shows a viscoelastic behaviour (Navi and Stanzl-Tschegg 2009). In the following, only the measurable creep in the transverse directions (R and T) is presented. The linearity of the curves for all of the samples in R and T-directions is shown in Fig. 8. Please note that the plots are intentionally plotted as strain (dependent value) as a function of stress (independent value), since the stress is constant and the measured property is the resulting strain, with the slope being the creep compliance. Also note that for some of the samples, the data for time 30 and 40 days were missing. That is why two sets of time points are presented in Fig. 7. However, a sensitivity analysis by choosing different time points showed no effect on the output parameter, i.e. stress threshold evaluation. As a result, a combination of time points was chosen to demonstrate the independence of the stress threshold from the selection of time points, as was expected. It can be seen from Figure 8 that all the samples, with the exception of sample 65900-01, show a linear relationship between stress and strain. By linear fitting of the data, the threshold stress for each orthogonal direction was calculated. The threshold stress in T direction is between 0.05 and 0.2 MPa and in R direction between 0.3 and 0.5 MPa. As expected, R direction has a higher threshold stress due to relatively higher stiffness in comparison to T direction, which can be attributed to the presence of stiffening rays in the R direction (de Borst et al. 2012). Although subject to a relatively large scatter in these ranges, average values of the threshold can be calculated, which were found to be 0.125 MPa in the T direction and 0.4 MPa in the R direction. If creep should be minimised, one can compare stresses from finite element simulations with these values for different relieving support systems.

Viscoelastic parameters
A summary of the viscoelastic parameters, i.e., relative modulus α and relaxation time τ for all of the samples can be found in Table 3. If R and T values are pooled together, the average and standard deviations of these parameters are given in Table 4.
First, some of the values of α and τ in Table 3 are outliers and should be interpreted with care. These are indicated by asterisk symbol. The main reason is that a negligible creep behaviour was observed in those samples (for example, in the case of samples 65901-09, 65904-06 and 65900-01).
Examining the values of α and τ with increasing stress level for each sample, a trend can be perceived. With increasing stress level, the relative modulus α becomes higher and the relaxation time τ becomes shorter. In Table 3, this trend is quite evident for sample 65901-11, 65904-06 (excluding the first load step), for sample 65900-01 (excluding the first load step and only in terms of τ) and for 65903-07 (excluding the last load step).
When the creep results from the various relevant samples are pooled together for each direction and stress level, the same trends can be observed in Table 4, with the exception of the highest stress level in R-direction, where α decreases. However, the scatter is considerable, and the trends cannot be ascertained to any high level of statistical confidence. Overall values for the entire stress range are also presented in Table 4. Although the variation of the average values of the viscoelastic parameters with respect to the stress level may lead to believe that the material cannot be considered to be linearly viscoelastic, the similitude of the isochronous stress-strain relations and the creep compliance (an engineering property) shown in Fig. 8 implies that the assumption of linearity is justified in regard to the generally high variability in material properties. Therefore, the average values of α and τ are calculated and given for two directions of the wood (R and T) in Table 4, which are candidates for material parameters in a finite element model to predict structural creep deformation.
Finally, the results in Table 4 can be compared with some experimental data for other kinds of wood in the literature. Viscoelastic characterization in R and T directions is not widely reported, in particular for archaeological wood. In a study on viscoelastic characterization of wood, Ozyhar et al. (2013) investigated creep on samples made from European  Hu and Guan (2019) for unaged beech wood (Fagus orientalis) in L, R and T directions. From their data, an average relative modulus α of 0.33 and 0.51 for R and T directions, respectively, was obtained. This is only in rough agreement with the present value of α = 0.62 in R-direction but in good agreement with α = 0.58 in T-direction. The relaxation times τ were, however, not comparable for these studies. Polymer materials are known to show multiple relaxation times, which can be attributed to different deformation mechanisms (Ward and Sweeney 2004). The prominence of relaxation time for polymer material depends on the timeframe of the creep test. The chosen test times differed over several orders of magnitude, which means that comparison of relaxation times were not meaningful.

Conclusion
There is a general lack of available creep data for wood materials which could be used for predictions of long-term deformations of structures, although neglecting the timedependent deformations is not justified (Ozyhar et al. 2013;Huč and Svensson 2018). This is particularly the case for archaeological wood materials, where there is little creep data available, despite the creep deformation in large wooden structures in cultural heritage. The results from compressive creep tests of oak samples from the Vasa ship essentially showed linearity in the viscoelastic material model (standard linear solid) according to the Boltzmann superposition principle for small cubic samples loaded in T and R directions. In addition, the viscoelastic properties of the wood as well as a stress threshold, in which negligible creep is expected, were obtained. The creep behaviour of the samples loaded in longitudinal L direction was found to be negligible compared to that in the transverse R and T directions. The obtained rheological parameters may serve as first input parameters in a user-defined material model in a finiteelement analysis for creep deformation predictions of the entire ship. They would however need to be complemented by the coupling between the main material axes, i.e. Poisson numbers, and shear moduli.