A hybrid macro-modelling strategy with multi-objective calibration for accurate simulation of multi-ring masonry arches and bridges

This paper presents an efﬁcient hybrid continuum-discrete macro-modelling strategy with an enhanced multiscale calibration procedure for realistic simulations of brick/block-masonry bridges. The response of these structures is affected by the intrinsic nonlinearity of the masonry material, which in turn depends upon the mechanical properties of units and mortar joints and the bond characteristics. Finite element approaches based upon homogenised representations are widely employed to assess the nonlinear behaviour up to collapse, as they are generally associated with a limited computational demand. However, such models require an accurate calibration of model material parameters to properly allow for masonry bond. According to the proposed approach, the macroscale material parameters are determined by an advanced multi-objective strategy with genetic algorithms from the results of mesoscale ‘‘virtual” tests through the minimisation of appropriate functionals of the scale transition error. The developed continuum-discrete ﬁnite element macroscale description and the calibration procedure are applied to simulate the nonlinear behaviour up to collapse of multi-ring arch-bridge specimens focusing on the 2D planar response. The results obtained are compared to those achieved using detailed mesoscale models conﬁrming the effectiveness and accuracy of the proposed approach for realistic nonlinear simulations of masonry arch bridges.


a b s t r a c t
This paper presents an efficient hybrid continuum-discrete macro-modelling strategy with an enhanced multiscale calibration procedure for realistic simulations of brick/block-masonry bridges. The response of these structures is affected by the intrinsic nonlinearity of the masonry material, which in turn depends upon the mechanical properties of units and mortar joints and the bond characteristics. Finite element approaches based upon homogenised representations are widely employed to assess the nonlinear behaviour up to collapse, as they are generally associated with a limited computational demand. However, such models require an accurate calibration of model material parameters to properly allow for masonry bond. According to the proposed approach, the macroscale material parameters are determined by an advanced multi-objective strategy with genetic algorithms from the results of mesoscale ''virtual" tests through the minimisation of appropriate functionals of the scale transition error. The developed continuum-discrete finite element macroscale description and the calibration procedure are applied to simulate the nonlinear behaviour up to collapse of multi-ring arch-bridge specimens focusing on the 2D planar response. The results obtained are compared to those achieved using detailed mesoscale models confirming the effectiveness and accuracy of the proposed approach for realistic nonlinear simulations of masonry arch bridges.

Introduction
Old masonry arch bridges belong to the cultural and architectural heritage and still play a critical role within railway and roadway networks in Europe and worldwide. These structures were built following empirical rules and were not designed to resist current traffic loading and the loads induced by extreme events, such as earthquakes. An accurate assessment of the ultimate performance of these complex structural systems represents a crucial step to prevent future failures and preserve such historical structures for the next generations.
Masonry arch barrels are the key structural components of masonry arch bridges. Their nonlinear behaviour is strongly influenced by the mechanical properties of the two constituents, masonry units and mortar joints, and their arrangement to form the brick/blockwork of the arch (i.e. masonry bond). Two main categories of masonry arch bridges can be identified: stone masonry and brick masonry bridges [1]. In the first group, the arches are built from large voussoirs organised in a single arch ring. Conversely, in the case of brick masonry bridges, a multi-ring arrangement is usually utilised, where the number of rings depends on the span length of the arch. The rings are typically bonded together using the stretcher method, where the connection between adjoining rings is guaranteed by continuous mortar joints. To date, numerous laboratory and in-situ tests have been performed to investigate the failure mechanisms of masonry arches and bridges, considering also the influence of backfill, under monotonic and cyclic loading conditions [2]. Specific studies on multi-ring arches showed how ring separation and shear sliding generally affect the ultimate strength and failure mode [3][4][5][6], where weak circumferential mortar joints have been found to lead to an ultimate strength reduction of about 30% for short spans and up to 70% in the case of longer span arches.
In previous research, different numerical strategies have been proposed to simulate the nonlinear behaviour of masonry arches and bridges [2]. Generally, approaches based on limit analysis principles can be effectively used to estimate the ultimate load capacity [7][8][9]. However, such strategies do not provide information about https://doi.org/10.1016/j.compstruc.2022.106769 0045-7949/Ó 2022 Elsevier Ltd. All rights reserved. the nonlinear response before collapse, and they are often based upon crude assumptions, e.g. the representation of masonry as a no-tension material, which may lead to underestimating the ultimate resistance of masonry arches. Previous studies also comprised simplified 2D finite element (FE) limit-analysis descriptions to simulate the arch-backfill interaction [10,11] and 3D nonlinear FE strategies with elasto-plastic solid elements [12][13][14], where masonry is assumed as a homogeneous isotropic material disregarding its anisotropic nature. Isotropic modelling approaches are widely employed in engineering practice due to their computational efficiency, especially for the analysis of large bridges. However, they may lead to an unrealistic representations of typical failure modes not directly associated with longitudinal bending. Also, their application to masonry may require complex calibration procedures to account even in a simplified way for its anisotropic nature, as shown in [15] with reference to masonry walls. Furthermore, the use of more complex damage/plasticity orthotropic models [16][17][18][19][20], based on damage and/or plasticity formulations where tensile, shear and compressive failure mechanisms are described, are still largely applied for research and not yet considered for the practical assessment of realistic structures.
More recent numerical models for masonry arched structures and bridges include the micro-model strategy proposed by Milani et al. [21] using triangular rigid elements and nonlinear links, the discrete macro-element method (DMEM) [22][23][24] and the distinct element method (DEM) [25,26]. A detailed 3D mesoscale modelling strategy for masonry arch bridges has been developed at Imperial College London [27,28], which is used as the reference solution for the calibration of the proposed macroscale approach hereinafter. According to this strategy, the masonry parts of the bridge are simulated by using linear solid elements and 2D nonlinear interface elements to explicitly allow for the masonry bond [29]. The backfill is modelled by elasto-plastic solid elements, and the connection between the masonry components and the backfill is represented through nonlinear interfaces allowing for the actual frictional interaction. This approach generally leads to accurate response predictions, including under extreme loading, but it is associated with significant computational cost which can hinder its use for the practical assessment of real large structures.
With the aim of achieving a suitable compromise between accuracy and efficiency, this paper proposes a hybrid continuumdiscrete macroscale description for multi-ring masonry arches and masonry arches bridges. Elasto-plastic-damage continuum solid elements interacting with 2D nonlinear interfaces are employed to model a masonry arch, although, unlike the mesoscale strategy, mesh discretisation is not directly related to the dimensions of units and mortar joints. The damage-plasticity model proposed in [32] and a multi-surface cohesive-frictional model [15,30,32] are employed for solid and interface elements, respectively. Furthermore, an innovative multi-objective optimisation procedure, based on virtual tests developed adopting detailed mesoscale descriptions, is put forward and applied to evaluate the mechanical parameters of the hybrid model.
The proposed modelling strategy with the advanced calibration procedure is validated against mesoscale simulations considering multi-ring arches and masonry arch-backfill specimens with different geometrical and mechanical properties. The numerical results confirm the accuracy and high efficiency of the developed hybrid approach, which can be used for practical and accurate assessment of realistic masonry arch bridges.

The hybrid macro-modelling approach
In the proposed FE modelling strategy, the arch is discretised using a regular mesh of nonlinear continuum 20-noded solid ele-ments. In addition, 2D nonlinear zero-thickness interface elements are arranged along the circumferential mid-thickness surface of the arch to simulate damage associated with potential ring sliding and separation. In the simplest case where only one circumferential layer of interfaces is considered (Fig. 1), each interface lumps the linear deformability and non-linear behaviour of n-1 ring joints, with n being the number of rings of the physical arch. Importantly, the characteristics of the FE mesh with solid elements are not directly linked to the masonry bond. Thus, an arbitrary number of solid elements can be employed along the length of the arch, according to the desired level of response detail, but at least two solid elements should be arranged along the thickness of the arch to accommodate the mid-thickness nonlinear interfaces. The accuracy due to different discretisation along the circumferential direction is explored in the numerical applications described in the following sections.

The 3D damage-plasticity model
In the macroscale representation implemented in ADAPTIC [31], the isotropic plastic-damage material model presented in [15] is used for the 20-noded solid elements. A standard decomposition of total strains (e) in elastic (e e ) and plastic (e p ) components is considered, and the stress tensor (r) is obtained from the effective stress tensor ( r) and a scalar damage variable dð r; j t ; j c Þ. The latter variable depends on the stress state and two historical variables (j t , j c ) representing the evolution of plastic strains in tension and in compression. The material relationship is expressed analytically by: where E 0 is the initial fourth-order isotropic elastic tensor. The local plastic problem is solved at each integration point of the domain to evaluate the effective stress, adopting a nonassociated elasto-plastic constitutive law with Drucker-Pragerlike plastic flow potential, according to the approach proposed in [32].
The plastic behaviour is governed by the evolution of the yield surface: ð Þwith r i principal effective stress; 2 .
-f c j c ð Þ, f t j t ð Þ effective strength in compression and tension, respectively; -K c ratio of the second stress invariant on the tensile meridian to that on the compressive meridian at initial yield; -f b0 ratio between biaxial and uniaxial compressive strength.
To improve the computational robustness, both tensile and compressive strengths,  Fig. 2.
The global damage variable is obtained as: , scalar parameter ranging from 0 (all principal stresses are negative) to 1 (all principal stresses are positive) expressing the state stress; -w t ; w c are parameters governing the stiffness recovery from compression to tension and vice versa.
Since softening behaviour may lead to mesh sensitivity, a fracture-energy approach has been adopted to maintain objectivity in the results. In particular, the stress-strain constitutive relationship is defined at element level starting from a stress-crack opening curve based on fracture energy, assumed as material parameter and a characteristic length evaluated as a function of the element volume.
The model has been extensively used to simulate the mechanical behaviour of concrete [32,33] and masonry [34][35][36]. Some inherent model characteristics, however, hinder its use to represent specific shear failure modes typical of multi-ring masonry arches. More specifically, the adopted damage-plasticity continuum description does not enable the definition of the shear strength independently from the tension and compression strengths. It can be seen by applying Eq. (2) assuming a pure shear to the yield function: and, since r r ð Þ ¼ 0:5 for pure shear, the damage parameter becomes: Assuming that damage in compression has not developed, d c j c ð Þ ¼ 0, the equivalent damage parameter becomes: From Eq. (8), some considerations can be made on the shear behaviour of the model. Firstly, the initial shear strength (Eq. (5) with j ¼ 0) is governed by the initial compression and tension strengths and the parameter f b0 , which in practice is always in the range 1.12-1. 16. This confirms that it is not possible to define a specific shear strength independent from tension and compression strengths, as for instance a shear strength relating to the sliding of mortar joints, which is a typical shear failure mode for multiring masonry arches. A workaround to have some freedom in the definition of initial shear strength could be to calibrate f c 0 ð Þ ¼ f c0  appropriately and independently from the observed compressive behaviour, while both f t0 and f c;max would still be determined based on their specific failure modes. The second consideration is that the evolution of nominal shear strength depends on the parameter w c (see Eq. (8)) which is defined based on the expected cyclic response (stiffness recovery). A typical value, leading to complete stiffness recovery from tension to compression, is w c ¼ 1:0 [37,38]. Inserting this in Eq. (7), the expression for damage in pure shear in the absence of compression damage is obtained d j t ð Þ ¼ 0:5d t j t ð Þ. The conclusion is that the evolution of nominal shear strength is completely governed by damage in tension, without the possibility for specifying an alternative more realistic constitutive relationship.
Finally, it is worth mentioning that the macroscale damageplasticity continuum representation is not capable of distinguishing failure due to shear parallel to the mortar bed joint s zx from that due to shear orthogonal to the mortar bed joint s xz , as in the Cauchy solid these two stresses are equal, and the yield surface cannot consider separate contributions. In reality, while the former failure mode is governed by sliding of the units on the weak planes represented by the mortar joints, the latter is governed by the internal rotation of bricks depending on their geometric shape ratio and brick interlocking, as schematically shown in (Fig. 3) [39]. To allow for these different phenomena enriched continuum representations, e.g. Cosserat continuum [39], would need to be employed.
To overcome these intrinsic limitations of typical continuum damage-plasticity constitutive models, an alternative hybrid macroscale representation is proposed, in which, as outlined before, shear sliding along the continuous circumferential mortar joints of multi-ring arches is described by introducing nonlinear interfaces whose material characteristics are defined based on the calibration strategy described in Section 3.

Constitutive model for nonlinear interface elements
2D 16-noded interface elements [29] are employed for the midthickness circumferential interfaces using the plasticity-damage constitutive model proposed in [30]. According to this description, interface tractions and relative displacements describing the static and kinematics of the element, are composed of a normal component in the direction orthogonal to the interface and two shear components on the plane of the interface. The effective stresses are evaluated at each Gauss point by solving a linear hardening elasto-plastic problem considering multi-surface plasticity. Then, the nominal stresses are obtained by multiplying the effective stresses by the damage matrix D, containing the damage index in tension, shear and compression ranging from 0 (no-damage) to 1 (complete damage).
Similarly to the solid elements, a standard decomposition between elastic and plastic deformations is considered and the concept of effective stress s ¼ K 0 e À e p À Á is introduced, where K 0 ¼ diag k n k t k t f gis the diagonal initial stiffness matrix with k n and k t the normal and shear stiffness, s ¼ ½ r s 1 s 2 , e ¼ ½e c 1 c 2 and e p ¼ ½e p c p1 c p2 the effective stress, the total strains and the plastic strains, respectively. The nominal stresses are evaluated from the effective stress according to: where D represents an anisotropic damage tensor, containing distinct variables for the normal (D n ) and the tangential (D t ) directions. A tri-linear plastic yield domain is considered to simulate the tensile (Mode I), shear (Mode II) and crushing (Mode III) failure modes. Three distinct plastic works, corresponding to each fracture mode rule the evolution of the damage variables. The plastic yield domain (Fig. 4) is composed of three surfaces, F t , F c , and F s respectively, associated with the tensile (mode I), compression and shear (mode II) failure modes, as defined by: where f t and f c are the tensile and compression material strengths, / the friction angle and q a linear hardening variable, ranging from 0 (initial value) to the limit value q lim ¼ c With the increase of q the surface F t reduces until becoming a point when q reaches the value q lim . On the other hand, F s increases with the increase of q (Fig. 4). Two associated plastic flows are related to F t and F c , while a plastic potential G s in shear, obtained from F s substituting / to / g , is considered to take into account the effects of masonry dilatancy.
Following the solution of the plastic problem, the damage evolution is evaluated as a function of the three ratios r c ¼ W pc =G c , r t ¼ W pt =G t and r s ¼ W ps =G s where W pc ; W pt and W ps are the plastic works in compression, tension and shear, respectively, and G c ; G t ; G s the corresponding fracture energies. Finally, the nominal stresses are given by Eq. (9). More details on the model formulation can be found in [30].

Calibration procedure
The mechanical calibration of the proposed model requires the determination of several material parameters defining the linear and nonlinear behaviour of the 3D solid elements and the 2D interfaces, as described Sections 2.1 and 2.2, and reported in [15,29,30] For this reason, an objective and robust calibration procedure represents a fundamental step to guarantee the model accuracy and applicability.  This work presents an original multi-objective calibration procedure based on the multiscale approach proposed in [15] which considers the representation of a structure under suitable boundary conditions according to two scales: mesoscale, indicated hereinafter by the superscript m, and macroscale, with the superscript M. The considered setup is called virtual test, and it is assumed there exists a mapping M : X m ! X M between the mesoscale and the macroscale domains. As elaborated subsequently, different to the procedure proposed in [15] which consisted of a singleobjective optimisation algorithm, the newly proposed procedure leads to a multi-objective optimisation problem allowing for a set of optimum solutions (Pareto Front) which improves the robustness and accuracy of the model calibration procedure.
According to the original formulation put forward in [15], stress power equivalence between the two scales is approximately enforced on the entire domain of the virtual test. The stress power equivalence reads: where _ represents the error rate due to the approximations induced by the specific macromodel utilised. Considering pseudostatic stress states, the equality between internal and external work implies: where t are the surface forces on the boundary C, while b are volume forces. Neglecting the contribution of these latter and considering the chain rule of differentiation, Eq. (12) finally reads: Eq. (13) represents the error rate at time t due to the scale transition. In [15], a global non-negative monotonically increasing error function was defined: The extension of the original procedure, proposed in this paper, consists of partitioning the error defined as in Eq. (11) or in Eq. (13) as: The contributions _ 1 , _ 2 , . . . respectively refer to a volume partitioning in Eq. (15), with X Mjm 1 þ X Mjm 2 þ Á Á Á ¼ X Mjm , or load-based partitioning in Eq. (16). For the sake of simplicity, in Eq. (15) and (16) it is assumed that there is not any modification of volumes and surfaces in the scale transition, i.e., @C m In this case several error functions can be defined as: The solution of the calibration procedure is given by the solution of the multi-objective minimisation problem: The error partitioning defined in Eqs. (15) or (16) has two consequences. The first consequence is that it allows defining the granularity of the homogenisation, avoiding the possible error compensations given by different parts of the structure. For instance, if the nonlinearities in the mesoscale model are concentrated in one small region of the domain, it is possible to use volume partitioning in Eq. (15) to focus the calibration of the parameters governing the nonlinear behaviour of the macroscale representation in that region, while controlling the elastic parameters by matching the response in the remaining domain. The second consequence is that the calibration problem is turned into a multi-objective optimisation problem, in contrast to the original formulation [15] which was a single-objective optimisation procedure. As shown in [40,41], using multiple objectives in a calibration problem may strongly increase the robustness of the procedure. In the numerical applications reported in Sections 4 and 5, two partitions of the global error are considered to simplify the interpretation of the optimisation results. However, the use of a larger number of partitions may be considered.
The multi-objective optimisation problem is solved by means of a Non-Dominated Sorting Genetic Algorithm [42], implemented in TOSCA-TS software [43]. The optimum is given by the Pareto Front (PF), which represents the set of non-dominated solutions. A careful investigation on the features of the Pareto Front may highlight possible inconsistencies of the model to calibrate [40] and represents a key part of the calibration procedure towards the definition of the most representative solutions and a significant improvement of the original procedure presented in [15].
Finally, it is worth noting that the calibration strategy considers the evolution of the stress power over time, and thus it cannot properly allow for the additional work contribution due to initial loading. Thus, it is preferable to avoid initial loading in the virtual test. However, this does not limit the applicability of the procedure, as multiple loads with independent loading paths can be introduced without any modifications in the methodology.
The proposed strategy enables an objective evaluation of the macroscale model parameters given the masonry mesoscale properties which can be obtained directly from simple in-situ or laboratory tests performed on masonry units and mortar (or tests performed on small assemblages of units such as triplets), following consolidated methodologies already reported in the literature [44].

Numerical simulations of medium span masonry arches and bridges
Two medium-span masonry arch specimens, one interacting with backfill as found in typical masonry arch bridges, are analysed by means of detailed mesoscale models [29]. The results of the mesoscale analyses are then used as reference solutions to highlight some limits of a typical continuum macroscale description, and to investigate the improved accuracy guaranteed by the proposed continuum-discrete hybrid representation for multi-ring masonry arches.

Masonry arch and bridge specimens
The first specimen (Fig. 5a) consists of a 5 m span three-ring brick-masonry arch. The arch is characterised by 1250 mm rise, 330 mm thickness and 675 mm width. Adjacent rings with 215 Â 102.5 Â 65 mm 3 bricks are connected according to the stretcher method by continuous circumferential 10 mm thick mortar joints [45]. The second specimen (confined arch) comprises a brick-masonry arch with the same geometrical characteristics of the bare arch interacting with backfill material, which, extends 2460 mm horizontally from the two supports of the arch and 300 mm vertically above the crown, according to the experimental layout considered in [4], (Fig. 5b). Full supports are assumed at the base of the arch and the backfill, while simple supports against the horizontal longitudinal displacements are applied on the two vertical sides of the backfill. Moreover, the horizontal transverse displacements on the two lateral faces of the arch and backfill are restrained to represent a plane strain condition (Fig. 5b).

Mesoscale simulations
In the numerical mesoscale description, 20-noded elastic solid elements are used to simulate the brick units and 16-noded interfaces [29] are employed to represent both the radial and the circumferential mortar joints. As the focus is on the 2D response, a mesh with only one element along the representative 1 m width of the arch specimens is considered. The mesoscale description of the arch requires 240 3D solid elements, 403 2D interface elements and 6453 nodes to which correspond 19,359 DOFs. The backfill is modelled adopting a FE mesh with 15-noded tetrahedral elements. Finally, nonlinear interface elements are utilised to model the physical interface connecting the arch to the backfill where the generally nonmatching meshes of the two domains are coupled using a mesh tying method [50].
Two masonry types have been considered in the analyses: a strong masonry to represent modern good quality brickwork, and a weak masonry to represent historical masonry [46]. The mesoscale mechanical parameters are reported in Tables 1 and 2. These parameters have been selected based on previous numerical studies where the adopted mesoscale description for masonry arches and bridges was validated against experimental tests. In particular, the parameters of the strong masonry have been adopted in [27,28] to reproduce the response of a two-ring 3 m span arch, while the parameters for weak masonry were adopted in [47] in further validations against physical experiments.
Following [28], an elasto-plastic material model with a modified Drucker-Prager yield criterion is employed for the backfill, assuming a Young's modulus E b ¼ 500 MPa, a cohesion c b = 0.001-MPa, a friction and a dilatancy coefficient tan/ b ¼ 0:95 and tanw b ¼ 0:45. The nonlinear interfaces simulating the interaction between the arch and the backfill at the extrados of the arch have tensile strength f fi ¼ 0:002 MPa, cohesion c fi ¼ 0:0029 MPa, friction coefficient tan/ fi ¼ 0:6 and zero dilatancy.
In the numerical simulations of the bare arch, two initial vertical forces F 0 = 22.5 kN are applied at the quarter and three-quarter span and then maintained constant during the subsequent loading stage, when a vertical force F is applied at quarter span and monotonically increased up to collapse. Both forces F 0 and F are uniformly distributed on a patch area of 210 Â 675 mm 2 .
When the arch interacts with the backfill, the initial load corresponds to the weight of the arch and the backfill both with a specific weight of 22 kN/m 3 , while the force F is applied on the top surface of the backfill on a patch area of 400 Â 675 mm 2 centred at the quarter span of the arch (Fig. 5b). To improve the numerical stability, nonlinear dynamic analysis is performed by imposing an initial velocity of 0.1 mm/s at the loaded nodes, which is maintained constant during the simulation up to collapse. Zero viscous damping is considered in the analyses. Fig. 6 shows the load-displacement curves of the bare arch (Fig. 6a) and the arch interacting with backfill ( Fig. 6b), where the force F is plotted against the vertical deflection at the quarter span of the arch. In Fig. 6a, the ultimate load estimated by limit analysis, based on the classic Heyman's hypotheses [48], is reported for comparison.
A significant influence of the masonry typology on the global response is observed both in the case of the bare arch, where the ultimate load ranges from 31kN to 66kN, and for the arch with backfill, where the peak force varies from 74kN to 137kN. As  expected, the initial stiffness of the bare arch is significantly affected by the masonry characteristics. Conversely, the confined arch shows almost the same initial stiffness for weak and strong masonry. Considering the specific weight of the strong masonry (Table 1), standard limit analysis provides a prediction of the peak-load (39.43 kN) significantly lower compared to the strongmasonry model due to the hypothesis of no-tension material. At the same time, it provides an overestimated peak-load compared to the weak-masonry model because it neglects the sliding between the rings. The failure mechanisms and the equivalent von-Mises stress contours of the two specimens, with both masonry typologies, are reported in Fig. 7. Finally, Fig. 8 shows the tensile damage contours at the last step of the analysis obtained by the different models. The failure mechanism of the models with weak masonry is characterised by shear sliding along the circumferential interfaces, mainly concentrated in the zone between the left support of the arch and the loading area at quarter span, and close to the three-quarter span of the arch. This mechanism prevents the activation of flexural plastic hinges. In the case of strong masonry, flexural failure is observed with the opening of four radial cracks in both the bare arch (Fig. 7) and the arch with backfill (Fig. 7d). In the models with weak masonry, significant damage in the radial joints is observed close to the load. In the Table 2 Interface mechanical parameters adopted in the analyses.

Masonry
kn À kt ½N=mm 3 f t À f c À c ½MPa Gt À Gs À Gc ½N=mm tg/ À tg/ g ½À   models with backfill, large portions at the extrados of the arch are affected by shear-sliding damage, both at the radial and circumferential interfaces. This damage develops also at the frame-backfill interfaces. In the following, the mesoscale solutions are assumed as the baseline results for the assessment of more efficient macroscale models and the proposed hybrid continuum-discrete descriptions for multi-ring arches.

Macroscale simulations
The two specimens presented in Section 4.1 have been analysed by a continuum macroscale description for masonry, using the isotropic damage-plasticity constitutive law described before. Since the model is developed employing quadratic elements, a relatively coarse mesh can be used to improve computational efficiency. More specifically, a mesh with two elements along the thickness of the arch with a length in the circumferential direction approximately equal to half the thickness of the arch has been considered in the numerical simulations. Moreover, as for the mesoscale model, only one solid element is arranged along the 1 m width of the arch. As a result of this, the masonry arch is represented by 80 3D solid elements, 40 2D interface elements and 981 nodes corresponding to 2943 DOFs. It can be observed that the macro-modelling description allows a reduction of 85% of DOFs compared to the mesoscale description demonstrating the potential for considerable reduction in computational demands with the proposed model.
The aim of this investigation is to explore the accuracy and potential limitations of a standard continuum isotropic macroscale approach to predict the response of multi-ring masonry arches, where the model material parameters are calibrated according to two alternative simplified and advanced procedures.

Simplified calibration procedure
In initial macroscale simulations, the material model parameters for masonry have been evaluated through a simplified calibration procedure considering the mesoscale material properties reported in Tables 1 and 2. The macroscopic Young's modulus for the masonry material E has been determined by combining in series the stiffness of brick units with that of the mortar interfaces along the direction of the arch. The tensile strength and fracture energy f t ; G t and the compressive strength f c are assumed coincident to the corresponding values of the mesoscale interfaces. The remaining parameters for the damage plasticity model are assumed equal to standard values used in previous studies for modelling masonry materials.
More specifically: The ratio between initial and maximum compressive strength f y ¼ f c0 f c;max is assumed equal to 0.3 according to [32,33,51]; The dilatancy angle w is taken equal to 35°which is consistent with the value adopted for modelling quasi-brittle material as concrete [32,33] and corresponds approximately to the median of the values (ranging from 10°to 50°) typically used for masonry [34][35][36]; The eccentricity of the plastic flow potential is taken as ¼ 0:1 to improve computational robustness as suggested in [37]; l governing the relative influence of damage and plasticity in tension (l ¼ 0 for fully damage material) is assumed equal to 0.2; The plastic strain at maximum compression stress k c;fc is taken as 0.002 following [51]; The ratio between the plastic strain at damage onset in compression and the plastic strain at maximum compression q c is considered equal to 1.0, as damage is assumed to develop in the softening branch of the stress-strain response.
As noted in Section 2.1, preliminary numerical simulations showed a significant influence of the parameter governing the stiffness recovery in compression w c on the global response of the arch. Thus, two limit values (0, 1) are considered, while parameter w t determining the stiffness recovery in tension is assumed as zero. The complete set of mechanical properties for the continuum macro-modelling description are reported in Table 3.  Fig. 9 shows the load-displacement responses predicted by the macroscale descriptions which are compared against the reference mesoscale curves. In the main, the macromodels predict the initial stiffness of the masonry arches accurately, yet significantly overestimating the peak strengths without providing a realistic representation of the post-peak behaviour as given by the reference mesoscale models. Furthermore, very different macromodel curves are obtained depending on the adopted value for w c . In particular, the largest differences between the mesoscale models and the corresponding macromodels are achieved when w c = 1. It should be noted that this value is recommended by most software implementations [38,39] to model the cyclic response of quasi-brittle materials.
In Fig. 10, the influence of the dilatancy angle on the global response of the weak masonry bare arch is shown. Since this parameter governs the normal plastic deformation due to shear, it is expected that by increasing w the global behaviour becomes more ductile due to the confinement effects exerted by the surrounding elements. Given its high influence on the global beha- Table 3 Macroscopic mechanical parameters resulting from the simplified calibration procedure. [-] Kc   viour, it is apparent that more accurate calibration is needed for such critical parameter. The deformed shapes at failure (at the last step of the analysis) with the tensile damage contour distributions are depicted in Fig. 11. The models with strong masonry exhibit flexural failure (Fig. 11c, d) which is in agreement with the failure mode predicted by the mesoscale models (Fig. 7c, d) and mostly characterised by damage in tension concentrated at the intrados and extrados of the arch corresponding to the opening of plastic hinges. The models with weak masonry show featuring a local shear failure developing underneath the area where the external load is applied and flexural damage at the extrados of the arch at the side opposite to the loaded area (Fig. 11a, b). This is not predicted by the mesoscale model, which shows ring separation at failure (Fig. 7a, b). This main difference confirms the inability of the continuum isotropic damage-plasticity model to represent shear sliding between adjacent rings, which is a characteristic failure mechanism of multiring arches well captured by detailed mesoscale models. Moreover for the arches with strong masonry, the use of the macroscale continuum isotropic model leads to a significant overestimation of the ultimate strength and ductility, where the numerical predictions are affected significantly by some model parameters (e.g. w c and w) which cannot be determined via simplified calibrations.

Advanced calibration procedure
To improve the accuracy of the macroscale predictions, the advanced calibration procedure described in Section 3 has been applied to determine the macromodel material parameters, focusing first on the specimens with weak masonry where the initial macroscale predictions, based on a simplified calibration of the model material parameters, were not in good agreement with the mesoscale results.
The bare arch in Fig. 5a, subjected to two constant initial forces at the quarter span (L/4) and three-quarter span (3/4L), both equal to 16 kN, and to a patch load applied at L/4 and increased up to collapse, is used as the virtual test for the calibration of the model parameters. The specific loading condition ignoring the selfweight contribution of masonry has been chosen to activate both flexural damage and shear sliding between adjacent the rings, thus providing suitable information to the optimisation algorithm.
It should be pointed out that the selection of appropriate virtual tests which should activate the most critical failure modes of the investigated masonry specimens is the critical step for a successful application of the proposed calibration strategy. For instance, in the case of multi-ring arches, the failure mechanism of the virtual test should be characterised by flexural damage, namely the activation of one or more plastic hinges and shear sliding along the rings. However, if it is not possible to identify a virtual test with these characteristics, multiple virtual tests may be utilised, and the multi-objective optimisation procedure should consider error functions for each of these.
Some parametric analyses, not included in the paper for the sake of brevity, have been performed to identify the parameters that affect most significantly the arch response. As a result of these parametric analyses, and to limit the computing time associated with the model calibration, six model parameters are considered as unknowns in the optimisation procedure: the Young modulus (E), the tensile strength (f t ) and fracture energy (G t ), the ratio (f y ) between the yielding and ultimate compression strength, the parameter governing stiffness recovery from tension to compression (w c ) and the angle of dilatancy (w). The range for each parameter is reported in the Table 4. The remaining parameters of the solid elements are assumed equal to the default values in Table 3.
A load-based partitioning strategy is used for the solution of the calibration problem based on two objectives as defined in Eqs (16), (18) with x 1 , x 2 the errors due to the loads at L/4 (F 1 ) and 3/4L (F 2 ), normalised with respect to a reference value with the same units (final squared strain energy, divided by the time interval of the virtual test, [J 2 /s]).
The evaluated PF (Fig. 12) appears discontinuous; the minimum of x 2 (4 Á 10 À3 ) is reached for x 1 ¼ 0:10. Conversely, the minimum of x 1 (6 Á 10 À4 ) corresponds to a much larger value x 2 ¼ 1: 19. This implies that calibrating the response on the force-displacement curve of the variable load at L/4 (minimum x 1 ) may entail significant error in the total energy. In Fig. 12a, the solutions on the PF are shown using the normalised errors Þ ranging from 0 to 1. In the following, the solution corresponding to the minimum global error  Fig. 12b shows the displacement d 1 at L/4 against the load F 1 associated with all the solutions of the PF compared to the response of the mesoscale model. Three families of curves can be observed: the curves that minimise the error related to F 1 , which fit very well the response of the mesoscale response; the curves that minimise x 2 , i.e., the error related to F 2 , which provide a sig-nificant overestimation of the peak-strength of the arch (from 35 kN to 45 kN) and the curves that minimise the global dimensionless error allowing for the two objectives of the optimisation procedure. These solutions and in particular the solution associated with the minimum error provide a good prediction in terms of initial stiffness and peak load, but they show some differences regarding the post-peak stage. In Fig. 12c, the vertical displacement at three-quarter d 2 span is plotted against the load F 1 . It is possible to see that in no case the final uplift shown by the mesoscale model  is attained by the macroscale models of the PF. However, the dashed line, identifying the compromise solution of the multiobjective optimisation, shows a general good agreement with the mesoscale curve. It is important to point out that even when the minimum error solution provides a satisfactory approximation of the load-displacement curves, it may lead to a failure mechanism inconsistent with that predicted by the virtual test. Therefore, a comparison in terms of failure mechanism is crucial to choose the best solution from the PF.
The ultimate deformed shapes corresponding to the minimum error solution and the two ends of the PF, each providing the minimum of one error function, are reported in Fig. 12d. It can be noticed that the solution corresponding to x* 2,min (best fit of the load-displacement at the loaded point) does not indicate shear failure of the arch, leading to a wrong failure mode prediction.
Conversely, the two other solutions predict different shear failure mechanisms. In particular, the minimum error solution (x* min ) represents well the failure mode of the virtual test by predicting circumferential shear deformations, in accordance with the shear sliding among adjacent rings determined by the mesoscale model.
In Table 4, the calibrated parameters for the minimum error solution are shown, together with the range identified by the solutions in its surroundings (x Ã < 2:5x Ã min ). These latter values allow investigating the sensitivity of the calibrated response to each parameter, compared to the initial variation range. Analysing the  results, all parameters seem rather univocally determined as a significant reduction of the ranges is observed, with the only exception of f y . It is interesting to see that the calibrated w c is close to zero, unlikely the typical assumption considered in the simplified calibration.
The results of the model calibration are validated by analysing the confined and bare arches subjected to the loading condition described in the initial mesoscale simulation in Section 4.2.
Moreover, two additional load conditions for the bare arch are investigated in which the arch is subjected to a patch load alternatively applied at the mid span and at one eight span without initial symmetric forces. The load-displacement curves for the weak masonry models are shown in Fig. 13, where the continuum calibrated model is compared against the mesoscale model and the continuum model with simplified calibration procedure for w c = 0.  Generally, the model calibrated by the advanced procedure exhibits a much-improved agreement with the mesoscale model. However, in the case of the bare arch loaded at the quarter span, the continuum model shows a premature shear failure which leads to a significant underestimation of the maximum load and displacement capacity of the arch (Fig. 13c). In the case of bare arch loaded at the mid span (Fig. 13a), the macromodel calibrated by the advanced procedure underestimates the peak-load value. It provides, however, an adequate prediction of the residual strength and pre/post peak response. Finally, a satisfactory comparison can be observed in the cases of bare arch loaded at one eighth span and the confined arch specimen (Fig. 13b and d).
The failure mechanism of the bare arch loaded at quarter span, displayed in Fig. 14a, is characterised by an evident punching effect due to the shear failure of the element of the arch underneath the load, which is not observed in Fig. 7a. On the contrary, the failure mechanism of the confined arch (Fig. 14b) is rather consistent with the mechanism obtained by the mesoscale model (Fig. 7b). The sliding between adjacent rings is represented by shear failure of the solid elements; however, unlikely the continuum model calibrated by means of the simplified procedure (Fig. 11b), the punching effect is not observed.
In conclusion, it can be affirmed that the rigorous calibration procedure allows for a substantial improvement of the continuum model predictions in representing the confined arch specimen. However, the calibrated continuum model still shows evident limits in simulating the bare arch response, as it provides an unrealistic failure mode underestimating both the strength and ductility of the arch.
It is worth pointing out that despite the fact that the adopted calibration strategy is based on the use of the entire arch specimen for the virtual test, its applicability is still computationally efficient as it applied on a simple 2D strip model neglecting the backfill and its interaction with the arch.

Hybrid model simulations
In this section, the proposed hybrid model described in Section 2 is employed to simulate the bare and the confined arch of Fig. 5. The same mesh with quadratic solid elements considered in the previous macroscale continuum simulations has been employed, but circumferential nonlinear interfaces (Section 2) have also been introduced to connect each pair of adjacent solid elements along the thickness of the arch according to the proposed hybrid macroscale representation (Fig. 1b). The model material parameters of the solid elements and the circumferential interfaces are calibrated through the rigorous procedure described in Section 3 and applied before to determine the material properties for the continuum The variation ranges of the unknown parameters are reported in Table 5. The calibration was performed for both weak and strong masonry by using the same procedures and objective functions as for the continuum model. Fig. 15a displays the solutions belonging to the PF (x Ã 1 À x Ã 2 ) for the calibration of the weak masonry model, while the corresponding load-displacement curves are reported in Fig. 15b, c. Comparing these results to the solutions of the optimisation for the continuum model in Section 4.3.2, it can be observed that: -The solutions of the PF are uniformly distributed in the space x 1 -x 2 and are associated with much lower errors.
-The load-displacement curves are less dispersed and much closer to the mesoscale predictions.
These remarks confirm that the further free parameters included in the optimisation algorithm (k Ã ; c Ã ; G M t ) effectively improve the quality of the results. In this case, the absolute minimum of x 2 (4:88 Á 10 À3 ) is reached for a value of x 1 ¼ 4:97 Á 10 À3 while the minimum of x 1 (4:66 Á 10 À5 ) corresponds to x 2 ¼ 7:57 Á 10 À3 . Following the procedure in Section 4.3.2, the two errors are normalised leading to a minimum solution error  Table 5. Contrary to what is observed in the case of the continuum model (Fig. 12d), here the minimum and PF solutions provide the same failure mechanism (Fig. 15d), confirming the ability of the hybrid model to represent the sliding mechanism among adjacent rings.
The PF for the model with strong masonry is shown in Fig. 16a.
The minimum error of x 2 ¼ 0:011 is reached for x 1 ¼ 0:0056, while the minimum error x 1 ¼ 0:0012 is associated with x 2 ¼ 0:337. The minimum solution error corresponds to x Ã min ¼ 0:3207, (A in Fig. 16a) and the matching set of parameters are indicated in Table 5. It can be noticed that the latter solution is characterised by a low value of the dimensionless cohesion (c*=0.69). This circumstance may potentially lead to a response characterised by sliding between the ring which is in disagreement with the mesoscale results. For this reason, another solution (B in Fig. 16a), corresponding to x Ã B ¼ 0:4673, is also considered. This solution has been chosen as the solution with minimum error among those characterised by c Ã > 1.
The response curves of the PF solutions are shown in Fig. 16b and c with the two reference solutions reported with dashed lines. Two families of curves, which tend to minimise the two objectives separately are visible. Although the PF is less regular than the previous case, the selected solutions confirm a good match with the mesoscale curves. Finally, Fig. 16d depicts the failure mechanisms corresponding to the two selected solutions, A and B. It is observed that the two failure mechanisms are relatively consistent with each other. However, the failure mechanism associated with solution A comprises sliding along the circumferential interface, evidencing ring separation not observed in the virtual test. On the contrary, Solution B provides a better approximation of the flexural failure mechanism of the virtual test, thus is corresponds to the best solution to calibrate the macroscale hybrid description for the strong-masonry arch.
Analogously to the procedure for the continuum model in Section 4.3.2, the results of the calibration are validated considering the bare and confined specimens, plus two further models representing the bare arch subjected to a concentrated force at mid span and at one eighth span. The load-displacement curves of the hybrid model are displayed in Fig. 17, where they are compared against the mesoscale curves and the predictions of the continuum macroscale model calibrated by the advanced procedure in Section 4.3.2.
The results of the hybrid macromodel are in a good agreement with those obtained by the mesoscale model confirming a generally improved prediction compared to the results provided by the continuum macromodel. The only exception is represented by the load condition with the force at one eight span, where the continuum model provides a better prediction of the peak loads. However, the curve of the hybrid model, also in this case, is more consistent to the mesoscale response in the pre-and post-peak stages. It appears that the presence of the backfill reduces the differences between the results, with mesoscale and hybrid model predictions almost coincident.
The failure modes of the bare arch and the arch interacting with backfill are displayed in Fig. 18, where the von-Mises equivalent stress distribution in the solid elements and the damage contours on the interface elements are also shown. A good agreement between the failure modes of the hybrid model and those obtained by the mesoscale model, both in terms of stress distribution (Fig. 7) and damage index distribution along the ring-to-ring and arch-tobackfill interfaces (Fig. 8). Importantly, the ring sliding mechanism, which could not be predicted by the continuum model, is well described by the proposed hybrid macroscale representation. The results of the calibration analyses in term of load-displacement curves and failure mechanisms for the strong masonry model are shown in Figs. 19 and in 20, respectively. The curves obtained using the solution B parameters indicate a good agreement with the mesoscale results for all the considered models and loading conditions. Conversely, solution A leads to a significant underestimation of the ultimate strength of the structure in the cases of the bare arch loaded at mid-span (Fig. 19a) and the arch interacting with backfill (Fig. 19d). The failure mechanisms obtained by the macromodel calibrated with the solution B (Fig. 20) result in a good agreement with those obtained by the mesoscale model (Fig. 7). In conclusion, it can be stated that the advanced calibration procedure leads to a realistic set of mechanical parameters to describe the global response of weak and strong masonry arches, under a wide range of boundary and loading conditions. The adopted approach to select the reference solution from the results of the multi-objective optimisation procedure, based on the analysis of the Pareto Front, appears to be sufficiently accurate and robust. However, the circumstance by which the minimum error solution (A) provided less satisfactory results compared to those obtained by another solution belonging to the PF solution (solution B) denotes that further improvements to the selection strategy from the PF and/or the definition of the multi-objective optimisation problem may be needed. This open issue will be investigated by the authors in future studies.

Numerical simulations of long span masonry arches and bridges
In this section, the hybrid model and the advanced calibration procedure are employed to analyse a masonry arch with seven rings and a 16 m span. The numerical results, obtained using ADAPTIC [31], enable investigating the computational efficiency of the proposed macroscale strategy and the improved accuracy due to the use of various nonlinear circumferential interface layers for modelling large masonry arches typical of realistic masonry bridges. Similar to the discussion in Section 4, the performance of the macroscale hybrid description is assessed based on loaddisplacements curves and predicted failure mechanisms which are compared against computationally expensive mesoscale simulations.
The considered masonry arch specimen with stretcher bond is characterised by a 770 mm thickness and a span-to-rise ratio of 4.0. The considered material properties for the masonry constituents are those for the weak masonry material reported in Tables 1 and 2. Fig. 21 shows the geometrical layout of the mesoscale model for the specimen with a representative 1 m width. The vertical loads are applied at the top of the backfill at quarter span of the arch, assuming the same boundary conditions described in Section 4, where the bases of the arch and the backfill are fixed, and the vertical ends of the backfill are restrained against horizon-  The virtual test requires loading the bare arch by the masonry self-weight and two constant symmetric forces of 22 kN applied at one quarter and three quarter span, after which an increasing force is applied at one quarter span until the failure of the arch. Three hybrid macroscale models are adopted, utilising one, two and three circumferential interface layers equally spaced along the thickness of the arch. Table 6 reports the number of elements and DOFs for the three arch macromodels (1-Interface, 2-Interface and 3-Interface), which require 3.3%, 7.7% and 10.3% of the total DOFs of the detailed mesoscale model.
The macroscale mechanical parameters resulting from the calibration procedure are summarised in Table 7. The minimum error solution is considered for each model, while the parameters corresponding to the two ends of the Pareto Front are also reported only for the basic mesh with one circumferential interface layer (1-   Interface). It can be observed that the model parameters for the three solutions of the 1-Interface model are very close highlighting the robustness of the calibration algorithm. The load-displacement capacity curves and the deformed shape at failure predicted by macroscale 1-Interface model are shown in Fig. 22a-b, where these are compared against the mesoscale results. Fig. 22a presents the response curves of the virtual test used for model calibration, whereas Fig. 22b shows the results for the arch confined by backfill (bridge specimen) which have been considered for model validation. In the graphs, the response curves of the continuum model without circumferential interface layers and calibrated by the practical/empirical procedure described in Section 4.3.1 are also reported for comparison.
In the case of the bare arch, the response obtained using the proposed macroscale hybrid model is very close to the mesoscale results confirming the effectiveness of the developed calibration procedure. Conversely, the continuum model leads to a significant overestimation of the arch load-bearing capacity. The failure mechanisms of the macroscale and mesoscale descriptions for the bare arch are shown in Fig. 23, where a good agreement can be observed. Both models predict a mixed failure mechanism characterised by the activation of three flexural plastic hinges and ring sliding in the two portions of the arch close to the skewbacks (Figs. 23a and c). The damage contours at the mesoscale interfaces and those at the circumferential interface layer of the macroscale 1-Interface model are shown in Figs. 23b and d, respectively. Comparing the load-displacement curves of the arch interacting with backfill (Fig. 22b), it is observed that the hybrid model predicts well the response of the system until about 10 mm. For larger vertical displacements, the mesoscale model shows softening beha-     viour, though the response predicted by the macroscale 1-Interface models with selected calibrated material properties exhibits slight hardening. Such discrepancy translates to some differences in the characteristics of the failure mode shown by the deformed shapes and stress contours in Fig. 24. More specifically, the mesoscale model predicts distributed ring sliding leading to local shear failure associated with a drastic reduction of the arch resistance. On the other hand, the macroscale model shows sliding and flexural plastic deformations on the two continuum portions of the FE mesh separated by the single damaged circumferential interface layer which is not associated with sudden degradation of the loadbearing capacity of the arch. Fig. 25 shows the load-displacement mesoscale and macroscale curves for the bare arch and the bridge specimen. It can be seen that each calibrated macroscale model reproduces very well the arch response as determined by the baseline mesoscale model (Fig. 25a). Furthermore, the models with two or three circumferential interface layers (2-Interface and 3-Interface models) lead to improved predictions in terms of ultimate load and post-peak response compared to the model with a single interface layer.
The failure mechanisms depicted in Fig. 26 confirm that the use of multiple interface layers allows a more realistic representation of the distributed shear sliding failure mode, where no notable differences between 2-Interface and 3-Interface models can be observed.
Finally, the efficiency of the proposed hybrid modelling strategy is evaluated against the mesoscale results. In the case of the bare arch, the computing times required by the macroscale 1-Interface, 2-Interface and 3-Interface models are respectively 0.55%, 0.17% and 0.22% of the wall clock time required by the mesoscale model (Fig. 27a). In the case of the arch bridge specimen, the corresponding computing times of the hybrid models are respectively 3.4%, 9.8% and 10.0% of the mesoscale time. These results confirm the much enhanced efficiency guaranteed by the proposed modelling strategy for realistic simulations of large arch bridges, where the use of the mesoscale modelling approach may become computationally prohibitive. In the main, the introduction of multiple interface layers does not lead to a significant increase of the computing time which is, for the analysed cases, always less than 10% of that required by the mesoscale model.

Conclusions
In this study, a hybrid continuum-discrete macro-modelling description for brick-masonry multi-ring arches and arch bridges is proposed. According to this modelling approach, the arch and backfill domain are modelled by 3D continuum solid elements, while layers of 2D zero-thickness nonlinear interfaces arranged along the circumferential direction of the arch simulate potential ring separation and the interaction between the arch and backfill.
Two advanced damage-plasticity constitutive models are employed for the 3D solid and interface elements. An effective and robust multi-level calibration procedure, based on minimisation of stress power error solved by means of Genetic Algorithms is developed to evaluate the model mechanical parameters employing the results from detailed mesoscale models on virtual experiments. These parameters can be easily calibrated from non-destructive or low-destructive in-situ tests on masonry units and mortar joints, which renders the proposed calibration procedure suitable also for practical assessment of historical bridges.
The accuracy and potential of the proposed modelling strategy and calibration procedure is demonstrated by analysing 2D-strip medium and long span masonry arch specimens, also interacting with backfill and characterised by different failure mechanisms. The results of the hybrid model are compared to those obtained by detailed mesoscale and continuum finite element macroscale descriptions. It has been found that the proposed modelling strategy provides accurate predictions of the ultimate strength and displacement capacity of multi-ring masonry arches and the corresponding failure mechanisms. It allows for potential shear sliding between adjacent rings, where the use of only one circumferential interface layer is suitable for medium span arches but at least two layers are required for long span arches and bridges. Additionally, the proposed modelling strategy guarantees superior computational efficiency especially for the analysis of long span arches and bridges, where the use of detailed mesoscale modelling can become infeasible.
Importantly, the numerical results identified some drawbacks associated with the use of conventional isotropic finite element macromodels, which can be summarised as following: -The use of continuum finite element macromodels, without a rigorous calibration of the mechanical parameters, can lead to an inaccurate and non-objective prediction of the arch response. -Continuum finite element macromodels, even when calibrated by means of rigorous procedures, can fail in simulating the ultimate arch behaviour when it is driven by sliding between adjacent rings.
Both these limitations may significantly affect the results of safety assessments of masonry arch bridges, leading to a crudely approximated or completely misleading prediction of the effective safety level of the bridge and its mode of failure. In this regard, the proposed hybrid modelling strategy offers the possibility to significantly improve the accuracy of the numerical predictions without requiring a significant increase in the computational effort associated with the nonlinear analysis.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.