Finite element modelling of a historic church structure in the context of a masonry damage analysis

The paper includes a case study of modelling a real historic church using the finite element method (FEM) based on laser scans of its geometry. The main goal of the study was the analysis of the causes of cracking and crushing of masonry walls. An FEM model of the structure has been defined in ABAQUS. A non-linear dynamic explicit analysis with material model including damage plasticity has been performed. A homogenization procedure has been applied to obtain the material parameters used in the modelling of masonry. In the numerical analysis, the interactions between the church structure, the foundations and the ground were taken into account. The obtained results match well with the damaged area of the entire structure from the in situ survey, and it should be highlighted that the proposed FEM model allows for a rather precise identification of the causes and effects of cracking walls in a qualitative sense. Also a brief research summary is presented.


Introduction
As widely commented in the literature, masonry structures behave in a rather complex way [1], exhibiting, e.g., orthotropy along material axes, cracking under tensile stress, and plastic deformation. On the other hand, especially when referring to historic buildings, due to their historical value, restoration and rehabilitation are needed. This usually means that their technical conditions needs to be assessed based on accurate mechanical analyses of the old structures [2]. Numerical modelling by means of FEM is often necessary to perform static and dynamic analyses that would be useful to detect causes of damage [3][4][5][6] and to prepare efficient rehabilitation plans.
The main goal of this study is the numerical modelling and mechanical analysis of a XIVth century church in the delta of the Vistula River in the north of Poland. The church in Gnojewo near Malbork (see Figs. 1 and 2) is one of the oldest churches in Europe with preserved brick nogged timber wall frame construction. It is also one of the best-preserved brick/timber structures in that area. Almost the entire roof, which was built before 1323, has been preserved in excellent condition to the present day. However, the condition of the church, which has been long unused, is worsening, and a number of large cracks have been observed. The technical conditions of the church as indicated by the report drawn up in 2013 [7] are very bad. Part of the components of the structure components are now in an emergency condition. All these components require immediate repair, replacement or reinforcement. Thus, this study is focused on the identification and assessment of the wall conditions and the soil to determine the causes of the cracked walls of the church.
The mechanical analysis of the church model will aid in assessing the main reasons for the damage and will be useful for some restoration and further preservation of the structure. Numerical modelling of masonry structures has been presented and discussed in the literature [8,9,1,10], but the implementation of efficient and ready material models including nonlocal damage of masonry structures [11][12][13] for application in commercial codes remains under development. Among the main problems observed when modelling historic masonry structures is their complex and irregular geometries, which are difficult to define, unknown building materials, mechanical properties and appropriate material models.
The geometry of the Gnojewo church has been defined basing on a scanned points cloud. This was performed using a 3D terrestrial laser scanner (TLS) [3], and the 3D geometry was then defined to be used in FEM-based software for further numerical modelling of its structure. With respect to other geomatic techniques, TLS technology offers a number of advantages, e.g., a capacity to acquire the three-dimensional geometry of an object's surface without direct contact with the structure. This avoids alteration of the material and allows for wide access to structural elements. It also provides high-density data as well as high measurement accuracy. A similar methodology has been already used for historic structures as described in the literature [14][15][16]9].
The properties of the building material were assessed on the basis of laboratory tests. Considering masonry to be a brittle material that has a high strength in compression and residual tensile strength [1], the authors focus their attention on explaining the cracking process based on stress analysis. Assuming that the cracks appear in the area of a tensile stress, a mechanical analysis was performed. The obtained results are discussed, and some possible causes of the damage are indicated by analysing the area of high tensile stress in the computational model of the church.
The main original aspect of the paper is the proposal of numerical detection of damage causes of large historical masonry structures, in nonlinear range, considering foundation (soil) settlement of the masonry structure. This type of approach hardly appears in literature since most of the analyses do not include soil structure interaction in the context of damage detection [17,5,18,11].

Church structure
The analysed church is a brick nogged timber wall frame construction on a stone foundation of boulders. After 1323, the church was rebuilt several times. The main parts of structure and phases of the church reconstruction are shown in Figs. 3 and 4. In the first phase of redevelopment at the end of the fourteenth century, the church was reinforced by an outside solid brick wall. In the second phase, in the fifteenth century, the church was enlarged by an additional five chapels, which have been roofed vaulted using bricks. In the Baroque period after the year 1708, the windows were enlarged, and some doorways were bricked up. Finally, in the fourth phase, after 1819, a brick tower was built. At that time, the new owner tore down the walls between the main nave and chapels located on the north wall, creating a new second nave. The church was used for its intended purpose until 1945. Currently, the church is not in use, but it is included in the register of monuments.
More information related to the history of church can be found in a study undertaken by the staff of the Institute of Monuments and Conservation of the Nicolaus Copernicus University in Toruń, Poland [19]. Nowadays, the church is a two-nave structure. All the walls are brick masonry with lime mortar. The northern wall is preserved rather well, and only small vertical and diagonal cracks are observed. The aisle located at the north wall shall be shuttered by brick vaults. In the western wall, a sloping cracking can be seen above the arch of a window (see Fig. 14). The southern wall is comprised theoretically of two layers (see Figs. 5 and 6). The inner part of the wall is a kind of half-timbered wall with a timber frame and bricks (see Fig. 6b). However, the wood components are very damaged by fungi and insects (see Fig. 6). The southern external wall is made of brick. The wall is very cracked, and the cracks are mainly vertical (see Fig. 5). On the southern wall, one can see some local deformations. Some of the bricks are deeply corroded.
In the XIX century, two brick buttresses were built in the western corner of the south wall and in the east wall (see Fig. 4). In 2009, during the renovation, three columns of concrete were built into the southern inner wall (see Figs. 3 and 4). These three pillars are not founded on any foundation. The eastern wall is also comprised of two layers, and some vertical cracks of significant width are observed on it (see Figs. 2 and 15).
The roof of the church is made of wood (see Figs. 5-7). The main destroyed roof timber truss is considered in analysis here (see Fig. 7). The roof of the tower is omitted. The material of the truss is modelled in a simplified manner as elastic and isotropic even if wood is known as an anisotropic material. Because the church is a structure that was built centuries ago, the material may also have properties that significantly differ from the contemporary wood. Thus, the wood material parameters have been identified on the basis of a series of standard laboratory tests [20] performed on some specimens that were available in the building, and an average value was taken to the analysis.

Laser scanning and geometrical data preparation
The numerical models of the structures were created based on a standard architectural inventory with support from laser scanning (see Fig. 8). This technique and other similar techniques are not new and have been used for similar purposes, as presented in detail in [9,21]. A Leica ScanStation P20 3D laser scanner was used for the inventory. The data acquisition was based on the resolution  necessary to correctly document the structural elements of the church and to point to the damage of the structure since a preliminary study showed some cracks and incorrect structural geometry that was to be considered in the numerical model. A cloud of points has been obtained. Its alignment was carried out by means of common targets between successive scans combined with a system of adjustments between surfaces based on the iterative closest point (ICP) algorithm [16]. Afterwards, cleaning and segmentation were performed. Isolated points were eliminated, the mesh was regularized, and the cloud of points was divided into two parts referring to the church structure and its surroundings. Next, triangulation and texturing of the surface model were performed. A 2D Delaunay flat triangulation was applied to each surface of the structure with respect to the visualization plane. The texturing allowed for obtaining a photorealistic 3D model of the church, from which some ortho-images were generated. The 3D laser scanner data were processed in the 3D point cloud software suite, which consists of Leica Cyclone stand-alone software and Leica CloudWorx plug-in tools for CAD systems.

Geotechnical investigations
The church is resting on the foundation approximately 1-1.5 meters below ground level (m bgl). As a part of field research, the locations of test points have been specified. The soil samples were retrieved from a depth of 4 m. Within the testing soil classification, water content, specific weight, cohesion, friction and soil consolidation properties were performed. On that basis, two main geotechnical layers were specified. The first layer of the soil under the church foundation is a humus clay that goes down to 2.7 m bgl and then up to 4 m bgl, where the second layer of dust clay is placed. Groundwater occurs below the foundation of the building (1.7-1.8 m bgl) in the form of filtration. The geotechnical investigation results indicate that the soil conditions are rather good for supporting the    Fig. 9). In this paper, the support substrate on which the structure rests is treated in a simplified manner as a homogeneous and elastic material, without modelling the details of the support. The constrained modulus from oedometer testing is E ode = 29700 kPa [22]. Assuming that the Poisson's ratio of the tested soil substrate is v = 0.25, we can determine the Young's modulus of the soil according to the formula = [23]. Furthermore, it is assumed that the average value of the soil mass is 1775 kg/m 3 .

Homogenization
The masonry church wall material modulus of elasticity was calculated using a homogenization algorithm described in papers [24,25]. In this method, an elementary cell is analyzed, where based on strain compatibility and stress balance equations of its components, the basic material properties of the entire element can be calculated. The cell is comprised of a brick and mortar attached along its two edges (see Fig. 10). The method assumes that the width and height of an elementary cell is greater than its thickness, which makes it possible to assume a plane stress state in the wall; the layout of bricks and joints in the wall is regular, and the stress gradient in the construction is small, with no concentrations.
The horizontal stress balance equations σ x can be presented in the form where ε ij represents strains in the i-th direction of the j-th component, and l j is the width of the j-th component. The material parameters of the brick wall components are presented in Table 1. The homogenized wall characteristics in two directions perpendicular to each other taken to the analysis are shown in Table 2.
Additionally, the compressive strength of a masonry wall f k was determined based on empirical relationship (see e.g. [26]) where χ is a coefficient accounting for the characteristics of a historic wall, f b is the compressive strength of bricks, and f m is the compressive strength of mortar. Based on compressive stress parameters the tensile stress parameters were determined by means of parametric analysis according to the crack image (10% of the maximum compressive stresses were appropriately adopted for concrete, 4.5% for materials in phases  At this point it should be noted that masonry in its regular texture is always orthotropic. Consequently, the proposed numerical modelling method of masonry due to S. Oller and his co-workers is simplified but very useful in large masonry structure modelling in   nonlinear range. Additional important information on masonry structure orthotropic or isotropic modelling can be found in [27][28][29][30][31]. At the point it should be emphasized that homogenization procedure should be used rather only in the elastic and isotropic range of material behaviour. On the other hand, the use of the isotropic damage models is in some sense acceptable due to the fact that, when regularity of the structure is more and more lost, the masonry structure becomes progressively isotropic. Furthermore, it should be noted that in the inelastic range, at failure analysis, it is still not possible to account for the orthotropic materials with a commercial code (e.g. ABAQUS). The reader will find a detailed description of this issue in the papers [32][33][34].

Finite element modelling and numerical analysis
The model of the church in Gnojewo is based on 3D geometry obtained from a laser scanned and post-processed cloud of points imported into the computational system (see, e.g., [14]) (see Fig. 3). The numerical model of the church and the mechanical analysis were performed using ABAQUS code. The whole numerical model (see Fig. 13) is built of 4-node linear tetrahedral volumetric finite elements with one integration point and three degrees of freedom in each node-type C3D4 [35]. The total number of finite elements in all the analysis cases equals approximately 3.5 million. The FE mesh density was determined in a parametric study based on convergence analysis of the cracked form. The numerical model is loaded only by its self-weight. The geometrically nonlinear analysis was performed using a dynamic explicit procedure [36], as it is computationally efficient for the analysis of large models with relatively short dynamic response times and for the analysis of extremely discontinuous events or processes. The dynamic explicit analysis can also be used to perform a quasistatic analysis. In the proposed numerical solution (quasi-static technique), the load is applied smoothly; consequently, slow deformation triggers a low strain rate. This type of analysis allows for using a consistent, large-deformation theory. In view of the above, the proposed method of analysis is highly effective for studying the damage and degradation of masonry structures. Due to this approach, one can track very carefully the propagation of gradually occurring damage (like cracking and crushing) under an increasing load value, regardless of the type of loads. The ABAQUS package apparently offers several other methods performing dynamic analysis of nonlinear problems, such as implicit dynamic analysis. However, implicit dynamic analysis is not as effective as the explicit approach. The authors' experience shows that while modelling large complex nonlinear problems, like damage of masonry structures, a solution convergence problem often emerges at the time of structural degradation start while applying implicit dynamic analysis.
It should be noted that such an analysis may sometimes be time-consuming. One disadvantage in an explicit dynamic analysis is that a few very small elements will force the entire model to be integrated with a small time increment. One can use variable integration steps or subcycling methods to reduce this problem [36]. In all the numerical analysis cases, the computational time for a single simulation was no longer than 24 h. The analysis was carried out on a supercomputer using 24 computing cores at Tricity Academic Supercomputer & networK center of informatics (http://task.gda.pl).
The support of the structure is modelled as an elastic foundation, without modelling the details of the support. The soil under the structure is modelled using the same C3D4 type of element. In the case of dynamic analysis, unfortunately, one cannot use the classic foundation element available in ABAQUS [36] (using * ELASTIC FOUNDATION command). The necessary soil volume in numerical models can be identified individually, e.g. according to crack growth of the masonry. That may be implemented due to the crack pattern change accompanying the change of soil volume in a modelled structure. Next, given a certain soil volume the main image of cracks remains constant and independent of soil volume. As the volume of soil increases, structural settlement stabilizes. Numerical analysis considers soil weakening under the south-eastern corner of the church. Considering the material behaviour of the specific parts of the structure, the soil (basic and degraded), stone foundations and timber roof are modelled as linearly elastic (see Table 3), whereas the masonry walls (see Table 3 and Fig. 11) and concrete columns (see Fig. 12) are modelled as continuum plastic with damage and stiffness degradation in the sense of continuum damage mechanics. Plasticity-based damage constitutive modelling is widely used to describe the nonlinear properties of masonry structures and concrete. Such a model was first proposed for concrete analysis [37][38][39] and reinforced concrete (RC) structures [40] but also can be used to represent other brittle material behaviours such as for ceramic masonry structures [41,18,42]. The applied material model is described extensively in [36]. It assumes that the main two failure mechanisms are tensile cracking and compressive crushing of the masonry or concrete. The evolution of the yield (or failure) surface is controlled by two hardening variables ∼ ε t pl , which is a tensile equivalent plastic strain, and ∼ ε c pl , which is a compressive plastic strain. Both refer to failure mechanisms under tension and compression loading. It assumes that uniaxial tensile and compressive response of the brittle-like material is characterized by damage plasticity, as shown in Figs. 11 and 12. Under uniaxial tension, the stress-strain response follows a linear elastic relationship until the value of the failure stress σ t0 is reached. The failure stress corresponds to the onset of micro-cracking in the brittle material. Beyond the failure stress, the formation of micro-cracks is represented macroscopically with a softening stress-strain response, which induces strain localization in the concrete structure. Under uniaxial compression, the response is linear until the value of initial yield, σ c0 . In the plastic regime, the response is typically characterized by stress hardening followed by strain softening beyond the ultimate stress. It is assumed that the stress-strain curves can be converted into stress-plastic strain curves. If E 0 is the initial (undamaged) elastic stiffness of the material, the stress-strain curves under uniaxial tension and compression  loading are, respectively, The effective tensile and compressive cohesion stresses are The concrete damage plasticity model assumes that the reduction of the elastic modulus is given in terms of a scalar degradation variable d as = − E d E (1 ) 0 . The stress-strain relationship in the multiaxial conditions takes the following form of the scalar damage elasticity equation , where D el 0 is the initial (undamaged) elasticity matrix, d is the stiffness degradation value, and ⩽ ≤ d 0 1 . The effective stress is defined as = − σ D ε ε : ( ) The model uses the yield function of Lubliner et al. [38], with the modifications by Lee and Fenves [39], to account for different evolutions of strength under tension and compression. The evolution of the yield surface is controlled by the hardening variables, ∼ ε t pl and ∼ ε c pl . In terms of effective stresses, the yield function takes the form= where Here,σ max is the maximum principal effective stress, q and p are stress invariants, σ σ / b c 0 0 is the ratio of initial equibiaxial compressive yield stress to initial uniaxial compressive yield stress, K c is the ratio of the second stress invariant on the tensile meridian to that on the compressive meridian at initial yield for any given value of the pressure invariant such that the maximum principal stress is negative. It must satisfy the condition pl is the effective tensile cohesion stress, and ∼ σ ε ( ) c c pl is the effective compressive cohesion stress. To avoid the convergence problems of models with softening behaviour and stiffness degradation, a viscoplastic Duvaut-Lions regularization of the constitutive equation is used. That implies using a small viscosity parameter. This technique causes the consistent tangent stiffness of the softening material to become positive for sufficiently small time increments. The phases (Phases: 1/2, 3, 4 and 5) of the elaboration of the church structure are considered using different material properties of the masonry appropriate to their building times (see Table 3 and Figs. 11 and 12). Key parameters such as modulus of elasticity, Poisson's ratio, density, stress-strain compression/tension curves and damage variables adopted for masonry walls and concrete columns have been defined on the basis of following a review of the existing literature [43,17,44,45] and our own research [22,46]. Other material model damage plasticity parameters such as deviatoric out-of-roundness, biaxial strength ratio, vertex rounding (eccentricity), which are determined by the initial shape of the failure surface or the dilation angle that specifies the direction of the plastic flow, and viscosity, which affects the numerical stability of the constitutive law, adopted in the numerical models have been defined in an approximate way on the basis of data available in the literature [40][41][42]47].

Results and conclusions
The main results of the numerical analysis are shown in Figs. 14-16. All the main cracks of the real structure of the church are reflected by the numerical results. Fig. 14 shows the main cracks of the west-south walls in the real structure and the resulting Tension damage d t is expressed by means of a tension damage factor, i.e. when using ABAQUS and the DAMAGET function [36]. This parameter is a scalar variable that describes stiffness degradation and takes values from 0 (undamaged material) to 1 (material completely destroyed). Damage to the structure due to crushing (ABAQUS function DAMAGEC) is practically not observed. Thus the problem of structure crushing has not been described in the paper in detail.
This paper presents a numerical analysis of a real historic church structure using FEM based on the laser scans of its geometry. Particular attention has been given to identify the causes of cracking and crushing of masonry walls. A non-linear dynamic explicit analysis has been performed. A crack damage-plasticity isotropic material was adopted to model the masonry, which typically exhibits very low tensile strength in tension. Parameters for the materials used in the masonry modelling on the basis of a homogenization procedure were proposed. In the numerical analysis, the interactions between the church structure, foundations and soil were taken into account.
From the overall analysis of the results obtained, the following conclusions can be drawn.
• The first and main cause of the cracking damage is a result of improper techniques used in the successive reconstructions of the church. Due to frequent church rebuilding in many parts of the structure, ceramic materials (bricks) with different strength parameters are combined (see e.g. Fig. 16), which is unfortunately contrary to the known principle that "No man also seweth a piece of new cloth on an old garment: else the new piece that filled it up taketh away from the old, and the rent is made worse. And no man putteth new wine into old bottles: else the new wine doth burst the bottles, and the wine is spilled, and the bottles • The second but also very important cause of the masonry cracking is varied settlement of soil under the church structure (see Fig. 4). The numerical study results clearly indicate that the soil conditions, although seemingly good enough for supporting the structure, unfortunately are the cause of cracks in masonry. The geotechnical investigation carried out did not show problems with the bearing capacity of the soil merely because of the inadequate number of geotechnical test drillings (see Fig. 9). An additional geotechnical drilling should be conducted at the south-east corner of the church, where the soil seems to be weaker. The soil weakness can be confirmed by the reconstruction of the church in the 18th century (see Figs. 4 and 15 -phase 4) because the buttresses probably built in at that time were supposed to prevent the structure from collapsing. Moreover, a sacristy built later at the same corner made the situation even worse because it caused the overloading of the soil and the east masonry wall cracking (see Fig. 15).
• It should be emphasized that the soil model is an important element of the whole analysis. In particular, modelling of massive masonry structures without taking into account the interaction of structure and soil can lead to incorrect assessment of the building's behaviour. Inclusion of the soil model in the presented analysis, even as a simple isotropic continuum, allowed to obtain simulation results corresponding to real structural damage, in particular the place of the appearance of cracks (see Figs. [17][18][19][20]. This would not be possible only by applying boundary conditions for the foundations of the church, as it is often found in the literature. • The proposed FEM model allows for the rather precise identification of the causes and effects of cracking walls in a qualitative sense but not very well in a quantitative sense. The main reason for this is that a classical FE analysis within the continuum mechanics still used in commercial software (e.g. ABAQUS) does not properly describe the problem of localized zones. Thus, numerical models suffer from a pathological mesh-dependency. This problem is described very well by Tejchman and Bobiński, in [48]. As mentioned in that paper, some researchers try to solve the problem of localization zones using 2D instead of 3D finite elements when modelling masonry structures.
• Furthermore, it is desirable to perform various numerical investigation to analyse variation of structural performance in different phases of its reconstruction. This is a complex task due to the lack of sufficient data concerning the object. Also the sensitivity analysis of selected material parameters to assess their influence on numerical results could be interesting. Nevertheless, these issues were deliberately neglected due to the problem complexity and the volume of the paper.

Declaration of Competing Interest
The authors declared that there is no conflict of interest.