A digital tool to design structurally feasible semi-circular masonry arches composed of interlocking blocks

This work deals with a digital tool to design stable semi-circular masonry arches composed of interlocking blocks which are kept together by interlocking connectors on their faces. These blocks, comparing to conventional blocks, increase the sliding resistance and reduce the workmanship. However, the digital tools were developed mostly to design arches with conventional blocks. The proposed tool tries to fill this gap by addressing the work in three stages. First, a heuristic method is developed to define the relationships between the geometry of an interlocking face and the sliding resistance. Then, a structural analysis procedure is developed based on the limit analysis and the heuristic method to define the stability condition of the arch. Finally, optimization algorithms are developed to find the thinnest arch by means of two minimization strategies dealing with the relationship between the sliding resistance of the blocks and the geometry of the interlocking faces, differently. The algorithms consider some control points on a given thrust line and automatically adjust them to minimize the thickness, while the stability condition checks the structural feasibility during the geometry adjustment. To evaluate the accuracy of the proposed heuristic method, the results obtained with nonlinear FE analysis are used for comparison.


INTRODUCTION
This work develops a digital tool to design structurally sound semi-circular masonry arches composed of interlocking blocks. The interlocking blocks are rigid block units which, on their faces, have connectors keeping the blocks together and preventing blocks from sliding.
The interlocking blocks, in comparison to the conventional blocks with flat faces, increase the shear resistance to external forces [1] and reduce the construction skills and instruments required [2]. Despite these advantages of interlocking blocks, the digital tools supporting designers for structurally informed architectural design were developed mostly for the conventional blocks. Most of those efforts (reviewed by Rippmann [3]) were focused on designing single-layer masonry vaults. To design masonry assemblages with diverse topologies, Whiting [4] proposed a method in which, given a block assemblage, the block geometry changes automatically to obtain the structurally feasible form. In that work, however, interlocking the blocks during the geometry modification was avoided. There are a few studies [5,6] that applied simple geometric grammars but no structural constraints to design vaults that are composed of interlocking blocks with limited geometries.
In the literature there are different methods of structural analysis to develop the introduced digital tools. Most of them are based on the limit analysis theory aimed at calculating the ultimate load factor satisfying static and kinematic conditions. According to this theory, internal forces are distributed at the interfaces between rigid blocks. These interfaces have infinite compressive, no tensile strength, and infinite [7] or finite sliding resistance including associative [8,9] or non-associative [10][11][12][13] frictional resistance. In this framework, it is worthy to mention a recent numerical limit analysis procedure, named discontinuity layout optimisation (DLO) [14], which has been developed to obtain accurate upper-bound solutions for plane-strain collapse problems. With reference to vaulted single-layered masonry structures, different methods find the stress states at interfaces: line of thrust methods [7,15,16]; membrane theory [7,17,18]; force network method [19,20]; convex and concave contact formulations [9,21], non-smooth contact dynamics [22]. Interesting are also some recent experimental and/or analytical works on the collapse failure of vaulted masonry structures under horizontal loading [23,24] and after actual displacements of the supports [25][26][27]. However, all these methods were applied to analyze structures composed of conventional block with isotropic friction.
Other distinctive methods of masonry analysis, such as discrete and finite element analysis, rarely have been used to develop tools to design structurally feasible masonry assemblages, mostly due to their complex computational techniques. In the former method, masonry is modelled as an assemblage of distinct blocks. For such a block system, the equations of motion are solved using an explicit time integration method [28][29][30]. Instead, in finite element analysis (FEA), masonry is modelled as connected or distinct elements with different material properties for bricks and mortar joints (detailed and simplified micro models) or as a continuum by smearing the bricks and mortar joints into an isotropic or anisotropic homogeneous continuum (macro models) [31][32][33].
This work aims at filling the gaps described above about design of structurally sound assemblages of interlocking blocks and proposes a digital framework for the structurally informed design of a semi-circular arch composed of interlocking blocks and subjected to its own weight. To develop this framework, limit state analysis is adopted and extended to interlocking blocks. This structural typology with symmetrical loading and geometry belongs to a special class of non-associative friction problems for which provably unique solutions within limit analysis approach exist [16]. Therefore, the application to this structure allows focusing the attention to the structural behaviour of interlocking blocks rather than to the issue of associated or non-associated flow rules, which is beyond the scope of this work.
Three main stages are herein introduced to help the designer to evaluate the structural feasibility of a model and automatically update the model to make it structurally optimal.
1. The first stage of this work is to develop a novel heuristic method, assisting the designers to estimate the interlocking block resistance to sliding. The proposed method allows equating an interlocking block to a conventional block with different sliding resistances along two orthogonal directions (normal and parallel to the connectors).
2. The goal of the second stage is to present a new structural analysis procedure that, applying the heuristic method proposed above, checks if the semi-circular arches composed of interlocking blocks are stable. To achieve this goal, this work extends the limit analysis developed by [16] which applied the line of thrust method to analyze the stability of semi-circular arches composed of blocks with isotropic finite friction. The extension deals with semi-circular arches composed of interlocking blocks which, using the heuristic method, are modelled as conventional blocks with orthotropic finite friction (in which the sliding resistances are the minimum and maximum in two orthogonal directions, respectively), in order to analyze their structural stability.
3. Finally, the third stage of this research is to develop an optimization method minimizing the structurally feasible semi-circular arch thickness. The structural feasibility of the model during optimization is checked by the introduced structural analysis procedure.
Considering the relationship between the block resistance to sliding and the geometry of interlocking faces, the optimization method offers two strategies to minimize the arch thickness: (1) updating the sliding resistance due to the connector geometry change during optimization; (2) changing the number of the block connectors (initially given by the designer) during the optimization process, to keep the sliding resistant constant.
The paper is organized as follows. The structural modelling along with the heuristic method introduced above are presented in Section 2, while Section 3 performs the extension of limit state analysis to interlocking blocks. Section 4 develops two optimization strategies to design semi-circular arches with interlocking blocks. Results are presented, discussed and validated in Section 5. Finally, the conclusions are summarized in Section 6.

Structural model of the interlocking block arch
The structural model adopted in this work is based on the assumption that masonry structures are composed of assemblages of discrete rigid interlocking blocks in contact interaction, to be constitutively defined, where the displacements of each block should be considered separately. How a semi-circular arch and each of its discrete interlocking blocks are modelled in the developed framework is described as follows.
Within the arch modelling, the designer assigns the radius R for the semi-circular arch centreline, the arch width b, and the number of blocks m composing the arch (Figure 1a). The latter number should be odd to model the arch keystone. All arch blocks have same sizes.
On the other hand, the block modelling is based on the assumption that each interlocking block has two corrugated and two flat faces. Interlocking faces lock blocks to form the arch and the corrugated faces of the connectors are assumed to have rectangular cross sections ( Figure 1b). Different shapes of the connectors could be investigated, but this further analysis will be addressed in future work. To model the interlocking faces, the designer should specify the total number n of projections and depressions of the block. This number should be an odd number in order to guarantee the symmetry of the block shape with respect to the centreline of the arch.
Within the introduced discrete approach, the block interfaces are treated as the elements of the problem while the blocks are simply defining the geometry of the problem. Therefore, the analysis is fully related to the behaviour of the interfaces, which can then be regarded as systems of stresses. Considering finite friction, these stresses are governed by unilateral contacts and constraints, due to the rugged projections on their faces, termed asperities, they can be regarded as connectors of two interlocked blocks.
Using corrugated interlocking faces minimum along the connectors and a maximum along the direction normal to the connectors ( Figure 1b). In fact, the friction on the equivalent flat face can be regarde friction, in which the sliding resistances are the minimum and maximum in two orthogonal directions, respectively. The next section presents a method to analyse the sliding resistance of an interlocking block and to relate the sliding resi frictional resistance of conventional blocks. For example, the interlocking sliding resistance respectively section is aimed at defining a heuristic relationship between the sliding resistance of interlocking blocks and that of the Finding such a relationship, limit state analysis method structural behaviour of a semi The first step for searching out direction orthogonal to the connectors of interlocking reference to a single interface between a fixe given normal force N and a lateral load 6 Considering conventional blocks with rough flat faces in contact , these stresses are governed by unilateral contacts and constraints, due to the rugged projections on their faces, termed asperities. asperities, they can be regarded as connectors of two interlocked blocks.
Using corrugated interlocking faces, the block resistance to sliding at the interface is a minimum along the connectors and a maximum along the direction normal to the connectors ). In fact, the friction on the equivalent flat face can be regarde in which the sliding resistances are the minimum and maximum in two orthogonal directions, respectively. The next section presents a method to analyse the sliding resistance of an interlocking block and to relate the sliding resistance of interlocking blocks to the frictional resistance of conventional blocks. e interlocking block adopted in this work has the maximum respectively along the direction normal and parallel to the section is aimed at defining a heuristic relationship between the sliding resistance of interlocking blocks and that of the equivalent conventional blocks with rough relationship, limit state analysis method is then extended structural behaviour of a semi-circular arch composed of interlocking blocks ( The first step for searching out this relation is to study the sliding resistance along the direction orthogonal to the connectors of interlocking blocks. The analysis is carried out with reference to a single interface between a fixed lower block and the upper one subjected to a and a lateral load T applied to its centre of gravity ( Figure 2 conventional blocks with rough flat faces in contact with , these stresses are governed by unilateral contacts and isotropic frictional asperities. By scaling up the , the block resistance to sliding at the interface is a minimum along the connectors and a maximum along the direction normal to the connectors ). In fact, the friction on the equivalent flat face can be regarded as orthotropic in which the sliding resistances are the minimum and maximum in two orthogonal directions, respectively. The next section presents a method to analyse the sliding resistance stance of interlocking blocks to the a composing interlocking block.

Heuristic formulations relating the sliding resistance of interlocking blocks to the
The interlocking blocks could have different sliding resistances in different directions, depending on the orientation and the mechanical and geometrical features of the connectors.   Since the failure mode is first governed by the shear collapse of the projections, the ultimate lateral force strongly depends on the amounts of the two resistances of Eq. (3). In fact, if T R > T C the ultimate lateral force will be frictional resistance is not enough to prevent sliding. This condition, implying non behaviour of the interface, is represented by the black dot in Figure   constant shear resistance is superposed on the choesionless Coulomb yield condition. On the contrary, if T R < T C the lateral force will be governed by may occur because at the onset of the shear failure mode of projections the frictional resistances will also be activated and the sliding is prevented until the shear achieves the value T C (Figure 5b). This means that, to avoid sliding of an interlocking block connectors, the following inequality should be always met: 9 block. In this case, the shear resistances of each upper and lower T 0 , respectively. Whatever the position of projections i) or ii) in , the resultant limiting shear force of all connectors and the frictional resistance is the friction coefficient, N is the fixed normal force (including the own weight of . Sliding resistances of interlocking blocks with upper connectors higher than lower ones . Sliding resistances of interlocking blocks with lower connectors higher than upper ones (Case 3) Since the failure mode is first governed by the shear collapse of the projections, the ultimate trongly depends on the amounts of the two resistances of Eq. (3). In fact, if the ultimate lateral force will be T = T R because once the failure is activated the frictional resistance is not enough to prevent sliding. This condition, implying non behaviour of the interface, is represented by the black dot in Figure 6a, where the greater superposed on the choesionless Coulomb yield condition. On the the lateral force will be governed by T C (T = T C ) and a ductile behaviour may occur because at the onset of the shear failure mode of projections the frictional be activated and the sliding is prevented until the shear achieves the that, to avoid sliding of an interlocking block in the direction normal to the , the following inequality should be always met: In this case, the shear resistances of each upper and lower , respectively. Whatever the position of projections i) or ii) in , the resultant limiting shear force of all connectors and the frictional resistance of all (3) is the fixed normal force (including the own weight of higher than lower ones (Case 2). Besides, the sliding in the direction parallel to the connectors is only governed by the inequality ܶ ߤܰ. This means that the interlocking block proposed by this paper can be considered as a conventional block with orthotropic friction so that the frictional constraints are Ineq. (5) or ܶ ߤܰ, depending on the direction.
The heuristic formulation given by Ineq. (5) will be applied in the next section to the limit state analysis of a semi-circular arch composed of interlocking blocks and subjected to its own weight. In this case, only the sliding resistance along the direction orthogonal to the connectors of interlocking blocks is activated. Later in Section 4, this heuristic formulation will be implemented to develop the optimization algorithm for such structures.

INTERLOCKING BLOCKS
As introduced in Section 1, standard limit analysis of rigid block assemblages has been found to be a valuable computational tool for developing two groups of structural optimization problems: 1) to find the maximum load a structure can tolerate (predict the collapse load); and 2) to find the minimum material (i.e., minimum thickness) a structure can have to bear a specific amount of load. Both approaches lead to find the classical rocking failure mechanism of such structures.
To develop these optimization problems, either static or kinematic theorems can be applied.
The occurrence of Coulomb frictional sliding, which implies non-associative flow rule, is in general excluded from the analysis in order to ensure the validity of the normality condition, one basic hypothesis of classic limit state analysis. In fact, it is well-known that the bounding theorems of plastic limit analysis do not in general provide unique solutions for the collapse load factor if a non-associative flow rule is specified [34] and the classical procedure does not assure that the structure is safe.
However, Casapulla and Lauro [16] have identified a special class of non-associative friction problems for which provably unique solutions exist. The class comprises arches with symmetrical loading and geometry. The proposed procedure was applied to arches of this sort to both verify that the numerical and analytical solutions coincide and to investigate the convergence characteristics of the method.
Exploiting the uniqueness of the solution due to the symmetry, the present work is mainly aimed at extending this optimization method to design the minimum thickness of a semicircular arch composed of interlocking blocks and subjected to its own weight.
In fact, to minimize the thickness of such an arch, first the desired radius for the arch centreline R is specified and then the optimization problem tries to find the closest thrust line to this centreline which satisfies the equilibrium and yield conditions at contact interfaces.
This section explains how to find a thrust line meeting equilibrium and yield conditions. Next section introduces two optimization strategies to find the optimal conditions. The general thrust line is defined by a set of control points variable in function of two points A(0; Y A ) and B(X B ; 0), respectively at the crown and the springing joints ( Figure 7a). The X-coordinates of the control points are the same of the application points of ea while the Y-coordinates are found H tot can be found as: where W tot is the weight of the half mass. Being W i the weight of block of the control point i, ∆ܻ can be achieved by A valid thrust line should also meet the heuristic method described in Section 2.2, interlocking blocks composing a semi can be considered as conventional blocks whose sliding resistance can be evaluated by A valid thrust line should also meet the yield conditions at contact interfaces. A heuristic method described in Section 2.2, interlocking blocks composing a semi can be considered as conventional blocks whose sliding resistance can be evaluated by can be rewritten as: are the tangential and normal components of S i , respectively, According to Casapulla and Lauro's method [16], for the case of semi symmetric loads, the ultimate value for T i is independent of coordinates of the control points are the same of the application points of each block weight, , the horizontal reactive thrust (7) coordinate of the half-arch centre of coordinate both of its application point and the equilibrium condition: can be achieved by: yield conditions at contact interfaces. Applying the heuristic method described in Section 2.2, interlocking blocks composing a semi-circular arch can be considered as conventional blocks whose sliding resistance can be evaluated by Ineq.
, respectively, (Figure 7b After this substitution, Ineq. (10) can be splitted into the two following yield conditions for semi-circular arches composed of interlocking blocks: Lastly, it should be taken into consideration that since all block weights, the thrust line, and the limiting sliding forces for all blocks lie on the same plane, the problem can be considered as a 2d problem. In this case, the sliding resistance of all interfaces with orthotropic friction is only regarded to be the maximum sliding resistance value on those interfaces (orthogonal to the connectors in Fig. 1b). This means that an arch with interlocking blocks which embeds a thrust line whose control points are obtained by Eqs. (7) and (9) and meets Ineqs. (11) and (12), is a stable arch.

STRUCTURAL OPTIMIZATION
As explained earlier, the optimization problem of this paper is aimed at finding the closest thrust line to the specified arch centreline (minimum thickness) which satisfies the equilibrium and yield conditions at interfaces.
Unlike the conventional blocks whose frictional properties remain fixed during optimization, changing the arch thickness during optimization may change the sliding resistance at interfaces. In the following sub-sections, two strategies are proposed to minimize the arch thickness (Figure 8), based on different relationships between the geometry of the interlocking interfaces and the sliding resistance. In the first strategy, the sliding resistance changes during the optimization process; in the second strategy, the geometric properties of interfaces changes in a way that the sliding resistance remains fixed during optimization.
For both strategies the objective function is the minimization of the maximum value of the arch half-thickness (a min /2), as described in [16] for the arch with conventional blocks. The variables of the problem are X B and Y A of the generic thrust line (Figure 7). The optimal arch only embeds the obtained thrust line indeterminate structure into a mechanism Section 5.

Second strategy
In this strategy, shear resistance can be obtained. In should be at least three, and considering the arch initial /2) to be the maximum distance between the radius and control points on can be calculated as follows: where x i t0 is the same as x i t of the objective function (13).
This value remains fixed during optimization, while both arch thickness a and n change during optimization. Given n 0 , x i t0 and y i t0 , the relation between a (which equals 2ቚඥ(‫ݔ‬ ௧ ) ଶ ‫ݕ(‬ ௧ ) ଶ − ܴቚ) and n can be defined as follows: which should be an (a) integer, (b) odd, and (c) positive number during optimization. As a result, the optimization problem for the second strategy is as follows: This objective function of the optimization problem is the same as the previous strategy and the sliding resistance of the arch is met by applying the two first inequalities.

RESULTS AND EVALUATIONS
This section is composed of two main parts. First, a case study of a semi-circular arch of 10m centreline radius, containing 27 discrete blocks is analyzed by means of the two introduced optimization strategies. These results are investigated to compare the two proposed optimization strategies to each other and to the Casapulla and Lauro's optimization method developed for arches with conventional blocks [16]. Then, to evaluate the accuracy of the proposed heuristic method and the limit state analysis method, the interlocking block interface and the case study are analysed by Finite Element method for comparison.

Case study
In the following, after a brief description on the implemented process, presented to investigate the relation between the minimum two introduced optimization strategies Finding these relations, minimum thickness/Radius (t/R) ratio can be defined as a function of the parameters determining the block sliding resistance.

Implementation
For modelling the arch with calculations, Grasshopper's C programming language which runs within Rhinoceros 3D.
calculations and optimization are performed by MATLAB used as backend.
optimization, MATLAB's fminimax method was used. Multiple hard constraints in the second strategy imposed on variable constraints were used, weighted

Optimal arches with conventional
This sub-section uses the results obtained by the sake of comparison with applied to conventional blocks minimum arch thickness required for stability and the coefficient of friction introduced semi-circular arch It results that: 1) for µ > 0.395 centreline radius is approximately reproduced with this discretization (10.74%) failure mode is governed by a pattern of hinges at the intrados and extrados 0.395 ≥ µ ≥ 0.332, the failure mechanism is characterized by the sliding interfaces at the horizontal supports along with µ < 0.332 the equilibrium is no longer possible  The results obtained from the first strategy can also be represented from a larger perspective by Figure 11, for µ between 0.3 and 0.4

Optimal arches with interlocking blocks obtaine strategy
According to the second strategy proposed above, the initial value of the shear resistance is

Optimal arches with interlocking blocks obtained by the second optimization
According to the second strategy proposed above, the initial value of the shear resistance is arch thickness and the number of connectors change during optimization.

d by the second optimization
According to the second strategy proposed above, the initial value of the shear resistance is number of connectors change during optimization.
) is highly dependent on the initial due to the way that the final n is related to (17) is obtained when the initial value of shear resistance T R 0 , which should Considering the final n the initial thickness a 0 which is be underscored that the minimum thickness obtained by the second strategy might differ from that of the first strategy for the same values of friction coefficient is kept unchanged during the optimization of the second The right hand side of this equation is a known value while that this value is achieved. Both equilibrium condition with sliding constraints and result, a min satisfying Eq. (19) It is concluded that the thrust line which is found as the optimal result solution of the arch embedding arches cannot be precisely found can be found as follows.
For example, Figure 13 depicts the tangential force at interface by the first strategy when this interface reaches the limiting sliding value (blue curve). This  The right hand side of this equation is a known value while a min and n should be chosen so that this value is achieved. Both a min and n are constrained variables; a min equilibrium condition with sliding constraints and n should be a positive odd value. As a ) might be different from results of first optimization strategy.
It is concluded that the thrust line which is found as the optimal result may embedding it. In other words, the rocking failure mechanism of the cannot be precisely found. Still, the closest possible solution to that failure mechanism depicts the tangential force at interface m of the case study optimized by the first strategy when this interface reaches the limiting sliding value (blue curve). This depicts the tangential force at interface m of the case study ( 10.81 m) optimized by the second strategy (red curve) for the same . As shown, these two curves are quite close but not identical. In fact, blue curve represent the closest possible value to the limiting one satisfying Eq. (17) . Difference between tangential force at interface m for the optimal arch optimized by first and second optimization strategies two discussed issues, the relationship between the minimum arch thickness for the arch explained above, with n 0 = 3 and X  and Eq. (19), yellow continuous curve in Figure 1 different from that of Figure 12 (results of the first strategy) which is redisplayed by purple dash curve. Still, their overall configuration is similar, as explained above.

Validation and calibration via FEA
The heuristic method and the extension of limit analysis of this paper are validated by comparing the results obtained by the proposed methods with the results achieved by Finite Element Analysis (FEA) based on the simplified macro modelling technique proposed by [32]. In this method, blocks are morphed (controlled expansion) to eliminate mortar joints and are considered to be solid or shell elements, while mortar joints are zero-thickness interface elements. This validation includes three parts: in the two first parts the proposed heuristic method is studied by analysing two sets of interlocking blocks with different geometries; the third part studies the proposed limit analysis methods and optimization strategies by analysing the case study of the arch analysed above.
The FEA is developed in ABAQUS 6.11-1 by adopting the brittle cracking model. In the first two parts, to purely focus on the shear behaviour at the joint between the main body and the connectors of a block, the following model is designed and analysed: two stacked blocks with a shared interlocking interface are modelled. The connectors and the main body of each of the blocks are modelled with rigid body elements (rigid parts) while the joint between them is modelled with C3D8R elements (flexible part with brittle behaviour) (Figures 16b and 18b).
The material properties of elements are shown in Tables 1. The effective density of the expanded blocks is calculated using homogenization process proposed in [32]. The lower face of the lower block is bounded so that no degree of freedom is allowed. The upper face of the upper block is bounded so that one degree of freedom is allowed in the direction of lateral force application. The loads imposed on the model are the weight of the upper block and a lateral force distributed on a vertical face of this block (Figures 16c and   18c). The lateral force by which the principal stress reaches the tensile strength is found. This value is compared to the shear resistance obtained by the proposed heuristic method. The distribution of principal stresses and shear stresses τ xy in X direction on the plane whose normal vector is Y is studied as well.
model. Figure 16a shows the size of the each block.

Analysis of interlocking blocks with
The maximum lateral force obtained by FEA is 9.42 kN, when the maximum principal stress (0.36 MPa) is observed on the joint between the connector and main body of the lower block.

Analysis of the case study arch composed of
In the case of the arch with assigned radius are modelled with C3D8R elements and the mortar joints between them are modelled as shell interfaces with friction coefficient In fact, using the proposed limit state method, for this mechanism is observed, that is fully captured by the FE model, as shown in Figure 20.
In sum, the three cases examined shear resistance is less than the frictional resistance, the arch behaves as an arch with . When the shear resistance is larger than the frictional resistance, the structural behaviour of the arch is governed by the shear strength. For small values of shear strength, the combined rocking/sliding mechanism occurs for the optimal result. For large values of shear strength, pure rocking mechanism can be observed on the optimal result.
In the second method, the sliding resistance is kept fixed via changing the number of the connectors. The results are mostly dependent on the constraints keeping the number of the connectors an odd integer value and also on the initial value of sliding resistance defined by the designer. Due to these items, the optimal result for the minimum thickness might be thicker comparing to the result of the first optimization method, given the same shear strength, friction coefficient and geometric inputs. In this case, the mechanism does not occur on the optimal solution.
A good agreement of the results obtained by the heuristic method and optimization procedure was found in case of pure shear conditions, while some limits of validation were highlighted in case of mixed rocking/sliding failure modes. Future works will address these issues together with the extension of the approach to the 3d models. Furthermore, the sliding resistance of interlocking blocks will also be experimentally investigated together with the analysis of different shapes of the interfaces.

ACKNOWLEDGMENTS
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No.
791235. It reflects only the authors' view and the Agency is not responsible for any use that may be made of the information it contains.