ArchLab: a MATLAB tool for the Thrust Line Analysis of masonry arches**

According to Heyman’s safe theorem of the limit analysis of masonry structures, the safety of masonry arches can be verified by finding at least one line of thrust entirely laying within the masonry and in equilibrium with external loads. If such a solution does exist, two extreme configurations of the thrust line can be determined, respectively referred to as solutions of minimum and maximum thrust. In this paper it is presented a numerical procedure for determining both these solutions with reference to masonry arches of general shape, subjected to both vertical and horizontal loads. The algorithm takes advantage of a simplification of the equations underlying the Thrust Network Analysis. Actually, for the case of planar lines of thrust, the horizontal components of the reference thrusts can be computed in closed form at each iteration and for any arbitrary loading condition. The heights of the points of the thrust line are then computed by solving a constrained linear optimization problem bymeans of the Dual-Simplex algorithm. The MATLAB implementation of presented algorithm is described in detail and made freely available to interested users (https://bit.ly/3krlVxH). Two numerical examples regarding a pointed and a lowered circular arch are presented in order to show the performance of the method.


Introduction
The analysis of masonry arches is a classical problem of structural mechanics. It was from the intuitions by Hooke, hidden in a famous anagram, and by Gregory, explicitly stated, that started a series of researches employing the catenary analogy to model the equilibrium of arches by means of compressive internal actions that are funicular of the applied loads [11]. Alternative approaches, initiated by La Hire and developed by Coulomb, analyze the equilibrium of arches by studying collapse mechanisms of a series of voussoirs capable to transfer a limited set of actions [28]. A formal explanation of these approaches was given in two renowned papers by Heyman [25] and Kooharian [29], where the safe and unsafe theorems of the limit analysis of masonry structures are presented.
Nowadays several computational methods are available for the analysis of masonry structures. Some of them are based on the Finite Element Method (FEM) and are capable to take into account sophisticated material models [3,32,33,47]. However, an appropriate application of FEM based analyses requires an accurate knowledge of the value and spatial distribution of mechanical properties of materials and support settlements [26,27], which requires detailed survey techniques, unmotivated for ordinary structures.
A interesting alternative is represented by the methods that employ the No-Tension (NT) material model [7,17], which is characterized by a non-smooth behavior and requires suitable techniques to be successfully employed for the analysis of real structures [8,18]. Additionally, the NT model is embodied by the hypotheses of previously mentioned limit theorems. Some of their computational implementations include the Thrust Network Analysis (TNA) [11,12,34,36,44] and the Thrust Surface Analysis (TSA) [20], both based on the safe theorem, or the Rigid Block Analysis [13,24,31] and the Discrete Element Method (DEM) [30,49], which employ the unsafe theorem. Worth mentioning are recent proposals in which a fracture mechanics-based analytical method with elastic-softening of masonry is applied to analyse the structural behaviour of arch bridges and show how the arch thrust line is affected by crack formation [1,2].
Although these methods are all capable to analyze structures characterized by complicated geometries [35,46] and by unusual construction techniques [19,37], specific tools for the analysis of masonry arches are still of scientific interest [6,38,45]. Methods based on the kinematic approach are used to determine the collapse of arches undergoing spreading of supports [16,22]. Methods based on Heyman principles are employed to compute the minimum thickness [4,40] and the later load bearing capacity [5,14] of arches.
Interestingly, the thrust line analysis can also be used to determine the optimal shape of masonry arches subjected to vertical and horizontal loads [39,41]. In particular, Michiels and Adriaenssens [39] compute a unique thrust line by defining its span and height. The determination of this thrust line is obtained by iteratively adjusting the horizontal and vertical forces applied to the nodes of the thrust line. This thrust line is then mirrored with respect to the vertical axis so that the pair of mirrored thrust lines are used to define the profile of the arch, while an iterative procedure minimizes its total volume. The approach proposed by Nikolic [41], instead, employs an analytical modeling of a catenary arch of constant and finite thickness, for which the horizontal thrust is determined. Then a thrust line is analytically determined for this arch. Nikolic shows that this thrust line is not coincident with the catenary curve that defines the axis of the arch because of the altered position of the center of gravity of voissors due to their finite thickness and curvature.
Worth mentioning are the two MATLAB tools ArchNURBS [15] and FRS_Method [23]. The approach implemented in ArchNURBS is based on a nonuniform rational B-splines (NURBS) representation of arch geometry. A preliminary isogeometric finite-element elastic analysis of the arch is performed in order to determine the structural response under service loads, provided that the corresponding thrust line fulfills assumed geometric constraints. Successively a limit analysis based on the safe theorem by Heyman is carried on by considering equilibrium and yielding conditions of blocks interfaces. These equations and constraints are solved as a linear optimization problem that maximizes the load multiplier. The method employed in FRS_Method, insted, computes the line of thrust by solving an optimization problem that looks for the funicular polygon closest to the geometrical axis of the arch [48]. This means that a unique line of thrust is determined, which is neither one of the two limit configurations of minimum or maximum thrust, neither it obeys to the principle of least action. Additionally, the authors employed a new definition of geometric safety factor of the arch, which is alternative to the one given by Heyman.
The objective of this paper is the Thrust Line Analysis (TLA) of masonry arches, a specialization of the TNA to the two-dimensional case. As the TNA represents a dis-crete implementation of the no-tension membrane model [21] the TLA represents a discrete version of the funicular curve. The specialization to a planar line of thrust simplifies enormously the solution of the horizontal equilibrium equations that are used to compute the horizontal components of thrust. Actually, in the proposed version of the method, reference values of horizontal thrust are computed in closed form at each iteration.
The proposed version of the TLA has been implemented in a MATLAB code, freely downloadable from https://bit.ly/ 3krlVxH. The method is implemented within the function ArchLab that takes as input a table containing the geometric and loading data of the arch. Also some ready-to-run example files are provided together with a plotting function, useful to graphically visualize the solutions.
After describing the general method, this paper illustrates in detail how each portion of the algorithm has been implemented in the MATLAB function ArchLab. Results regarding the analysis of a pointed and a lowered arch are also reported in order to show its performance.

Thrust Line Analysis of masonry arches
A discretized thrust line, or funicular polygon, is represented in Figure 1 Borrowing the nomenclature from the TNA, segments and vertices that form the funicular polygon are called branches and nodes. Branches and nodes are ordered form left to right; hence, branch b j connects nodes n j and n j+1 , while branches b j−1 and b j share the same node n j .
The geometry of the thrust line is described by the coordinates x j = (x j , y j ) of nodes, defined in a two-dimensional Cartesian reference frame. According to Heyman's principles, the thrust line is required to be contained within the thickness of the arch, hence the vertical coordinates of nodes are required to fulfill the inequalities y j,min and y j,max , where y j,min and y j,max are the heights of the arch intrados and extrados. However, different choices can be done when setting y j,min ≤ y j ≤ y j,max , depending on the specific needs of the user. For instance they can be set equal to the thirds of the arch thickness if the limit condition implies a fully compressed arch.
The generic j-th node is loaded by external forces f j = (f jx , f jy ), while the first and last nodes are also subjected to the reactions r l and rr of left and right springers of the arch. Nodes are also subjected to the thrust forces transmitted  by the adjacent branches. In order to model these internal forces, a thrust value t j is associated to the generic j-th branch of the thrust line. Since it represents a compressive force, t j must be positive, i.e. t j > 0, j = 1, ..., B.
Within the TLA framework, the horizontal position of nodes and the external forces are known. Thrusts associated to branches, vertical position of all nodes and the reactions of springers are unknown. In the spirit of the Heyman's safe theorem, these unknowns are computed by employing equilibrium equations only.
The equilibrium of the generic j-th node of the thrust line is where t j−1 j and t j j are the thrust forces that branches b j−1 and b j transmit to node n j . These forces have modulus equal to the thrust values t j−1 and t j associated to branches b j−1 and b j and they have direction parallel to the corresponding branches. Thus, they are computed as Employing previous formulas, equilibrium of node n j becomes Two exceptions shall be considered for the generic equilibrium equation (3). Actually, the first and last nodes of the thrust line are loaded by the thrust forces associated to the first and last branches, respectively, and by the springers' reactions r l and rr. Accordingly, the equilibrium equations of these two nodes read These equations are easily solved for r l and rr as Thus, springers' reactions are easily computed after the thrust values of all branches and vertical coordinates of all nodes are determined. To this end, the set of equilibrium equations associated to internal nodes is solved first.
In order to simplify this set of equations, it is useful to introduce the so-called reference thrustst j = Rt j . They are proportional to the horizontal projectiont j of t j by means of the positive scalar parameter R > 0. Such reference thrusts are introduced within equation (3) by setting where h j = x j+1 − x j . Previous equation has been obtained by invoking the similitude of triangles in Figure 2.
Employing previous formula within Eq. (3), the equilibrium of the j-th internal node becomeŝ This is the generic element of the set of equations used to perform the thrust line analysis.
Notice that the set of equations that is obtained by writing down (7) for all internal nodes is composed of 2(N − 1) = 2N − 2 equations. The number of unknowns, that is represented by B reference thrusts, N vertical coordinates of nodes, and the value of R, exceeds by 2 the number of available equations. Actually, B + N + 1 = 2N. However, these unknowns are also subjected to a set of inequalities introduced before, i.e. t j > 0, y j,min ≤ y j ≤ y j,max and R > 0 .
This means that the set of equations and inequalities can admit infinite solutions, one unique solution or even no solutions. In the first case one is interested in determining two extreme thrust line configurations, respectively associated to a minimum and maximum value of the horizontal component of thrust at springers. Accordingly, these solutions are called solutions of minimum and maximum thrust. Recall that horizontal components of thrusts are inversely proportional to the parameter R, see, e.g., Eq. (6). Hence, these solutions are also associated to a maximum and minimum value of R, respectively. Furthermore, these two solutions correspond to the deepest and shallowest geometry of the thrust line. Hence, they are also called deepest and shallowest solutions.
In case the solution is unique, one can imagine that the interval between the deepest and the shallowest solutions is reduced in such a way that the two extreme cases coincide. This situation represents a limit condition for the equilibrium of the arch since only one equilibrated distribution of internal forces is possible. A small modification of the arch geometry or a modest alteration of the loading condition can result in an unsafe structure.
Finally, in case no solution exists, then there is no distribution of internal forces that equilibrate applied loads and fulfill the material assumptions. This means that the arch is unsafe.
The procedure for solving the set of equations (7), for j = 2, ..., N − 1, and the corresponding inequalities, is particularly simplified in case all nodal forces are vertical. Actually, in this case it is possible to set up a procedure in which thrusts and vertical coordinates of nodes are uncoupled. This simpler case is described first, while the more general case of thrust lines subjected to both vertical and horizontal loads is presented later on.

Solving procedure for thrust lines subjected to vertical loads
In case all nodes of the thrust line are loaded exclusively by vertical forces, the determination of the thrust line configurations of minimum and maximum thrust is relatively simple. It is implemented by the few lines of MATLAB code reported in Figure 3.
In absence of horizontal forces, i.e. when f jx = 0, j = 1, ..., N, the x component of equation (7) readŝ that leads to uniform distribution of reference thrusts in all branches. Accordingly, all reference thrusts can be set equal to an arbitrary small positive valuet 0 . This is done at line 4 in Figure 3. After the reference thrusts are assigned, the algorithm solves twice the vertical equilibrium equations to compute the solutions of minimum and maximum thrusts, respectively. This is done by invoking the function VertEq at lines 7 and 10 of Figure 3. The MATLAB code of this function is reported in Figure 4.
The vertical equilibrium of the generic j-th node is obtained by selecting the y component of equation (7), which readst The set of N −1 vertical equilibrium equations of all internal nodes are written in matrix form as where y is the vector collecting the y coordinates of all nodes, fy collects the vertical component of the forces applied at nodes n 2 ..., n N−1 and D is the N − 2 × N matrix collecting the reference thrust densitiest j /h j of all branches. It  is obtained by assembling the contribution of each branch, that is It is assembled in rows j and j + 1 and columns j and j + 1 of D (lines 8 − 18 of Figure 4). The unknowns y and R are determined by solving Eq. (11) together with the set of inequalities that express the corresponding constrains. Hence, y and R are computed by solving the following constrained linear optimization problem: Notice that the objective function in (13) is either −R or +R, depending whether the optimization problem is used to maximize or minimize R. Either one of the two cases is selected by setting pm = -1 or pm = +1 to the last input variable of VertEq (lines 7 and 10 in Figure 3). Actually, from equation (6) it is clear that higher values of R are associated to lower values of thrust and vice versa. Accordingly, when the objective function −R is used, the value of R is maximized and the optimization procedure returns the so-called solution of minimum thrust or deepest solution. Similarly, the solution of maximum thrust or shallowest solution is obtained when R is minimized, i.e. when the objective function is set as +R.

Solution procedure for thrust lines subjected to horizontal and vertical loads
When nodes are loaded by horizontal and vertical forces, both the horizontal and vertical components of the equilibrium equations of nodes contain the unknown R. Hence, the horizontal and vertical equations are coupled and cannot be solved in sequence. Thus an iterative procedure that solves repeatedly the horizontal and vertical equilibrium of the thrust line needs to be adopted. The corresponding MATLAB code is reported in Figure 5 : lines 3 − 18 compute the deepest solution and lines 20 − 35 the shallowest one.
Similarly to the case of null horizontal forces the only difference between these two portions of code regards the sign of the objective function used to solve the vertical equilibrium equations of nodes. The mentioned iterative procedures start by choosing an arbitrary tentative value of R so that the horizontal equilibrium equation can be solved for the reference thrustst j of nodes. This is done at lines 4 and 21.
For any iterative estimate of R, including this first guess, the horizontal equilibrium equation of the generic j-th node becomeŝ Accordingly, the reference thrustt j associated to a generic branch can be computed as a function of the reference thrustt j−1 of the previous branch. Hence, after choosing an arbitrary value of reference thrust for the very first branch, all successive thrusts can be evaluated in sequence. Of course, this solution does not necessarily fulfill the constraint t i > 0, i = 1, ..., B. This issue is easily solved by modifying the tentative distribution of reference thrusts by adding t i min to all thrust values. This set of positive reference thrusts that fulfill the horizontal equilibrium equations and the relevant constraints are computed by invoking the function RefThr (lines 8 and 25 in Figure 5). The corresponding MATLAB code is reported in Figure 6.
Once reference thrusts are assigned, it is possible to evaluate the nodal heights y and a new estimate of R by solving the constrained linear optimization (13). This is done by invoking the function VertEq (lines 12 and 29 of Figure 5), which has been already described in the previous section.
This new estimate of R is then used again to evaluate a new estimate of reference thrusts by repeating the same set of operations. The procedure is iterated until the relative difference between two successive estimates of R is lower than a given tolerance, that is Such a convergence check is coded at lines 15 and 32 of Figure 5.

Numerical examples
Two numerical examples are reported hereafter, regarding the thrust line analysis of both a pointed arch and a lowered circular arch. Both problems are solved by the code contained in the MATLAB scripts ExampleArch, also available from the link reported earlier. The function ArchLab implements the TLA by the algorithm described in section 2, the two table files PointedArch.csv and LoweredArch.csv contain the input data for the examples regarding the analysis of the lowered and pointed arches and the plotting function PlotArch is used to graphically represent the solutions.
The function ArchLab takes as input the horizontal position x of the nodes of the thrust line, their height limits, y min and ymax, and the applied forces, fx and fy. These data are included in the two table files PointedArch.csv and LoweredArch.csv, respectively. The geometry and the loading condition of these arches is determined as a function of a few geometric and loading parameters that are (Figure 7): the arch span S, rise H and thickness t; its unitary weight Fsw; uniformly distributed load Fu; the height Y f and unitary weight of filling F f .
Geometric and loading parameters of the two models are reported in Table 1. Being planar models, unitary weights and distributed loads are assigned as force per unit area and force per unit length, respectively. In order to stress out the algorithm, the thrust lines of both arches has been generated by using a horizontal spacing between nodes equal to one hundredth of the span. This resulted in a thrust model composed of 122 nodes and 121 branches.
The intrados of the pointed arch is the composition of two circles of radius D/2 = H 2 /S + S/4 ≈ 2.08 m. Their centers are positioned at the height of the arch springers and distant D/2 − S/2 ≈ 0.58 m on both sides of the axis of symmetry. The arch is loaded by horizontal forces which are pro-    Figure 9. Also in this case, the even lower value of the ratio Rs /R d = 0.18 shows the capability of lowered arches to withstand higher values of horizontal forces.

Conclusions
The function ArchLab is used to determine thrust line configurations of minimum and maximum thrust of masonry arches subjected to both vertical and horizontal loads. It represents a specialization of the Thrust Network Analysis to the case of a planar line of thrust. The implemented algorithm has been described in full detail in section 2 where a precise reference to the implemented MATLAB code is also reported. When only vertical loads are applied, the algorithm is capable to determine the thrust line configurations without the need of any recursive method. Conversely, when also horizontal loads are applied, the procedure re-F. Marmo quires iterations. In both cases, the algorithm computes the horizontal components of reference thrust in closed form at each iteration, while the height of nodes is determined by solving a constrained linear optimization problem by using the MATLAB builtin implementation of the Dual-Simplex algorithm.
By implementing additional few lines of code, not commented in this paper for brevity, it is possible to employ the same functions described earlier to solve a form-finding problem. It is sufficient to set y min = 0 and ymax = H for all nodes of the model and compute the deepest configuration of the thrust line. Here H represents the arch axis height. An iterative procedure can be used to update the self-weight of each branch of the thrust line and compute updated nodal forces accordingly. An example of such an iterative algorithm is coded within the script ExampleFF, which finds the optimal shape of an arch of span S = 8.0 m and rise H = 3.0 m, subjected to its own self-weight. The solution is obtained after only 4 iterations and is then plotted by using the function PlotFF. It is reported in Figure 10, where it is compared with the catenary curve y(x) = H + a − Cosh(x/a), with a = 3.0668 m, showing perfect agreement.

Conflict of interest:
The author states no conflict of interest