Numerical Simulation of a Masonry Arch Bridge with Initial Defects Based on Cohesive Elements

: Most of the existing masonry bridges have been in service for a signi ﬁ cant duration, and as a result of construction limitations, these structures often exhibit intricate geometric defects. Furthermore, under prolonged loading conditions, the rheological behavior of rock can induce deformation in masonry bridges, leading to a continuously evolving stress state. Employing an idealized model for safety assessment frequently results in an overestimation of their load-bearing capacity. To accurately evaluate the load-bearing performance and remaining service life of masonry bridges, as well as to prevent safety incidents, this study employs a parametric approach to establish a two-phase numerical model of masonry bridges. In this model, cohesive elements are introduced to simulate the bonding relationship, while the distribution pa tt ern of geometric initial defects is determined based on the theory of conditional random ﬁ elds. Additionally, the rheological behavior of rock is incorporated through a custom-wri tt en Abaqus user subroutine. Building upon this foundation, the probability distribution of the load-bearing capacity of masonry bridges is reconstructed using the maximum entropy method with fractional moment constraints. The resulting outcomes are compared and validated against those obtained using the decomposition conditional correlation matrix. Finally, the e ﬀ ectiveness and applicability of the proposed method are demonstrated through numerical simulations and ﬁ eld measurements conducted on an actual bridge. The ﬁ ndings reveal that the method introduced in this paper adequately accounts for the stochastic nature of geometric initial defects, objectively re ﬂ ects the operational performance of masonry bridges, and e ﬀ ectively simulates the complete failure process of such structures. Consequently, this method provides a solid basis for the safety assessment of masonry bridges.


Introduction
Stone masonry arch bridges have a global presence, as evidenced by numerous studies [1][2][3].In the southern and southwestern mountainous regions of China, where stone resources are abundant but transportation is challenging, the proportion of stone masonry arch bridges exceeds 90%.Most of these bridges were constructed based on the practical knowledge and expertise of local craftsmen.However, as these bridges age, various issues such as deformation, cracking, and damage have become apparent [4].The ability of these bridges to withstand current transportation demands is a significant concern.To ensure the preservation of stone masonry arch bridges, it is imperative to propose a method that can effectively account for their initial defects, as well as the impacts of existing damage and deformation.Such an approach would enable an accurate assessment of their loadbearing performance and provide a solid theoretical foundation for their reinforcement and maintenance.
The finite element method (FEM) is widely utilized for simulating the load-bearing performance of stone masonry structures [5][6][7][8].Currently, there are several mainstream modeling methods available for studying masonry materials.For instance, macroscale modeling involves the introduction of equivalent nonlinear constitutive models to describe the material's failure process.Karaton.M et al. [9] obtained material properties for homogenized modeling through experiments and evaluated the nonlinear seismic response of the Malabadi Bridge, built in the 12th century, using homogenized modeling methods.Aydin et al. [5] evaluated the structural performance of the Salpeder Bridge using homogenized modeling methods.However, these methods do not fully consider the heterogeneity of masonry materials, leading to discrepancies between the analysis results and the actual behavior.
On the other hand, microscale modeling [10][11][12][13] treats the material as a discrete system composed of particles and investigates the mechanical properties, fracture behavior, particle arrangement, etc., by simulating the interactions and movements between particles.For example, Micaela Mercuri et al. [12] studied the fracture behavior and size effect of masonry dome structures through microscale modeling.However, microscale modeling requires significant computational resources, which may pose challenges when addressing engineering problems.In summary, while macroscale modeling with equivalent nonlinear constitutive models has been used to assess the behavior of masonry structures, it overlooks the material's heterogeneity.Microscale modeling, on the other hand, considers the discrete nature of masonry materials but demands substantial computational resources.
To address this limitation, Zhao [14] developed a two-phase numerical model for stone arch bridges that incorporates cohesive elements.Experimental verification has demonstrated the effectiveness of this two-phase numerical model in simulating the mesoscale failure process of the structure.The mesoscale simulation also has certain limitations, such as the inability to capture crack propagation.It captures defects under extreme conditions.This paper adopts the mesoscale modeling approach due to its ability to integrate macroscopic and microscopic features and provide higher efficiency required for engineering problems.Additionally, due to construction-related factors, geometric initial defects are inevitably present in stone masonry arch bridges [15][16][17][18][19][20].
The presence of geometric initial defects in stone masonry arch bridges can significantly influence their mechanical performance [21,22].These defects can lead to various changes in the behavior of the bridge.First, geometric asymmetry and nonuniform deformation caused by defects can result in nonlinear responses in the arch bridge.This can lead to stress concentration and uneven stress distribution, ultimately affecting the loadbearing capacity and stability of the bridge [23].Second, defects can exacerbate the deformation and displacement of the bridge, further reducing its structural stability and potentially leading to structural damage or collapse [24].Currently, scholars primarily consider initial defects using low-order modal methods.However, due to the complexity and highly random nature of actual defects, low-order modal methods are unable to accurately capture the load-bearing capacity and deformation development trends of stone arch bridges [25].It is also important to note that rocks exhibit viscoelastic behavior [26,27].The scouring of water flow can also have an impact on the bridge [28].During the extended service life of stone masonry arch bridges, intricate stress fields develop within the structure, and the mechanical properties of rock materials may undergo deformation.Hence, it is essential to analyze the load-bearing capacity and deformation evolution of these bridges over time and assess their future load-bearing and deformation conditions.This study aims to address the material heterogeneity of stone masonry arch bridges by employing a mesoscale modeling approach.In this approach, a two-phase model is developed, distinguishing between the masonry blocks and mortar joints.To describe the mechanical behavior of the masonry blocks, solid elements are utilized.These elements capture the response of individual blocks, considering their material properties and geometric characteristics.According to Heyman's work [29], it is considered that the compressive strength is infinitely large.On the other hand, cohesive elements are introduced between adjacent blocks to represent the failure mode of the mortar joints.These cohesive elements simulate the bond between the blocks and account for the behavior of the joints under different loading conditions.By adopting this two-phase modeling approach, the study aims to capture the material heterogeneity of stone masonry arch bridges and simulate the mesoscale failure process of the structure more accurately.This can provide insights into the load-bearing capacity, deformation behavior, and failure mechanisms of the bridges, ultimately leading to improved understanding and analysis of their structural performance.In addition, the adopted model can accurately reflect the size effect, leading to more precise and reliable results.By considering the influence of size [11,12] on structural behavior, the model allows for a more comprehensive analysis, ensuring that the obtained results are both accurate and safe.This consideration is crucial for engineering applications, as it enables a better understanding of the structural response under different size conditions and aids in making informed design decisions.Furthermore, to accurately consider various potential initial defects in stone arch bridges, a defect simulation method based on the theory of conditional random fields is introduced.Subsequently, a parameterized modeling approach is employed to conveniently establish finite element models with different initial defects for analysis.The load-bearing capacity of the structure is simulated based on the risk method [30].Finally, the method proposed in this paper is applied to analyze an actual bridge, accurately retracing its deformation process, and based on this, analyzing its load-bearing capacity.The method presented in this paper can provide a certain reference for the analysis of the load-bearing capacity and deformation of stone arch bridges.

The Establishment of a Two-Phase Numerical Model
This study utilizes Abaqus to establish the model.The specific modeling methodology and the selection of relevant parameters are outlined as follows.

Parameterized Modeling
Parameterized modeling refers to the process of creating digital models through the utilization of preprogrammed algorithms [31,32].This modeling approach enables the automatic generation of models by employing internal logical parameters, thereby circumventing the need for manual construction.By leveraging this methodology, the accuracy and efficiency of modeling can be significantly enhanced.In the present study, a parameterized modeling approach was employed to represent the material heterogeneity inherent in stone masonry arch bridges.To this end, the Grasshopper software, which operates on the Rhino platform, was adopted.Grasshopper, functioning as a visual programming software, encapsulates functions into discrete "components" and employs "wires" to establish logical connections among these components.This software offers considerable flexibility, allowing for the integration of custom scripts to address intricate problems.The modeling workflow implemented in this study encompassed the following stages: (1).Creation of the initial geometric model: A preliminary two-dimensional geometric model was constructed based on the provided parameters, such as the span and rise of the bridge.(2).Calculation of initial geometric defects: Custom Python scripts were developed to calculate the initial geometric defects, which were subsequently incorporated into the initial geometric model.(3).Model export: A three-dimensional model was derived from the two-dimensional representation, and a bespoke ABAQUS scripting program was employed to facilitate the seamless integration of Grasshopper with ABAQUS for conducting batch computations.
By employing this parameterized modeling approach, the primary objective of this study was to augment the accuracy and efficiency of modeling stone masonry arch bridges.This, in turn, enables a more comprehensive analysis of their load-bearing capacity, deformation behavior, and future performance.

Rock Constitutive Models
The rock undergoes continuous deformation under long-term loading due to the viscoelastic effect [33][34][35], which is particularly significant for long-serving stone arch bridges and needs to be considered.This paper employs the Burgers model [36] to simulate the rheological behavior of rocks.The Burgers model is a viscoelastic model that consists of a series connection of the Maxwell model and the Kelvin model, as shown in Figure 1: where  and  are the shear modulus and viscosity coefficient of the Kelvin body, and  ,  ,  are the bulk modulus, shear modulus, and viscosity coefficient of the Maxwell body, and   is the spherical stress tensor.The present study develops a corresponding Abaqus subroutine based on the aforementioned principles to capture the rheological effects of rocks.

Cohesive Elements Constitutive Models
In this study, cohesive elements are employed to simulate the filling material between the masonry units and capture the bonding and sliding behavior between them.Different cohesive constitutive models can significantly affect the numerical simulation results.Camanho [37] proposed a bilinear constitutive model.To more accurately describe the mechanical properties of cohesive elements, this paper adopts a linear-exponential constitutive model, as shown in Figure 2, to describe the stress-displacement relationship.σ n and shear stresses σ s and  t satisfy Equation ( 2), the cohesive element experiences damage.At this point, the normal displacement and shear displacements of the element are denoted as  t and δ s, t 0 , respectively.
After the initiation of damage, the stress-displacement relationship enters the damage evolution stage.The damage level is represented by the introduction of a damage index D, which is defined as follows: The equation represents the stress linearly corresponding to the current displacement in the undamaged state.After the element enters the damage evolution stage, its failure is primarily determined based on energy criteria or displacement criteria.In this study, the B-K energy criterion is chosen as the criterion for failure assessment, and its specific form is as follows: Under this criterion, the damage index D associated with the linear-exponential constitutive model is defined as follows, where G n , G s , and G t represent the work performed by the normal and tangential tractions, and G n C and G s C represent the critical fracture energies in the normal and tangential directions, respectively.The parameter  represents the total energy consumed during unit failure.

Model Validation
The validity and accuracy of the proposed cohesive element in this paper are verified by comparing the theoretical solution with the finite element solution of the ultimate buckling load for hingeless circular arches.Initially, the theoretical solution formula for the ultimate buckling load of hingeless arches is derived.Pi et al. [38] employed the virtual work principle method to establish the equilibrium equation and obtained the bending moment of the hingeless circular arch under hydrostatic pressure as follows: where 2Θ is the arch opening angle;  is the angular coordinate; R is the curvature radius of the circular arch.µ e is a dimensionless axial force parameter, defined as: P is a dimensionless parameter that measures the nominal axial pressure qR and the actual axial pressure N, and is expressed as: The nonlinear equilibrium equation between the dimensionless axial force parameter and the nominal axial pressure is: where: where  is the modified slenderness ratio of the arch ring, defined as When the arch ring joints experience cracking due to tension, an increase in load can lead to rapid crack propagation and eventual structural failure.Therefore, it is assumed that the structure reaches its maximum bearing capacity when cracks occur.At this point, the tensile stress induced by bending and the compressive stress caused by compression at the cracking location of the structure should satisfy the following equation: where h is the thickness of the arch ring.By Equations ( 6) and ( 7), we have: By combining the above two equations, we obtain where: The values of P and  at joint failure can be determined, allowing for the subsequent calculation of the external hydrostatic pressure.In this paper, models with various slenderness ratios and different rise-span ratios, all with a span of 10 m, are selected for analysis and comparison.
The key parameters of these models are presented in Table 1.The parameter  , is the ratio of peak stress to the corresponding displacement as shown in Figure 2.
The size of computational elements significantly impacts the accuracy and reliability of numerical simulations.In this study, we investigated the influence of varying mesh sizes on the numerical simulation results for a non-articulated circular arch characterized by a rise-span ratio of 0.2 and a slenderness ratio of 100.The results, as illustrated in Figure 3, provide insights into the correlation between mesh size and the load-bearing capacity of the arch.It is evident from the graph that a reduction in mesh size leads to a decrease in the arch's load-bearing capacity.This decline can be attributed to the premature attainment of the predetermined stress state by the edge elements due to the decreased mesh size, resulting in their removal and subsequently compromising the overall load-bearing capacity of the arch.
However, it is noteworthy that as the number of elements exceeds 8, the load-bearing capacity tends to stabilize.Based on this observation, the present study selects a mesh quantity of 10 to ensure a reliable and accurate representation of the arch's load-bearing behavior.
Figure 4 illustrates the theoretical and numerical solutions for the bearing capacity of hingeless circular arches under different rise-span ratios and slenderness ratios.These curves represent the theoretical values and finite element numerical solutions of the ultimate buckling load.According to Figure 4, it can be observed that the load-carrying capacity of the hingeless circular arch decreases with an increase in the slenderness ratio, while it increases with an increase in the rise-to-span ratio.These findings are consistent with existing research results [39,40].Across different rise-to-span ratios and slenderness ratios, the theoretical and numerical solutions for the load-carrying capacity of the hingeless circular arch structure are closely aligned, with an average deviation of 8.5%.This high level of agreement demonstrates the effectiveness of the modeling method used in this study.

Initial Defect Simulation Based on Conditional Random Field (CRF) Theory
In the parameter space ΩR n , for any position vector ω∈Ω in the space, X(ω) is a random variable [41][42][43][44][45].In the actual structure, the defect distribution at different positions is correlated, and the correlation coefficient function ρ(ω 1 ,ω 2 ) is used to express it.The boundary conditions of the arch structure have a great influence on the overall mechanical performance of the structure.For hingeless arch structures, when the geometric defects at other parts of the structure remain random, the randomness at the structure supports disappears.A conditional random field is introduced to reflect this phenomenon.For the n-dimensional random variable in the parameter space, the covariance matrix of the variable  = [ ,  , … ,  ] is: In the equation: µ i and σ i 2 are the mathematical expectation and variance of ω i , respectively; C(ω i ,ω j ) represents the covariance of ω i and ω .If there are k dimensions known in the n-dimensional random variable u, then u can be divided into two parts of k dimensions and n-k dimensions, respectively denoted as  and  , and their expectations µ and covariances C are also expressed accordingly as: When () is a normal distribution random field and  is known, the normal conditional distribution probability density function of  is: where: The boundary conditions of the actual existing arch structure are known, assuming that the boundary conditions are 0 and 1.In this case, the conditional random variable  () satisfies: Given values ū1 and ū2, let ω and ω ' represent any two arbitrary random points within the arch structure, excluding the boundaries.The corresponding mathematical expectation and covariance matrix in Equation ( 16) can be expressed as follows: ,  ′ =   ,  ′ − (0, ) 0,  ′ + (, 1)  ′ , 1 1 −  (0,1) + (0,1)  0,  ′ (, 1) + (0, )  ′ , 1 1 −  (0,1) According to Equations ( 18) and ( 20), the conditional correlation function can be obtained as follows: ,  ′ =   ,  ′ − (0, ) 0,  ′ +(, 1)  ′ , 1 The degree of correlation between node deviations in an arch structure depends on the distance between the nodes.If two nodes are close to each other, their deviations are highly correlated; whereas, if they are farther apart, the degree of correlation is smaller.According to the literature references, the correlation between two nodes in the model follows an exponential correlation coefficient function where the correlation decays exponentially with increasing distance.The correlation coefficient function can be expressed as: With  being the correlation length parameter, substituting the correlation function into Equation ( 22) yields the corresponding conditional correlation matrix, denoted as C ωω ' .Decomposing the obtained conditional correlation matrix allows for the characterization of defect distribution when the geometric initial defects are treated as random variables.
For the case of a hinged circular arch, with a span of  = 10 m, a slenderness ratio of  = 0.3 an aspect ratio of  = 100 , and divided into 50 equal segments along the arc length, we assume that the deviation at each node within a segment follows a normal distribution with a mean of 0 and a standard deviation of  = 0.01 m.The correlation length in the correlation function is  = 1 .We will now decompose the conditional correlation matrix and present the defect distribution for the first three orders, as illustrated in Figure 5.

Initial Defect Verification Based on Conditional Random Field Theory
The introduction of a probability distribution model for the ultimate bearing capacity of a structure, considering geometric initial defects based on the conditional random field (CRF) theory, can effectively reflect the actual geometric initial defects in the structure.To validate the accuracy of the method proposed in this paper, which applies initial defects based on CRF theory, a probability distribution model for the ultimate bearing capacity of the structure incorporating geometric initial defects is constructed using the Latin hypercube sampling (LHS) [46] and fractional moment maximum entropy (FMME) [47,48] methods.
Assuming there is a unit hypercube in a D-dimensional space, denoted as  , it is decomposed into M orthogonal subspaces of dimension  where ( = 1,2, ⋯ ).These subspaces satisfy the following conditions:  =   ⊥  ,  ≠  (24) In each of the lower-dimensional subspaces Θ , ̄ samples points are randomly selected.They are denoted as: The obtained sample points are randomly paired, resulting in sample points in the original space as: where Ξ represents the random pairing function.
The obtained covariance matrix is subjected to Cholesky decomposition.The resulting lower triangular matrix is then multiplied by the sample points obtained from Equation (26) to obtain the samples of the geometric defect random field.
To reconstruct the probability distribution of the ultimate bearing capacity, the thorder fractional moment matrix [18] of the random variable  is defined as: = () represents the transformation function applied to the input random variable  to obtain the output random variable .The function  () represents the probability density function of the output variable.Combining with Equation ( 26), we have: Based on the principle of maximum entropy [49], combined with Equation ( 27 According to the key condition for obtaining the optimal solution from the above equation, the estimated probability density function  () is obtained as: We can obtain: After obtaining the estimated probability density function  , the Kullback-Leibler (K-L) divergence based on the true probability density  and the estimated probability density  is expressed as: Due to [] being the entropy of the true probability density  , it is independent of  and .Therefore, minimizing the Kullback-Leibler (K-L) divergence implies minimizing the following expression: Based on the above equation, the mathematical representation of the maximum entropy problem based on the fractional moment matrix is: In the specific calculation process, the original high-dimensional space is divided into LPSS (low-dimensional projection subspace) forms.In each subspace, 1000 sample points are extracted and randomly combined to generate initial geometric defects.The parametric modeling method is used to establish a geometric model, and the obtained initial geometric defects are applied.Nonlinear analysis is performed to obtain the corresponding extreme buckling load.
The Lagrange multipliers  and the order of the fractional moment matrix  obtained by solving the problem using the downhill simplex method based on Equation ( 36) with m = 3 constraints are presented in Table 2.The probability density is reconstructed by substituting them into Equation ( 31) as shown in Figure 6.From Figure 6 and Table 3, it is evident that the probability density reconstructed using the Latin hypercube sampling method based on the fractional moment matrix predominantly falls between the stable bearing capacity of structures with random secondorder distribution defects and those exhibiting perfect structures.Furthermore, the reliability index calculated from the probability distribution obtained via the LPSS method, corresponding to β = 3.7, closely corresponds to the bearing capacity associated with the second-order antisymmetric geometric initial defects of the random distribution.This substantiates that the imposition of random second-order distribution defects can, with the use of fewer computational resources, reconstruct the most adverse defect distribution for structures in most cases, thereby providing a certain guaranteed bearing capacity.

Experiments
The lower-order defects obtained from the conditional random field method are the most probable forms of defects, and lower-order defects have a greater impact on the structure.The initial defect distribution in the actual structure theoretically approximates the lower-order defects obtained from the conditional random field method.To validate the applicability of the method, we take an actual stone arch bridge as an example.The Yucheng Bridge is located on the eastern side of the Dongjiang River, northwest of Futou Village, Gulai Town, Shengzhou City, Zhejiang Province.It is a single-span elliptical stone arch bridge constructed in the 16th year of the Qing Dynasty (1836).It is designated as a provincial-level cultural heritage site in Zhejiang Province.The appearance of the bridge is depicted in Figure 7. Based on Table 4, it can be observed that the best-fitting result is achieved using an elliptical arc.Therefore, in this study, the arch ring is approximated using an elliptical arc.Through analysis, the equation of the ellipse can be obtained as follows: Based on the aforementioned equation, a geometric model was constructed with the geometric parameters as shown in Table 5. Cohesive elements were employed to simulate the contact of masonry, while the Mohr-Coulomb model was utilized to simulate the behavior of the backfill [6,50,51].Using the aforementioned theory of conditional random fields, the stone masonry nodes are regarded as random points, and a parameterized modeling method is employed to apply second-order antisymmetric geometric initial defects, as illustrated in Figure 8:  Based on experimental data illustrated in Figure 10 and the relevant literature, the constitutive parameters of rock considering the rheological effect are presented in Table 6:   The parameters for cohesive elements are set as shown in Table 7.According to the study by Mercuri et al. [11,12], the self-weight effect of arch structures has a significant impact on size effects.For large-scale structures, as the size increases, the structure experiences a decrease in equivalent load-carrying capacity due to the non-negligible self-weight effect.The model proposed in this paper also reflects this effect.To verify the size effect, the model shown in Figure 11a was subjected to the load shown in Figure 11b, and the load-carrying capacity is illustrated in Figure 12.According to Figure 12, it can be observed that as the structure size is scaled up proportionally, the equivalent load-carrying capacity decreases.Furthermore, the relationship between the two variables exhibits a strong linear correlation.This demonstrates that the model used in this study can reflect the size effect.Therefore, when simulating the construction of the Yucheng Bridge, it is important to consider the effect of gravity.
The material parameters according to experiments and a range of typical values are shown in Table 8: It should be noted that in the normal constitutive relationship, to prevent intrusion between elements, the normal stiffness k n 0 is defined as a maximum value.In the tangential constitutive relationship, the normal stress perpendicular to the mortar joint affects its shear strength σ s,t 0 .Since specifying it as a constant cannot fully reflect the real situation, it is set as  , =  −    (37) In the formula,  is the cohesive force of the masonry joint;   Is the shear friction coefficient of the masonry joint.
The UTRACLOAD subroutine is used to apply a moving load with a magnitude of 0.003 MPa.The load spans 1 m and acts vertically downward on the top of the backfill area, repeating uniformly at a constant speed along the longitudinal direction.The CREEP analysis is set to a time step of 180 years.The numerical simulation results are depicted in Figures 13 and 14: As evident from Figure 13, the structure exhibited significant deformations corresponding to its own dimensions, displaying pronounced asymmetry.Extracting the finite element simulation results as depicted in Figure 13, the deformed arch axis of the Yucheng Bridge was obtained, as shown in Figure 14, and compared with the actual measured arch axis.To validate the accuracy of the finite element numerical simulation, a comparative analysis was conducted between the simulated results and the measured deformation of the arch axis.The Pearson correlation coefficient was used as the evaluation metric to measure the linear correlation between the two sets of data.The results showed that the Pearson correlation coefficient between the finite element numerical simulation results and the actual arch axis deformation reached 0.994.This indicates a very strong positive correlation between the two, demonstrating that the proposed method can effectively reflect the deformation characteristics of the arch structure.
During the on-site inspection, it was found that the area marked by the red box in Figure 15 experienced block detachment.As shown in Figure 16, the red box area can be divided into three subregions: I, II, and III.In region I, the blocks exhibit detachment at the lower part and compression at the upper part, while in regions II and III, the blocks detach at the upper part and are compressed at the lower part.Figure 15 shows the removal of cohesive elements in the finite element model.When the damage index D of a cohesive element reaches its maximum value, the cohesive element is deleted to simulate the cracking between the masonry units.It can be observed that the cracking pattern and distribution of the cohesive elements are consistent with the actual cracking situation of the arch ring.This indicates that the proposed method can effectively simulate the cracking and failure of the original structure.To predict the load-bearing performance of the Yu Cheng Bridge arch structure, this study further analyzed its capacity based on the results of long-term deformation analysis and compared it with the intact arch without considering initial defects and the intact arch considering initial defects.The load was applied with a width of 200 mm at one-quarter span of the bridge span.
Figure 17 presents the results of the bearing capacity analysis for different models.As depicted in the figure, the ultimate bearing capacity of the Yu Cheng Bridge is 9553 KN for the structure without any defects.However, after the introduction of geometric initial defects, the ultimate bearing capacity of the bridge decreases to 7566 KN, indicating a reduction of 20.8% compared with the defect-free structure.This demonstrates that geometric initial defects significantly diminish the structural load-carrying capacity, highlighting the potential danger of neglecting such defects.
Under the influence of long-term loads, both the material properties and geometric shape of the structure undergo changes.Considering the effect of long-term loads, the ultimate bearing capacity of the Yu Cheng Bridge decreases to 5440 KN.This represents an 18.1% reduction in load-carrying capacity compared with the arch ring structure considering only initial defects.Furthermore, when compared with the defect-free structure, the load-carrying capacity diminishes by 64.9%.These findings emphasize that long-term load effects also significantly diminish the structural load-carrying capacity.This reduction in load-carrying capacity is even more pronounced for ancient cultural heritage structures.The failure modes of the main arch ring considering the influence of long-term deformation are illustrated in Figure 18.With increasing load, cracks first appear at the bottom of the arch ring at the location of the applied load (one-quarter span).As the load further increases, cracks emerge at a three-quarters span, primarily in the upper portion of the arch ring.Eventually, cracks appear almost simultaneously at one-eighth and seven-eighths spans, causing the structure to transform into a mechanism and reaching its maximum load-bearing capacity.Analysis reveals that the load-bearing capacity of the Yu Cheng Bridge, considering the influence of long-term deformation, can reach 6205 kN.The current load-bearing capacity meets the requirements of a pedestrian bridge.In order to analyze the variation of the load-bearing capacity of the Yu Cheng Bridge with service life, this study considered different creep time steps and two types of pedestrian loads: uniformly distributed load over the full span and half-span.Based on previous research [52] and on-site measurements, it has been determined that the pedestrian load is a cyclic load with a magnitude of 0.03 MPa.Its influence extends 1 m along the length direction and spans the full width along the width direction.
The load-bearing capacity changes after different service times were analyzed and are shown in Figure 19.Performing linear regression on the data yields a good fit, and the results are as follows: In the given context,  ℎ represents the bearing capacity of the structure under a halfspan uniformly distributed load, while  represents the bearing capacity of the structure under a full-span uniformly distributed load.
Based on the above two equations, under the assumption of maintaining the current conditions and considering only pedestrian loads, the theoretical remaining service life of the Yu Cheng Bridge is estimated to be approximately 150 years.

Conclusions
This study establishes a two-phase model of masonry arch bridges based on the cohesive element using a parameterized modeling approach.The introduction of geometric initial defects is incorporated using the conditional random field theory.The rationality of the application of geometric initial defects is verified from the perspective of load-bearing capacity.Field measurements and numerical simulations are conducted to compare and validate the overall methodology, demonstrating its applicability and effectiveness.
1.The use of parameterized modeling methods greatly enhances the accuracy and efficiency of establishing geometric models.When combined with finite element software, it enables the rapid creation of numerical models for different geometric parameters.This provides a solid foundation for reconstructing the probability distribution of ultimate bearing capacity through a large number of samples.2. The two-phase numerical model based on cohesive elements allows for an intuitive and accurate simulation of the noncontinuous and nonlinear cracking process at the masonry contact surface of brick arch bridges.This provides a basis for evaluating the structural load-bearing capacity and failure modes.3. Using the conditional random field theory significantly reduces the structural load-bearing capacity and leads to noticeable asymmetrical deformation.Through comparison with actual structures, the effectiveness of the proposed method is verified.

Figure 1 .
Figure 1.Burgers model.The creep equation of the Burgers model under a three-dimensional stress state is:

Figure 2 .
Figure 2. Self-defined bilinear constitutive relation curve: (a) Normal constitutive model; (b) Tangential constitutive model.In Figure 2, σ n 0 and σ s, t 0 represent the maximum normal stress and shear stress, respectively.This study adopts the Quads damage initiation criterion, where the stress linearly increases with displacement before the initiation of damage.When the normal stress

Figure 3 .
Figure 3.The relationship between mesh size and load-bearing capacity can be represented by a curve, where n denotes the number of elements in the thickness direction of the arch ring.

Figure 4 .
Figure 4. Theoretical and numerical solutions of the bearing capacity of hingeless circular arches under different aspect ratios and slenderness ratios: (a) bridge aspect ratio f/l = 0.2; (b) bridge aspect ratio f/l = 0.3; (c) bridge aspect ratio f/l = 0.4; (d) bridge aspect ratio f/l = 0.5.

Figure 5 .
Figure 5. Geometric imperfections of the first three orders with random distribution in hingeless circular arches: (a) First-order defect distribution; (b) Second-order defect distribution; (c) Thirdorder defect distribution.

Figure 6 .
Figure 6.Statistical diagram of structural stability limit load under different initial imperfections.

Figure 8 .
Figure 8. Geometric initial defect of the arch ring.The three-dimensional finite element model of the Yucheng Bridge is depicted in Figure 9:

Figure 11 .
Figure 11.The validation of the size effect: (a) The structures of different sizes scaled up proportionally; (b) The load applied at the midspan (P represents the applied concentrated load).

Figure 12 .
Figure 12.The relationship between structural size and load-carrying capacity.

Figure 14 .
Figure 14.Comparison of finite element numerical simulation and actual arch axis.

Figure 15 .
Figure 15.Cracking diagram of arch bridge(The area inside red ovals is where the cracks appear): (a) Cracking distribution in the voussoir masonry; (b) Detachment of the lower part of the masonry; (c) Detachment of the upper part of the masonry.

Figure 16 .
Figure 16.Cohesive elements deletion area (The red dashed oval area is the zone where cohesive elements are deleted).

Figure 18 .
Figure 18.Structural failure characteristics(The location where cracks appear is inside the red dashed oval).

Figure 19 .
Figure 19.Relationship curve between bearing capacity and service time of the Yucheng Bridge: (a) uniformly distributed load on the half-span; (b) uniformly distributed load on the full span.

Table 1 .
Main parameters of extreme buckling load comparison model.

Table 2 .
Maximum entropy optimization parameter results.

Table 3 .
Statistical table of structural stability limit load under different initial imperfections.

Table 4 .
Correlation coefficients corresponding to different geometric curves.

Table 7 .
Mechanical properties for nonlinear interfaces.
The data presented in this study are available in the article.The authors declare no conflict of interest.