Static Upper/Lower Thrust and Kinematic Work Balance Stationarity for Least-Thickness Circular Masonry Arch Optimization

Abstract

This paper re-considers a recent analysis on the so-called Couplet–Heyman problem of least-thickness circular masonry arch structural form optimization and provides complementary and novel information and perspectives, specifically in terms of the optimization problem, and its implications in the general understanding of the Mechanics (statics) of masonry arches. First, typical underlying solutions are independently re-derived, by a static upper/lower horizontal thrust and a kinematic work balance, stationary approaches, based on a complete analytical treatment; then, illustrated and commented. Subsequently, a separate numerical validation treatment is developed, by the deployment of an original recursive solution strategy, the adoption of a discontinuous deformation analysis simulation tool and the operation of a new self-implemented Complementarity Problem/Mathematical Programming formulation, with a full matching of the achieved results, on all the arch characteristics in the critical condition of minimum thickness.

1 Introduction

This work further investigates the issue of (symmetric) circular masonry arch form optimization (Couplet–Heyman problem), in the quest of a least-thickness evaluation under uniform self-weight (Figs. 1 and 2). The modern framing of such a problem relies in the contemporary contributions by Jacques Heyman [1,2,3,4,5,6] and the recent revisitation in earlier companion work [7], with therein extensive references, also to various historical and development perspectives on the subject.

Fig. 1

Sketch of a symmetric circular masonry arch subjected to self-weight (of specific weight per unit volume γ) with all involved characteristic parameters (d: out of-plane depth of the arch)

Fig. 2

Five-hinge purely rotational collapse mode, with abutment-hinge reactions V and H (opposite of weight of half-arch and horizontal thrust)

The present investigation belongs to a research project by the authors on the statics of masonry arches [7,8,9,10,11,12,13,14], where the following treatments have been attempted: analytical [7, 14]; analytical–numerical, accompanied by a Discrete Element Method (DEM) investigation, through an available Discontinuous Deformation Analysis (DDA) tool [8, 9, 15,16,17], and including reducing friction effects and resulting mixed collapse modes [10, 11]; analytical–numerical, by an innovative Complementarity Problem/Mathematical Programming formulation, truly accounting for finite friction implications [12, 13]. This shall be framed within the relevant, updated literature, specifically considering the issue of minimum thickness evaluation [18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42], and dedicated attempts to reveal further and independent/complementary information on the analytical treatment, with additional separate validation in terms of numerical results.

Specifically, in the least-thickness collapse evaluation, classical Heyman solution [3] is shown to constitute a sort of approximation of the true solution (here labelled as “CCR” [7]), in Heyman assumption of self-weight distribution along the geometrical centreline of the arch, while Milankovitch solution [43,44,45], as a cornerstone of thrust-line-like analysis, in view of form optimization [46,47,48,49,50,51,52], may as well be derived, in the consideration of the real self-weight distribution along the arch, though at the price of a recorded increasing complexity in the explicit analytical handling of the governing equations (now analytically resolved to a very end in [14]).

The analysis makes reference to the classical three Heyman hypotheses of masonry structure behaviour (no tensile strength; infinite compressive strength; no sliding failure) and refers just to the potential development of a purely rotational collapse mode (as at infinite friction). Given the value of half-opening angle α of the (symmetric) circular masonry arch, the following three basic arch characteristics are sought, in the least-thickness condition of incipient collapse: angular inner-hinge position β from the crown; thickness t to radius r ratio η = t/r and non-dimensional horizontal thrust within the arch h = H/(w r), w = γ t d where H is the horizontal thrust, w = t d the specific weight per unit length of geometrical centreline of the arch, γ and d the constant specific weight per unit volume and out-of-plane depth of the arch. Alternatively [7, 14], one may also make reference to intrinsic non-dimensional horizontal thrust ĥ = η h = H/(γ d r²), defined at given material (γ) and geometrical properties (d, r) of the circular masonry arch to be optimized (critical η still to be sought).

The paper is organized as follows. First, Sect. 2 provides a basic framing, with all the main governing equations, specifically concerning equilibrium and tangency conditions, for the line of thrust (locus of pressure points) within the least-thickness arch. Then, Sect. 3 derives alternative analytical solution approaches to deliver “correct” CCR solution [7] vs. “approximate” Heyman solution [3], through a “Coulomb’s static approach” based on the so-called upper and lower horizontal thrusts, and a “Mascheroni’s kinematic approach”, based on the work/power balance at incipient collapse, with least-thickness condition consistently stated within that. A representation of “accurate” Milankovitch solution [43,44,45] is then also recalled, and the response of the mechanical system illustrated, in terms of the three achieved solutions. Subsequently, Sect. 4 develops a further validation part by a separate numerical treatment, where: first, a recursive determination of angular inner-hinge position is developed; second, a final DDA validation is deployed, completing that earlier presented in [8, 9]; third, an original self-implemented Complementarity Problem/Mathematical Programming computational strategy is adapted and operated, to deliver the final arch characteristics in terms of all kinematical (β), geometrical (η) and statical (h, $\hat{h}$) quantities. A full matching between numerical results and analytical outcomes is recorded, in showing the validity of the three solution instances in terms of least-thickness masonry arch form optimization and relevant collapse evaluation. Finally, some concluding remarks and perspectives of the study are outlined in closing Sect. 5.

2 Basic Mechanical System Equations

2.1 Equilibrium Condition

At incipient (assumed) purely rotational collapse, in the least-thickness condition, the equilibrium of the (symmetric) circular masonry arch (thus, of the half-arch) shall be imposed. From the rotational equilibrium of upper portion AB around inner-hinge B, one has (Fig. 3):

$$ H \cdot \left[ {r + \frac{t}{2} - \left( {r - \frac{t}{2}} \right)\cos \beta } \right] = W_{1} \cdot \left[ {\left( {r - \frac{t}{2}} \right)\sin \beta - x_{1} } \right] $$

(1)

where, in Heyman assumption of a uniformly distributed self-weight along the geometrical centreline of the arch, weight W₁ of the upper portion of the half-arch is obtained as

$$ W_{1} = \int_{\,\,0}^{\,\beta } {\underbrace {\gamma td}_{w}\;r\;d\varphi } = w\,r\,\left[ \varphi \right]_{\,0}^{\,\beta } = w\,r\,\beta $$

(2)

and the centre of gravity of the upper portion of the half-arch is located at following horizontal distance x₁ from the vertical axis of symmetry (and vertical distance y₁ from centre O):

$$ x_{1} = \frac{{\int_{{\;V_{1} }} {\gamma \,x\,dV} }}{{\int_{{\;V_{1} }} {\gamma \,dV} }} = \frac{{\int_{\,0}^{\,\beta } {\gamma td\,r\;rsin\,\varphi \,\,d\varphi } }}{{\int_{\;0}^{\,\beta } {\gamma td\,r\,d\varphi } }} = r\;\frac{{\left[ {\cos \varphi } \right]_{\beta }^{\,0} }}{{\left[ {\varphi } \right]_{\,0}^{\,\beta } }} = r\;\frac{1 - \cos \beta }{\beta } $$

(3)

$$ y_{1} = \frac{{\int_{{\;V_{1} }} {\gamma \,y\,dV} }}{{\int_{{\;V_{1} }} {\gamma \,dV} }} = \frac{{\int_{\;0}^{\,\beta } {\gamma td\,r\;r\cos \varphi \,\,d\varphi } }}{{\int_{\;0}^{\,\beta } {\gamma td\,r\,d\varphi } }} = r\;\frac{{\left[ {\sin \varphi } \right]_{\,0}^{\,\beta } }}{{\left[ \varphi \right]_{\,0}^{\,\beta } }} = r\;\frac{\sin \beta }{\beta } $$

(4)

By substituting Eqs. (2–3) into rotational equilibrium Eq. (1), and shifting to non-dimensional quantities η = t/r and h = H/(w r), one obtains a first equilibrium relation, (for instance) in terms of h:

$$ h = h_{\,1} = \frac{{\left( {2 - \eta } \right)\,\,\beta \;sin\,\beta - 2\left( {1 - \cos \beta } \right)}}{{2 + \eta - \left( {2 - \eta } \right)\cos \beta }} $$

(5)

Fig. 3

Statics and kinematics of a symmetric rotational collapse mechanism of a circular arch supporting only its own weight

Now, from the rotational equilibrium of total half-arch AC around hinge C at the shoulder extrados, one gets a second equilibrium condition (Fig. 3):

$$ H \cdot \left[ {r + \frac{t}{2} - \left( {r + \frac{t}{2}} \right)\cos \alpha } \right] = W \cdot \left[ {\left( {r + \frac{t}{2}} \right)sin\,\alpha - x_{W} } \right] $$

(6)

where W (W = W₁ + W₂) is the total weight of the half-arch, with half-angle of embrace α, acting at horizontal distance x_W from crown A. Given Eqs. (2)-(3), one has:

$$ W = w\,r\,\alpha \;,\;\;\;\;\;x_{W} = r\;\frac{1 - \cos \alpha }{\alpha } $$

(7)

and, by substituting Eq. (7) into rotational equilibrium Eq. (6) one gets a second equilibrium relation in terms of h and η:

$$ h\,\left( {2 + \eta } \right)\, = \,\left( {2\, + \,\eta } \right)\,\,\underbrace {{\alpha \,\frac{\,sin\,\alpha }{1\, - \,\cos \alpha }}}_{\displaystyle{{A = \alpha \cot \frac{\alpha }{2}}}}\, \, - \,2 = A\left( {2\, + \,\eta } \right) - 2\, $$

(8)

This equilibrium relation is again linear in h (and η), and also linear in group A(α) = αcot α/2, inserting the explicit dependence on half-opening angle α (see plot in Fig. 4). It may linearly be solved with respect to h (or η, or even A = αcot α/2), by obtaining

Fig. 4

Functional dependence of A = α cot(α/2) on α, with indication of stationary (on β) and limit parameters of present CCR solution

$$ h = h_{\,2} = A - \frac{2}{2 + \eta } = \frac{{\left( {2 + \eta } \right)\,\,\alpha \;sin\;\alpha - 2\left( {1 - \cos alpha } \right)}}{{\left( {2 + \eta } \right)\,\,\left( {1 - \cos \alpha } \right)}} $$

(9)

Equilibrium Eqs. (5) and (9) constitute a system of two static equations, which may be condensed in a single one, by eliminating h, as h₁ = h₂, to get the following single limit equilibrium equation:

$$ \frac{{\left( {2 - \eta } \right)\beta \sin \beta - 2\left( {1 - \cos \beta } \right)}}{{2 + \eta - \left( {2 - \eta } \right)\cos \beta }} = A - \frac{2}{2 + \eta } $$

(10)

2.2 Tangency Condition

Now, beside equilibrium, an optimality condition, in the least-thickness condition, shall also be set for the masonry arch, as the tangency condition of the line of thrust at haunch intrados B, according to the textual description in Heyman words [1,2,3,4,5].

However, Heyman [3] seems to have actually analytically stated this tangency condition in terms of the resultant thrust force [6], corresponding to a simplification, in the analysis, leading to a beautiful “linear algebraic problem” [7], as shown below.

Since the angle of inclination of the thrust force in B to the horizontal is such that its tangent is given by ratio W₁/H, which shall then coincide with the local inclination of the inner circle of the intrados profile, from Eq. (2) and relation H = w r h, one has

$$ \frac{{W_{1} }}{H} = \frac{\beta }{h} = tan\,\beta $$

(11)

leading to the following tangency equation for Heyman solution:

$$ {{h}}{ = }{{h}}_{H} = \beta \,\cot \beta $$

(12)

To get instead the true tangency condition of the line of thrust (locus of pressure points) at intrados B, one shall first derive its analytical representation, e.g. in terms of eccentricity e(β)= M/N of the centres of pressure with respect to the centreline of the arch (e taken positive from centreline towards centre O of the circular arch).

Towards that, one first gets a slightly modified and more general version of equilibrium relation (1), by the rotational equilibrium of any upper portion of the half-arch with respect to the centre of pressure at eccentricity e, at a general position β along the half-arch:

$$ H \cdot \left[ {r + \frac{t}{2} - \left( {r - e} \right)\cos \beta } \right] = W_{1} \cdot \left[ {\left( {r - e} \right)sin\,\beta - x_{1} } \right] $$

(13)

By solving this equilibrium relation with respect to eccentricity e (or to non-dimensional eccentricity −1 ≤ ê = 2e/t ≤ 1), again in terms of non-dimensional variables η = t/r and h = H/(wr), one derives the equation expressing the line of thrust as the locus of the centres of pressure of the resultant thrust force:

$$ \hat{e}\text{(}\beta \text{)} = \frac{{2e\text{(}\beta \text{)}}}{t} = \,\,\frac{{2\beta \sin \beta - 2\text{(}1 - \cos \beta \text{)} - h\,\text{(}2 + \eta - 2\cos \beta \text{)}}}{{\eta \;\text{(}\beta \sin \beta + h\cos \beta \text{)}}} = \frac{{num_{{\,\hat{e}}} \text{(}\beta \text{)}}}{{den_{{\,\hat{e}}} \text{(}\beta \text{)}}} $$

(14)

Eccentricity function ê(β) in fractional form in Eq. (14), displaying built-in property ê(0) = −1 at the crown (β = 0), depends on both η and h. By further posing h = h₂, Eq. (9), which would automatically set ê(α) = −1 also at the arch shoulder (β = α), one obtains the final expression of the eccentricity of the line of thrust passing from crown A and shoulder C, at any given value of α (thus of A = α cot α/2) and η:

$$ \hat{e}(\beta ) = \frac{{2e\text{(}\beta )}}{t} = \,\,\frac{{2\,(2 + \eta ) \,\left( {\beta \sin \beta - (1 - \cos \beta )A} \right) - \eta \,\left( {(2 + \eta )A - 2\cos \beta }\right)}}{{\eta \,(2 + \eta )\;\left( {\beta \sin \beta - (1 - \cos \beta )A} \right) + \eta \,\left( {(2 + \eta )A - 2\cos \beta } \right)}} $$

(15)

In the critical, least-thickness condition, the line of thrust touches intrados B where the hinge at the haunch forms:

$$ e\text{(}\beta \,\text{)} = \frac{t}{2}\,\,\,\,\,\,\,or\,\,\,\,\,\,\hat{e}\text{(}\beta \,\text{)} = 1 $$

(16)

and, at the same time, becomes tangent to the intrados of the arch. Thus, function e(β) has to display a stationary point at the haunch section (where e = t/2), i.e. the first-order derivative of eccentricity e(β) with respect to angular inner-hinge position β must (there) vanish:

$$ e^\prime\text{(}\beta \,\text{)} = 0\;\;\;\;\;or\;\;\;\;\;\hat{e}^\prime\text{(}\beta \,\text{)} = 0 $$

(17)

Thus, Eq. (17) has to hold, at the haunch, together with Eq. (16). By expanding, e.g. the derivative of the fractional non-dimensional form of e one has:

$$ \hat{e}^\prime(\beta ) = \frac{{num_{{\hat{e}}}^{\prime }\,den_{{\hat{e}}} - num_{{\hat{e}}} \,den_{{\hat{e}}}^{\prime }}}{{den_{{\hat{e}}}^{2} }} = \frac{{num_{{\hat{e}}}^{\prime } - \hat{e}\,\,den_{{\hat{e}}}^{\prime }}}{{den_{{\hat{e}}} }} $$

(18)

Thus, stationary condition (17) leads to:

$$ \hat{e}\text{(}\beta \,\text{)} = \frac{{num_{{\hat{e}}}^{\prime}(\beta \,\text{)}}}{{den_{{\hat{e}}}^{\prime}(\beta \,\text{)}}} = \frac{2}{\eta }\,\,\frac{{\text{(}\beta \,\sin \beta {)^\prime } - \text{(}1 + h\text{)}\sin \beta }}{{\text{(}\beta \,\sin \beta {)^\prime } - h\sin \beta }} = 1 $$

(19)

since the tangency condition has to hold at haunch B, where ê = 2e/t = 1, and one gets the tangency condition expressed in terms of h as follows:

$$ h = h_{e} = \,\,\frac{{\text{(}2 - \eta \text{)}\,\text{(}\sin \beta + \beta \cos \beta \text{)} - 2\sin \beta }}{{\text{(}2 - \eta \text{)}\sin \beta }} = \overbrace {\beta \,\cot \,\beta }^{{h_{H} }} - \frac{\eta }{2 - \eta } $$

(20)

where term f = (β sin β)′ = sin β + β cos β is involved (see Fig. 5 and later discussion).

Fig. 5

CCR solution. Functional dependence of functions f, g, f ² − 2 g S + S ² on β

Fig. 6

Plot with non-dimensional thrusts h_U, h_L, h₂ of the full semi-circular arch (α = π/2) according to Heyman solution: the three lines slightly intersect

This shows that assuming h_e = h_H, as taken by Heyman, leads to an approximation of true tangency condition h = h_e. This might look reasonable until η keeps small. All this makes the correct solution slightly more involved than the former, and leading to shift from a “linear” to a “quadratic algebraic problem” [7].

Finally, the governing system of the least-thickness masonry arch optimization problem for a self-weight distribution along geometrical centreline may be stated, in terms of h, as follows:

$$ \,\left\{ \begin{aligned} &h = h_{1} = \frac{{\left( {2 - \eta } \right)\,\beta \sin \beta - 2\left( {1 - \cos \beta } \right)}}{{2 + \eta - \left( {2 - \eta } \right)\cos \beta }}\,\,\, \hfill \\ & h = h_{2} = A - \frac{2}{2 + \eta }\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad{\text{Heyman}}/CCR \hfill \\ &h = h_{e} = \frac{{\text{(}2\, - \,\eta \text{)}\,\text{(}\sin \beta + \beta \cos \beta \text{)}\, - \,2\sin \beta }}{{\text{(}2 - \eta \text{)}\sin \beta }}\,\, = \,\,\overbrace {\beta \cot \beta \,}^{{h_{H} }} - \frac{\eta }{2\, - \,\eta }\,\,\delta_{CCR} \,\,\,\,\,\,\, \hfill \\ \end{aligned} \right. $$

(21)

where A = α cot α/2 (Fig. 4) and δ_CCR is an on/off control flag allowing to shift from Heyman (δ_CCR = 0) to CCR solution (δ_CCR = 1).

The first two equations, eliminating h, as h₁ = h₂, set the equilibrium relation, Eq. (10), the third equation, i.e. Eq. (20), sets the tangency (optimality) condition, in a shift between Heyman and CCR solutions. These equations will be the subject of a further separate analytical treatment in the following section.

3 Alternative Analytical Derivation and Interpretation

Now, further, independent and reconciling ways to derive the least-thickness condition, leading to same results, are here presented, namely: a “Coulomb’s static approach”, based on the so-called upper and lower horizontal thrusts and a “Mascheroni’s kinematic approach”, based on the balance of virtual work (or power).

3.1 Coulomb’s Static Approach

An alternative way of deriving the least-thickness condition could be based on a “Coulomb’s static approach”, according to the terminology adopted in Sinopoli et al. [22]. Specifically, reference is here made to the original derivations presented in Blasi and Foraboschi [19], which account for the explicit determination of the “upper” and “lower” horizontal thrusts. In practice, according to Coulomb’s view, it is stated that, to warrant equilibrium, the horizontal thrust should take values that are in between the minimum and maximum values that the horizontal thrust could assume. Such two limit values can be obtained as described below.

Lower horizontal thrust H_min(β) = H_L = H₁ is the minimum value of the horizontal thrust, applied to the extrados at crown A, which corresponds to impose the rotational equilibrium of any upper portion of the half-arch, of variable half-opening β, with respect to intrados inner-hinge B. In practice, this coincides with earlier-mentioned value H₁ in Sect. 2. Upper horizontal thrust H_max(β) = H_U is the maximum value of the horizontal thrust, applied to the extrados at shoulder C, which corresponds to impose the rotational equilibrium of any lower portion of the half-arch, of variable (α−β) opening, with respect to intrados inner-hinge B.

With reference to Fig. 3, and following the same type of earlier-explained derivations, the lower and upper horizontal thrusts may be determined by the above-described equations of rotational equilibrium as follows:

$$ H_{L} = H_{1} = \frac{{W_{1} \, \cdot \,\left[ {\left( {r - \dfrac{t}{2}} \right)sin\,\beta - x_{1} } \right]}}{{\left( {r + \dfrac{t}{2}} \right) - \left( {r - \dfrac{t}{2}} \right)\cos \beta }}\,\, $$

(22)

$$ H_{U} = \frac{{W\, \cdot \,\left[ {\left( {r + \frac{t}{2}} \right)sen\,\alpha - \left( {r - \frac{t}{2}} \right)sin\,\beta } \right]\,\, - \,\,W_{2} \cdot x_{2B} }}{{\left( {r - \frac{t}{2}} \right)\cos \beta - \left( {r + \frac{t}{2}} \right)\cos \alpha }} $$

(23)

where:

$$ \begin{aligned} W_{1} = wr\beta \,,\,\,\,\,W_{2} = wr\,\text{(}\alpha - \beta \text{)}\,,\,\,\,W = W_{1} + W_{2} = wr\,\alpha \,\text{;}\,\text{ } \hfill \\ x_{1} = r\,\,\frac{1 - \cos \beta }{\beta },\,\,\,\,x_{2} \, = r\,\,\frac{\cos \beta - \cos \alpha }{\alpha - \beta },\,\,\,\,x_{2B} \, = x_{2} - \left( {r - \frac{t}{2}} \right)\sin \beta \hfill \\ \end{aligned} $$

(24)

Thus, in usual non-dimensional form (η = t/r, h = H/(wr)), one has:

$$ h_{L} = h_{1} = \frac{{\left( {2 - \eta } \right)\,\,\beta \,sin\,\beta - 2\,\,\left( {1 - \cos \beta } \right)}}{{\left( {2 + \eta } \right) - \left( {2 - \eta } \right)\,\,\cos \beta }} = \frac{{num_{{\,h_{L} }} }}{{den_{{\,h_{L} }} }} $$

(25)

$$\begin{aligned} h_{\,U\,} &= \,\frac{{\left( {2 + \eta } \right)\,\alpha \,sin\,\alpha - \left( {2 - \eta } \right)\,\beta \,sin\,\beta - 2\left( {\cos \beta - \cos \alpha } \right)}}{{\left( {2 - \eta } \right)\,\,\cos \beta - \left( {2 + \eta } \right)\,\,\cos \alpha }}\,\\ & = \,\frac{{num_{{\,h_{\,U} }} }}{{den_{{\,h_{\,U} }} }}\, = \,\frac{{num_{{\,h_{2} }} - num_{{\,h_{L} }} }}{{den_{{\,h_{2} }} - den_{{\,h_{L} }} }}\end{aligned} $$

(26)

where the quantities at the numerators and denominators in the fractional forms have been introduced, and use has been made of Eq. (9), stating the rotational equilibrium of the total half-arch with respect to shoulder hinge C, also re-written as well below in fractional form:

$$ h_{\,2} = A - \frac{2}{2 + \eta } = \frac{{\left( {2 + \eta } \right)\,\,\alpha \;sin\;\alpha - 2\left( {1 - \cos \alpha } \right)}}{{\left( {2 + \eta } \right)\,\,\left( {1 - \cos \alpha } \right)}} = \frac{{num_{{\,h_{2} }} }}{{den_{{\,h_{2} }} }} $$

(27)

Equation (25) confirms that h_L corresponds to same value h₁ earlier derived in Eq. (5); thus h_L = h₁ and, through h₁, previously, use was already made of the lower thrust. Equation (26) represents an additional expression of the non-dimensional horizontal thrust, which is alternative to that of h₂, Eq. (27), earlier accounted for. This could be used to write the second ruling equilibrium equation in an alternative way.

Notice that, instead of what is happening for h₂, the dependence on α only through group A is not apparent in Eq. (26). However, it is straightforward to show that, based on the last relation that has been written in Eq. (26), stating the two equilibrium equations as h = h_L and h = h_U, i.e. stating equilibrium as h_U = h_L, is equivalent to state it as h_L = h₂, as previously done. Similarly, stating equilibrium as h_U = h₂ is also equivalent to state h_L = h₂. Thus, h_U could be used instead of h₂, which is dependent only on A, with same results. Finally, stating equilibrium as h_L = h₂, like previously done, is also equivalent to state h_U = h₂.

In practice, out of three equilibrium-linked horizontal thrusts h_L = h₁, h₂, h_U, any of the three equivalent equilibrium conditions below, based on two of them, could be employed to state the equilibrium equation:

$$ h_{L} = h_{2} \,\,\,\,\, \Leftrightarrow \,\,\,\,\,h_{L} = h_{U} \,\,\,\,\, \Leftrightarrow \,\,\,\,\,h_{2} = h_{U} $$

(28)

Notice that, here, all three thrusts h_L, h₂, h_U depend on η. Thrust h_L = h₁ depends only on angular position β; h₂ only on angular position α; h_U on both angular positions α and β.

As debated by Blasi and Foraboschi [19], see Figures 7 and 8 in their paper, h_L(β) and h_U(β), as a function of β, respectively, provide lower and upper bounds for h. Since such two curves depend on η, at a given value of α in h_U, they may turn out as follows: detached, providing, with a positive minimum relative distance (clearance), a measure of the “margin of safety” for the arch; intersecting, with a negative minimum relative distance, at same quote h_L = h_U = h₂, denouncing a sub-critical condition; tangent to each other at a zero relative distance, at a point where h_L = h_U = h₂, locating the true critical least-thickness condition.

Fig. 7

Plot with the non-dimensional thrusts h_U, h_L, h₂ of the full semi-circular arch (α = π/2) according to present CCR solution: the three lines are truly tangent

Fig. 8

Plot with the non-dimensional thrusts h_U, h_L, h₂ of the arch with α = 7 π/9 according to Heyman solution: the three lines clearly intersect

Also, since three curves h_U(β), h_L(β), h₂, the latter corresponding to a horizontal line at constant quote h₂, intersect, if it happens, at same quote h₂, disregarding the values of β corresponding to h_U = h_L, the intersection between h_U and h_L occurs at constant h, i.e. with Δh = 0. Thus, in the limit of the tangency condition at the minimum thickness, curves h_U(β) and h_L(β) are both stationary at the point where h_L = h_U, i.e. their local tangent is horizontal, as that of constant h₂. This states the critical condition as zero first-order derivative with respect to β, h_L′ = 0 or h_U′ = 0, where equilibrium Eqs. (28) hold. At such stationary points of h_U(β) and h_L(β), one has h_min = h_L = h₁ and h_max = h_U, with positive clearance h_max − h_min > 0 in over-safe condition h_max > h_min, negative clearance h_max − h_min < 0 in sub-critical condition h_max < h_min and no clearance h_max − h_min = 0 in critical condition h_max = h_min.

Now, it is quite straightforward to show, from the fractional forms reported in Eqs. (25)–(26), that the stationary conditions on either h_L or h_U, or even the condition of mutual tangency h_L′ = h_U′, as adopted by Blasi and Foraboschi [19], at h_L = h_U, lead to the same tangency condition as analysed earlier for the line of thrust, which was stated in terms of h = h_e, Eq. (20), through the definition of the line of thrust. Indeed, at the stationary points of h_L and h_U:

$$ {h_{L,U}}^{\prime } = \frac{{{num_{{\,h_{L,U} }}}^{\prime } - h_{L,U} \,{den_{{\,h_{L,U} }}}^{\prime } }}{{den_{{\,h_{L,U} }} }} = 0\;\,\,\,\;\, \Rightarrow \;\,\,\,\,h_{L,U} = \frac{{{num_{{\,h_{L,U} }}}^{\prime } \,}}{{{den_{{\,h_{L,U} }}}^{\prime } }}\, $$

(29)

where

$$ {\begin{aligned}num_{{\,h_{e} }} &= {num_{{\,h_{L} }}}^{\prime } = - \,{num_{{\,h_{U} }}}^{\prime } = \left( {2 - \eta } \right)\,\left( {\beta \,sin\,\beta } \right)^{\prime } - 2\,\sin \beta = \hfill \\ &= \left( {2 - \eta } \right)\,\left( {sin\,\beta \text{ + }\beta \,\cos \beta \,} \right) - 2\,\sin \beta = \left( {2 - \eta } \right)\,\beta \,\cos \beta \, - \eta \,\sin \beta \hfill \\ den_{{\,h_{e} }} &= {den_{{\,h_{L} }}}^{\prime } = - \,{den_{{\,h_{U} }}}^{\prime } = \left( {2 - \eta } \right)\,\sin \beta \hfill \\ \end{aligned}} $$

(30)

Also, alternative condition h_L′ = h_U′, at h_L = h_U, from (29) to (30), leads to:

$$ \begin{aligned}h_{L}^{\prime } &= h_{U}^{\prime } \;\,at\,\,\;h_{L}= \,h_{U} \\ & \,\,\,\,\, \Rightarrow \;\,\,\,\,h_{e} = \frac{{num_{{\,h_{e} }} \,}}{{den_{{\,h_{e} }} }} = \frac{{{num_{{\,h_{L} }}}^{\prime } \,}}{{{den_{{\,h_{L} }}}^{\prime } }} = \frac{{den_{{\,h_{U} }} h_{L} + \,den_{{\,h_{L} }} h_{U} }}{{den_{{\,h_{U} }} + den_{{\,h_{L} }} }} = h_{L} = h_{U} \,\end{aligned} $$

(31)

Thus, all the tangency conditions above are equivalent to h = h_e, as earlier directly stated on the line of thrust:

$$ \begin{aligned} &h_{L}^{\prime } = \,0\,\,\,\, \Leftrightarrow \,\,\,\,h_{U}^{\prime } = \,0\,\,\,\, \Leftrightarrow \,\,\,h_{L}^{\prime } = \,h_{U}^{\prime } \\ & \Leftrightarrow \,\,\,\,h\, = \,h_{e} \, = \,\,\frac{{\text{(}2\, - \,\eta \text{) (}\beta \sin \beta \text{)}^{\prime } - \,2\sin \beta }}{{\text{(}2\, - \,\eta \text{) }\sin \beta }}\, = \beta \cot \beta \, - \,\frac{\eta }{2\, - \,\eta } \hfill \\ \end{aligned} $$

(32)

This actually signals the easiest way to account for the tangency condition of the line of thrust at the intrados hinge. Once determined from equilibrium, the simpler expression of lower thrust h_L = h₁ (dependent just on β), its stationary condition h_L′ = 0 at h = h_L, immediately leads to tangency condition h = h_e, without even the need of passing through the definition of the line of thrust and of its e(β) eccentricity (Sect. 2.2).

In conclusion, given four thrust equations h = h_L = h₁, h = h_U, h = h₂, h = h_e, any of the four possible systems with three of them would equivalently lead to a correct least-thickness solution representation:

$$ \left\{ {\begin{array}{*{20}l} {h = h_{L} = h_{1} } \hfill \\ {h = h_{2} } \hfill \\ {h = h_{e} } \hfill \\ \end{array} } \right. \Leftrightarrow \left\{ {\begin{array}{*{20}l} {h = h_{L} = h_{1} } \hfill \\ {h = h_{U} } \hfill \\ {h = h_{e} } \hfill \\ \end{array} } \right. \Leftrightarrow \left\{ {\begin{array}{*{20}l} {h = h_{U} } \hfill \\ {h = h_{2} } \hfill \\ {h = h_{e} } \hfill \\ \end{array} } \right. \Leftrightarrow \left\{ {\begin{array}{*{20}l} {h = h_{L} = h_{1} } \hfill \\ {h = h_{U} } \hfill \\ {h = h_{2} } \hfill \\ \end{array} } \right. $$

(33)

These systems could be taken for a numerical solution of critical parameters β, η, h at given half-angle of embrace α. The solutions of the first three systems, involving tangency condition h = h_e, are expected to be more efficient than that of the last system, based just on equating three thrusts h_L = h₁, h_U, h₂. However, this latter system based just on three equilibrium relations, despite not explicitly accounting for the tangency condition, could be used as well for the final solution. The first system is still probably the simplest, also conceptually. As stated above, this corresponds to set the equilibrium relation as equilibrium condition h_L = h₂ and the stationary condition as either tangency condition e′ = 0 or h_L′ = 0.

All the above-outlined considerations can be inspected in the plots of horizontal thrusts h_L = h₁, h_U, h₂ depicted in Figs. 6, 7, 8, and 9, which, respectively, show the thrust functions depending on β, for the two reference cases of α = 90° and α = 140°, in comparison for Heyman and CCR solutions. As commented in [7], for the corresponding plots of the eccentricity of the line of thrust, this clearly shows Heyman solution to turn out sub-critical. Indeed, the three lines of h_L, h_U, h₂ cross within a fork of two values of β that surround the true value of β determined by CCR solution (which is almost in the middle), correctly leading, instead, to a true tangency condition. In Figs. 6, 7, 8, and 9, curve h_e(β) is as well reported, which is useful to represent the stationary, thus tangency, condition, by locating the stationary points of curves h_U(β) and h_L(β). Also, the plots of CCR solution show that the direct use of tangency condition h = h_e shall numerically become more effective, since the cutting over h_U(β), h_L(β), h₂ is sharper than that going through the quite flat stationary points at the solution of system h = h_L, h = h_U, h = h₂.

Fig. 9

Plot with the non-dimensional thrusts h_U, h_L, h₂ of the arch with α = 7 π/9 according to present CCR solution: the three lines are truly tangent

Also, it is confirmed that curves h_L(β) ed h_U(β) turn out quite flat near the common stationary point in the critical condition, which locates correct angular position β of the inner-hinge. Thus, as noted by Heyman in quoting Coulomb’s observations, it is clear that even approximate estimations of angular inner-hinge position β might lead to fairly correct values of h (and also of β). For this reason, despite missing the correct estimate of β, as instead obtained by CCR solution, Heyman solution still appears quite acceptable in engineering terms (at least for under-complete arches).

3.2 Mascheroni’s Kinematic Approach

An additional, independent way of deriving the least-thickness condition could be based on a “Mascheroni’s kinematic approach”, following again the terminology adopted in Sinopoli et al. [22]. Indeed, the below reported solution is derived from an alternative kinematic approach, based on the principle of virtual work (or power), which is written with reference to the purely rotational rigid-body five-hinge collapse mechanism of the arch (see Figs. 2 and 3).

Referring to the potential three-hinge rigid-body kinematic chain of the half-arch in Fig. 3, one takes the external virtual work (or power) equation, to state equilibrium at incipient collapse:

$$ \dot{L}_{e} = 0 $$

(34)

i.e., explicitly, for any nonzero angular rotations ψ and φ (or velocities $\dot{\psi }$ and $\dot{\varphi }$):

$$ \dot{L}_{e} = - \,\,W_{1} \,\, \cdot \,\,\,\left( {x_{1} - x_{{\varOmega_{1} }} } \right)\,\,\dot{\psi }\; - \,\,W_{2} \,\, \cdot \,\,\,\left[ {\left( {r + \frac{t}{2}} \right)\sin \alpha - x_{2} } \right]\,\,\dot{\varphi } = 0 $$

(35)

where, given Eqs. (2)–(3) and (24), W₁ = wrβ and W₂ = wr (α − β, are the weights of the two portions of the half-arch separated by inner haunch hinge B (thus, with total weight of the half-arch W = W₁ + W₂ = wrα), acting at abscissas x₁ = r (1 – cos β)/β and x₂ = r (cos β – cos α)(α − β) = x_2B+(r − t/2) sin β from the vertical axis of symmetry at the crown; x_Ω1 is the horizontal distance of the centre of rotation Ω₁ of the upper portion of the half-arch with respect to crown A and x_2B is the horizontal distance from the line of action of W₂ to point B. Equation (34), i.e. explicitly (35), is evidently a way to state equilibrium at the virtual rotational collapse mechanism that may develop.

The kinematic chain is a one-degree-of-freedom system. The relation between two angular velocities $\dot{\psi }$ and $\dot{\varphi }$ in Fig. 3 can be obtained by imposing that the horizontal velocity of inner-hinge B is the same for the two portions of the arch. Thus, one obtains the following kinematic link:

$$ \left[ {\left( {r - \frac{t}{2}} \right)\cos \beta - \left( {r + \frac{t}{2}} \right)\cos \alpha } \right]\;\dot{\varphi }\;\; = \;\;\left[ {\left( {r + \frac{t}{2}} \right) - \left( {r - \frac{t}{2}} \right)\cos \beta } \right]\dot{\psi }\; $$

(36)

Moreover, abscissa x_Ω1 can as well be determined, by imposing that the vertical velocity of inner-hinge B is the same for the two portions of the arch, by obtaining:

$$ \;\left[ {\left( {r + \frac{t}{2}} \right)\sin \alpha - \left( {r - \frac{t}{2}} \right)\sin \beta } \right]\;\dot{\varphi }\;\; = \;\;\left[ {\left( {r - \frac{t}{2}} \right)\sin \beta - x_{{\varOmega_{1} }} } \right]\,\,\dot{\psi }\; $$

(37)

namely:

$$ x_{{\varOmega_{1} }} \;\dot{\psi } = \; - \;\left[ {\left( {r + \frac{t}{2}} \right)\sin \alpha - \left( {r - \frac{t}{2}} \right)\sin \beta } \right]\,\,\dot{\varphi }\;\; + \;\;\left[ {\left( {r - \frac{t}{2}} \right)\sin \beta } \right]\,\,\dot{\psi }\; $$

(38)

By substituting Eq. (38) into virtual work (or power) Eq. (35) and by making use of relation (36) and explicitly of W = W₁ + W₂, x₂ = x_2B+(r − t/2) sin β, one has, in view of relations (22)-(23), for lower and upper thrusts H_L and H_U:

$$ \frac{{\dot{L}_{e} }}{{\;\left[ {\left( {r + \frac{t}{2}} \right) - \left( {r - \frac{t}{2}} \right)\cos \beta } \right]\,\,\dot{\psi }}} = \;\frac{{\dot{L}_{e} }}{{\left[ {\left( {r - \frac{t}{2}} \right)\cos \beta - \left( {r + \frac{t}{2}} \right)\cos \alpha } \right]\;\dot{\varphi }}} = \;H_{L} - H_{U} \; = 0 $$

(39)

Thus, stating $\dot{L}$_e = 0 as equilibrium condition is fully equivalent to state it as H_L = H_U, namely the same equilibrium condition found from previous Coulomb’s static approach. In other words, and in non-dimensional terms, $\dot{L}$_e = 0 is equivalent to assume an equilibrium condition in the form h_L = h_U:

$$ \dot{L}_{e} = 0\,\,\,\,\, \Leftrightarrow \,\,\,\,h_{L} = h_{U} $$

(40)

A second relation, which indirectly expresses the tangency condition of the line of thrust at intrados B, can be obtained by setting to zero the derivative with respect to β of the external virtual work (or power):

$$ \dot{L}_{e}^{\prime } = \frac{{\partial \,\dot{L}_{e} }}{\partial \beta } = 0\, $$

(41)

This comes from the following consideration: if the line of thrust has to turn out tangent to the intrados at haunch B, it should also not that much deviate from the circular intrados curve in the surroundings of B. As a consequence, equilibrium should be warranted also for small variations of point B, thus of angular position β, and this holds true if also the angular derivative of the external virtual work (or power), at intrados B, vanishes.

Given obtained Eq. (39), alternative stationary condition (41), becomes equivalent to:

$$\begin{aligned} \dot{L}_{e}^{\prime }& = \left( {H_{L}^{\prime } \, - H_{U}^{\prime } } \right)\, \cdot \left[ {\left( {r + \frac{t}{2}} \right) - \left( {r - \frac{t}{2}} \right)\cos \beta } \right]\,\dot{\psi }\\ &\quad + \left( {H_{L} \, - H_{U} } \right)\, \cdot \,\left( {r - \frac{t}{2}} \right)\sin \beta \,\,\dot{\psi } = 0 \end{aligned}$$

(42)

Thus, at $\dot{L}$_e = 0, i.e. at H_L = H_U, $\dot{L}$_e′ = 0 is equivalent to H_L′ = H_U′, i.e. one of the equivalent ways to set the tangency condition according to earlier Coulomb’s static approach. In other words, and in non-dimensional terms:

$$ \dot{L}_{e}^{\prime } = 0\,\,\,\,\, \Leftrightarrow \,\,\,\,\,h_{L}^{\prime } = h_{U}^{\prime } $$

(43)

so that system {$\dot{L}$_e = 0, $\dot{L}$_e′ = 0} formed by Eqs. (34) and (41) is wholly equivalent to system {h_L = h_U, h_L′ = h_U′} earlier obtained by Coulomb’s static approach. Thus, the two approaches are fully equivalent, and both lead to the same solution (as earlier derived), as it could be checked by independent numerical evaluations of the various solving systems that have been here derived.

3.3 Discussion on Heyman and CCR Solutions

Classical Heyman solution, directly solved from governing system (21), for δ_CCR = 0, finally leads to [3]:

$$ \,\,\,\left\{ \begin{aligned}& A_{H} = \beta \cot \beta \frac{{2\beta \cos \beta + \sin \beta \cos^{2} \beta + \sin \beta }}{{2\beta \cos \beta + \sin \beta \cos^{2} \beta - \sin \beta \cos \beta }} \hfill \\ & \hfill \\ &\eta_{H} = 2\frac{{\,\left( {\beta \; - sin\,\beta } \right)\left( {1 - \cos \beta } \right)}}{{\beta \;\left( {1 + \cos \beta } \right)}}\,\,\,\,\,\,\,\,\,\, \hfill \\ &\hfill \\ & h_{H} \, = \beta \cot \beta ,\,\,\,\,\,\,\,\,\,\,\,\hat{h}_{H} = \eta_{H} \,h_{H} = \,\,2\cot \beta \frac{{\,\left( {\beta \; - sin\,\beta } \right)\left( {1 - \cos \beta } \right)}}{{\left( {1 + \cos \beta } \right)}} \hfill \\ \end{aligned} \right.\;\;\;\;\;\; $$

(44)

corresponding to a “linear algebraic problem”, with achieved monotonic trends of solution triplet A(β), η(β), h(β) as a function of β.

Instead, for CCR solution (δ_CCR = 1), by eliminating in turn couples η and h, A and h, A and η, from original system (21), three quadratic equations in A(β), η(β), h(β) can be obtained. All this shows that correct CCR solution turns out “quadratic” in solution parameters A, η, h (or $\hat{h}$), as opposed to “linear” Heyman solution. As a consequence of that, parameters A(β), η(β), h(β) (or $\hat{h}$(β)) become double-valued functions of inner-hinge position β. The following “quadratic algebraic problem” in A, η and h (or $\hat{h}$) is then obtained [7]:

$$ \left\{ \begin{aligned}& S\text{(}2g - S\text{)}\,\,A^{2} - 2fg\,\,A + g^{2} = 0\,\,\text{;}\, \hfill \\ & \text{(}f + g\text{)}\,\,\,\eta^{2} - 4\,\text{(}g - S\text{)}\,\,\eta \; + \;4\,\text{(}g - f\text{)} = 0\,\,\text{;}\,\,\,\,\,\,\,\, \hfill \\ & 2S\,\,h^{2} - 2\text{(}f - S\text{)}\,\,h + g - f = 0,\\ &\text{(}f + g\text{)}\,S\,\,\hat{h}^{2} - 2\text{(}g - f\text{)}\;\text{(}f + S\text{)}\,\,\hat{h}\,\, + \,2\text{(}g - f\text{)}^{2} \, = 0\, \hfill \\ \end{aligned} \right. $$

(45)

where functions f(β), g(β) have been defined as follows:

$$ {\begin{aligned} \,\, & \hfill \,\,f = f\text{(}\beta \,\text{)} = \text{(}\beta \,\sin \beta \text{)}^{\prime = }\sin \beta + \beta \cos \beta \text{ = }S + \beta \,C\,\,, \hfill \\ &\,\,g = g\text{(}\beta \,\text{) = }f\text{(}\beta \,\text{) }\cos \beta \, + \beta \sin^{2} \beta = \beta + \sin \beta \cos \beta \,\text{ = }\beta + S\,C\,\, \hfill \\ & \,\,where\,\,\,\,S = S\text{(}\beta \,\text{)} = \sin \beta \,,\,\,\,C = C\text{(}\beta \,\text{)} = \cos \beta , \,\,with\,\,link\,\,\beta S^{2} + Cf - g = 0\,\, \hfill \\ \hfill \\ \end{aligned}} $$

(46)

As noticed, function f spontaneously comes from the derivative term inherited in the expression of h_e, Eq. (20). Functions f and g, with link βS ²+ C f – g =0, turn out to be useful to express the solution in compact form. Notice that g ≥ 0 for β ≥ 0; also, g – f =(1 – cos β)(β – sin β) ≥ 0 and f + g =(1+ cos β)(β+ sin β) ≥ 0 for β ≥ 0. Functions f(β) and g(β) both vanish at β = 0. They are jointly represented in Fig. 5.

Different ways to rewrite the first of Eqs. (45) similarly to Heyman formula (44)_a, with term A isolated on the right-hand side, could be the following:

$$ \begin{aligned} \frac{{2fg - g^{\text{2}} /A}}{{S\,\text{(}2g - S\text{)}}} = A\,\,,\,\,\,\,\,\frac{2fg}{{S\,\text{(}2g - S\text{) + }g^{\text{2}} /A^{\text{2}} }} = A\, \hfill \\ \frac{{g^{2} }}{{2fg - S\,\text{(}2g - S\text{)}A}} = A\,\,,\,\,\,\,\,\frac{{g^{2} + S\,\text{(}2g - S\text{)}A^{2} }}{2fg} = A \hfill \\ \end{aligned} $$

(47)

though obviously there now appears a dependence on A, thus on α, also on the left-hand side of the equation (thus inspiring possible recursive evaluations, as treated in the next section).

Similarly, for the sake of completeness, the following expressions could as well be written for η:

$$ \begin{aligned} \frac{{4\,\text{(}g - S\text{)}\; - 4\text{(}g - f\text{)}\,/\eta }}{f + g} = \eta \,\,,\,\,\,\,\,\frac{{4\,\text{(}g - S\text{)}}}{{\,f + g\text{ + }4\,\text{(}g - f\text{)}\,/\eta^{\text{2}} }} = \eta \, \hfill \\ \frac{{\,4\,\text{(}g - f\text{)}}}{{4\text{(}g - S\text{)} - \text{(}f + g\text{)}\,\eta }} = \eta \,\,,\,\,\,\,\frac{{4\,\text{(}g - f\text{)}\,\; + \text{(}f + g\text{) }\eta^{2} }}{{\text{ }4\;\text{(}g - S\text{)}}} = \eta \hfill \\ \end{aligned} $$

(48)

and for h:

$$ \begin{aligned} \frac{{2\,\text{(}f - S\text{)}\, - \,\text{(}g - f\text{)}\,/h}}{2S\,} = h\,\,,\,\,\,\,\,\frac{{2\,\text{(}f - S\text{)}}}{{2S\,\text{ + }\,\text{(}g - f\text{)}\,/h^{\text{2}} }} = h\, \hfill \\ \frac{g - f}{{2\,\text{(}f - S\text{)}\, - 2S\,h}} = h\,\,,\,\,\,\,\,\frac{{g - f\, + 2S\,h^{2} }}{{2\,\text{(}f - S\text{)}}} = h \hfill \\ \end{aligned} $$

(49)

or for ĥ:

$$ \begin{aligned} 2\,\text{(}g - f\text{)}\,\frac{{f + S\, - \,\text{(}g - f\text{)}\,/\hat{h}}}{{\text{(}f + g\text{)}\,S\,}} = \hat{h}\,\,,\,\,\,\,\,\frac{{2\,\text{(}g - f\text{)}\,\text{(}f + S\text{)}}}{{\text{(}f + g\text{)}\,S\,\text{ + }\,2\,\text{(}g - f\text{)}^{2} /\hat{h}^{2} }} = \hat{h}\, \hfill \\ \frac{{2\,\text{(}g - f\text{)}^{2} }}{{2\,\text{(}g - f\text{)}\,\text{(}f + S\text{)}\, - \,\text{(}f + g\text{)}\,S\,\hat{h}}} = \hat{h}\,\,,\,\,\,\,\,\frac{{2\,\text{(}g - f\text{)}^{2} \, + \text{(}f + g\text{)}\,\,S\,\hat{h}^{2} }}{{2\,\text{(}g - f\text{)}\,\text{(}f + S\text{)}}} = \hat{h} \hfill \\ \end{aligned} $$

(50)

The solution of CCR system (21), for δ_CCR = 1, or of quadratic system (45) can then be obtained in compact form, by explicitly solving for triplet (A, η, h):

$$ \,\left\{ \begin{aligned} & A_{CCR} = g\,\,\,\frac{{f \pm \sqrt {f^{2} - 2g\,S + S^{2} } }}{{S\text{(}2g - S\text{)}}}\,\,\text{;} \hfill \\ &\eta_{CCR} \, = 2\,\,\,\frac{{g - S \mp \sqrt {f^{2} - 2g\,S + S^{2} } }}{f + g}\,\,\text{;}\, \hfill \\ & h_{CCR} = \,\,\frac{{f - S \pm \sqrt {f^{2} - 2g\,S + S^{2} } }}{2\,S}\,\,\text{,} \\ &\hat{h}_{CCR} = \eta \,h = \text{(}g - f\text{)}\frac{{f + S \pm \sqrt {f^{2} - 2g\,S + S^{2} } }}{{S\text{(}f + g\text{)}}}\, \hfill \\ \end{aligned} \right. $$

(51)

where the signs in front of the square root terms are sorted out in triplet (+,–, +) or in triplet (–, + ,–), i.e. with sorting (A⁺,η⁻, h⁺) and (A⁻,η⁺, h⁻). Non-dimensional horizontal thrust ĥ can then be determined as well by the product of η and h (and sorted out in the same way as that of h). Notice that A(β), thus α(β), η(β) and h(β), or $\hat{h}$(β), are two-valued functions of β, i.e. there are two values of A, thus of α, η and h, or $\hat{h}$, that correspond to same inner-hinge position β. Term f² − 2 g S + S² is ≥ 0 for 0 ≤ β ≤ β^CCR_sβ (Fig. 5), where β^CCR_sβ is the root of f² − 2 g S + S² = 0, assuring that the solution turns out real-valued in that range of β.

A direct graphical confrontation between Heyman and CCR solutions can be appreciated in Fig. 10, where triplet A, η, h is represented by analytical plots as a function of angular inner-hinge position β. This representation in terms of β allows to highlighting the differences between the two solutions. This is mainly due to the dissimilar trends experienced on the estimated hinge position. Indeed, the trends of A, η, h are monotonic (single-valued) for Heyman solution and non-monotonic (double-valued) for CCR solution, with a very appreciable deviation in terms of β, especially at increasing half-opening angle α (decreasing A), already when approaching complete semicircular arch case α = A = π/2 and more and more when α goes beyond that (thus the greatest differences are revealed for over-complete arches with half-angles of embrace larger than 90°). Despite that, since as stated by Heyman, the hinge position is somehow an internal ingredient in the solution, in engineering terms, the final differences on η and h at variable angle of embrace are rather limited [7]. In the plots in Fig. 10, the trends for β small, corresponding to small angles of embrace α, which turn out the same for Heyman and CCR solutions, are as well represented.

Fig. 10

Heyman vs. CCR solution. Functional dependence of A, η, h on β

A further detailed representation of CCR solution is provided in Fig. 11, where triplet β, η, h is depicted as a function of A, by analytical parametric plots, with indication of typical characteristic values of the various involved parameters, with a nomenclature introduced in [7].

Fig. 11

CCR solution. Functional dependence of β, η, h on A = α cot(α/2)

The true appearance of the line of thrust that develops within the arch in the critical least-thickness condition is also analytically represented for CCR solution by the polar plots reported in Figs. 12, 13, 14 and 15, for some characteristic values of the half-angle of embrace, including: a reference case for α < 90°, i.e. α = 70°; the characteristic case of α^CCR_sβ that corresponds to the stationary condition of curve β^CCR(α) [7, 14]; another taken over-complete reference case, i.e. α = 140°; limit case α = α^CCR_l of CCR solution, leading to a vanishing horizontal thrust [7]. The plots make apparent the increasing thickness that the arch shall display to warrant self-standing equilibrium, at increasing opening angle, with a corresponding decrease in non-dimensional horizontal thrust h (and actually a bell-shaped trend of intrinsic non-dimensional horizontal thrust ĥ, going through a maximum [7, 14]). The lines of thrust start to “bend” for an α that is around that leading to the stationary condition on β and reach, in the limit configuration, a theoretical profile that gets to the intrados on the vertical axis of symmetry at the crown. This corresponds to a precarious, inverted pendulum equilibrium configuration, that is achieved by a resultant self-weight of the half-arch that is exactly vertically aligned on the underlying bearing hinge at the shoulder, with a zero transverse horizontal thrust (h = 0) and inner-hinge that disappears, by pulling-back to zero (β = 0), with a released section at the crown (giving rise to an overturning mechanism, at infinite friction [13]). In such a limit case, according to CCR solution, the thickness of the arch becomes equal to its radius (η = 1).

Fig. 12

Analytical plot of the line of thrust for CCR solution at α = 7π/18 = 70° (A = 1.74481)

Fig. 13

Analytical   plot   of   the   line   of   thrust   for   CCR   solution  at  α = α_sβ = 2.23031 = 127.788° (A = A _sβ = 1.09292). Widest angular inner-hinge position from the crown

Fig. 14

Analytical plot of the line of thrust for CCR solution at α = 7π/9 = 140° (A = 0.889347)

Fig. 15

Analytical   plot   of   the   line   of   thrust   for   CCR   solution   at   α = A_l = 2.64839 = 151.742° (A = A_l = 2/3). Limit case of zero horizontal thrust

3.4 Generalization to Milankovitch Solution

So far, the assumption of Heyman uniform self-weight distribution along the geometrical centreline of the arch was considered. Instead, due to the curvature of the circular arch and to the resulting wedged-shape of each ideal infinitesimal chunk of the continuous arch, it appears that its centre of gravity is slightly radially displaced, at radial distance r_G a bit larger than radius r:

$$ \begin{aligned} r_{G} \, &= \,\frac{{\int_{\,dV} {\gamma \rho \,dV} }}{{\int_{\,dV} {\gamma dV} }}\, = \,\frac{{\int_{\,0}^{\,d\varphi } {} \int_{\,r\, - \,t/2}^{\,r\, + \,t/2} {\gamma \rho \overbrace {\rho \,d\rho \,d\varphi \,d}^{dV}} }}{{\int_{\,0}^{\,d\varphi } \int_{\,r\, - \,t/2}^{\,r\, + \,t/2} {\gamma \underbrace {\rho d\rho d\varphi }_{dA} d} }}\, = \,\frac{{\gamma d\,d\varphi \,\int_{\,r\, - \,t/2}^{\,r\, + \,t/2} {\rho^{2} d\rho } }}{{\gamma d\,d\varphi \int_{\,r\, - \,t/2}^{\,r\, + \,t/2} {\rho \,d\rho \,} }}\\ &= \,\frac{{\gamma td\,rd\varphi \,r\,\left( {1\, + \,\eta^{2} /12} \right)}}{{\underbrace {\gamma td}_{w}\,\,r\,d\varphi }} \hfill \\ &\,\,\,\,\, = \,r\,\,\left( {1 + \eta^{2} /12} \right) \hfill \\ \end{aligned} $$

(52)

where Milankovitch [43,44,45] multiplicative correction factor (1 + η²/12) appears and comes then to affect the various governing equations (at growing resulting η).

Specifically [7], for equilibrium relations:

$$ h = h_{\,1} = h_{\,L} = \frac{{\left( {2 - \eta } \right)\,\,\beta \;sin\,\beta - 2\left( {1 - \cos \beta } \right)\left( {1 + \eta^{2} /12} \right)}}{{2 + \eta - \left( {2 - \eta } \right)\cos \beta }} $$

(53)

$$ h = h_{\,2} = A - \frac{2}{2 + \eta }\,\,\left( {1 + \eta^{2} /12} \right) = \frac{{\left( {2 + \eta } \right)\,\,\alpha \;sin\;\alpha - 2\left( {1 - \cos \alpha } \right)\left( {1 + \eta^{2} /12} \right)}}{{\left( {2 + \eta } \right)\,\,\left( {1 - \cos \alpha } \right)}} $$

(54)

eccentricity relations:

$$ \hat{e}\text{(}\beta \text{)} = \frac{{2e\text{(}\beta \text{)}}}{t} = \,\,\frac{{2\beta \sin \beta - 2\text{(}1 - \cos \beta \text{)}\,\left( {1 + \eta^{2} \text{/}12} \right) - h\,\text{(}2 + \eta - 2\cos \beta \text{)}}}{{\eta \;\text{(}\beta \sin \beta + h\cos \beta \text{)}}} $$

(55)

$$ \hat{e}\text{(}\beta \text{)}\, = \,\frac{{2e\text{(}\beta \text{)}}}{t}\, = \,\frac{{2\text{(}2 + \eta \text{)}\,\left( {\beta \sin \beta - (1 - \cos \beta )A} \right) - \eta \left( {\text{(}2 + \eta \text{)}A - 2\cos \beta \,\left( {1 + \eta^{2} \text{/}12} \right)} \right)}}{{\eta \text{(}2 + \eta \text{)}\,\left( {\beta \sin \beta - \text{(}1 - \cos \beta \text{)}A} \right) + \eta \left( {\text{(}2 + \eta \text{)}A - 2\cos \beta \,\left( {1 + \eta^{2} \text{/}12} \right)} \right)}} $$

(56)

and tangency condition:

$$ h = h_{e} = \,\,\frac{{\text{(}2 - \eta \text{)}\,\left( {\beta \sin \beta } \right)^{\prime } - 2\sin \beta \,\,\left( {1 + \eta^{2} \text{/}12} \right)}}{{\text{(}2 - \eta \text{)}\sin \beta }} = \overbrace {\beta \,\cot \,\beta }^{{h_{H} }} + \left( {1 - \frac{2}{2 - \eta }\,\,\left( {1 + \eta^{2} \text{/}12} \right)} \right) $$

(57)

leading to final Milankovitch governing system:

$$ \left\{ \begin{aligned} &h = h_{1} = \frac{{\left( {2 - \eta } \right)\,\beta \sin \beta - 2\left( {1 - \cos \beta } \right)\,\,\left( {1 + \eta^{2} /12} \right)}}{{2 + \eta - \left( {2 - \eta } \right)\cos \beta }}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,^{{\,\,\,\,\,\,\,\,\,^{{}} }} \hfill \\ & \, \hfill \\ &h = h_{2} = A - \frac{2}{2 + \eta }\,\,\left( {1 + \eta^{2} /12} \right)\qquad\qquad\quad\quad\quad Milankovitch\, \hfill \\ &\hfill \\ &h = h_{e} \, = \,\frac{{\text{(}2 - \eta \text{)}\,\text{(}\sin \beta + \beta \cos \beta \text{)} - \,2\sin \beta \,\,\left( {1 + \eta^{2} /12} \right)}}{{\text{(}2 - \eta \text{)}\sin \beta }}\\&\qquad\quad\quad\;= \overbrace {\beta \cot \beta }^{{h_{H} }} + \left( {1 - \frac{2}{2 - \eta }\,\,\left( {1 + \eta^{2} \text{/}12} \right)} \right)\, \hfill \\ \end{aligned} \right. $$

(58)

then configurating a more involved “cubic algebraic problem” [7]:

$$ \begin{aligned} & S^{2} (3g - 2S)\,\,A^{3} + 3\,(g - f)\,g\,S\,\,A^{2} - 3fg^{2} \,A + 2g^{3} = 0\,\,;\, \hfill \\ & S\,\,\eta^{3} + 3\,(f + g)\,\,\eta^{2} - 12\,(g - S)\,\,\eta + 12\,(g - f) = 0\,\,; \hfill \\ & 6S^{2} \,h^{3} - 3S\,(3f - g - 2S)\,\,h^{2} - 3\,(g - f)(f - 2S)\,\,h + 2\,(g - f)^{2} = 0, \hfill \\ & 2S^{3} \hat{h}^{3} - 3(g - f)(f + g + 2S)S\hat{h}^{2} + 6(g - f)^{2} (f + 2S)\hat{h} - 8(g - f)^{3} = 0 \hfill \\ \end{aligned} $$

(59)

Explicit analytical closed-form representations of the solution of such a cubic algebraic problem are now newly provided in [14]. Minimal differences between Milankovitch and CCR solutions may be appreciated at increasing opening angle of the arch, mainly for over-complete arches (Fig. 16).

Fig. 16

CCR vs. Milankovitch solution. Functional dependence of A, η, h, on β

Finally, similarly to previous Figs. 12, 13, 14 and 15, Fig. 17 represents a resuming analytical representation of the line of thrust on the true arch profile in the critical least-thickness condition, for a taken reference case of over-complete arches (α = 140°). Thereby, the salient differences between Heyman, CCR and Milankovitch solutions may be appreciated, all together, as: line of thrust actually going out of the arch profile for Heyman solution; line of thrust truly tangent to the intrados at the haunch for CCR and Milankovitch solutions; quite near representations for CCR and Milankovitch solutions; different angular inner-hinge positions, drastically dissimilar for Heyman solution; different resulting thickness estimates, with Heyman visibly being sub-critical and CCR just slightly under-conservative with respect to Milankovitch solution; positioning of the resulting self-weight resultant for CCR (and Heyman) solution, as opposed to the true location completely accounted for by Milankovitch solution.

Fig. 17

Analytical plot of the lines of thrust for Heyman, CCR and Milankovitch solutions at α = 7π/9 = 140°(A = 0.889347)

Further illustration on the characteristics of the mechanical system, for the three solutions, is available in [7, 14], including for the returning and bell-shaped trends of intrinsic non-dimensional horizontal thrust $\hat{h}$ = ηh, with diverging differences for Heyman solution in the dependencies on β and minimal differences in the dependencies on A and α. Specific aspects of the stationarity of these curves are analytically treated in [14], by explicit closed-form solutions, referring to the case of the symmetric circular masonry arch of “maximum horizontal thrust” and of “widest angular inner-hinge position”.

4 Numerical Validation

Independent, numerical-validation approaches are here outlined, for a confrontation to the analytical outcomes and for a further illustration of the symmetric circular masonry arch characteristics, in the optimality condition of minimum thickness.

4.1 Recursive Evaluation of Angular Inner-Hinge Position

4.1.1 Heyman Solution

The numerical solution of transcendental Eq. (44)_a for angular inner-hinge position β is quite straightforward (especially in inverse form A(β) = α cot α/2, at given β). However, a recursive procedure could promptly be devised, e.g. for a further root refinement of a guess that could be taken, for instance, from a proposed root fit of the solution [7].

Table 4 reports a possible convenient way of doing that, which requires origin solutions that are not much on the left of the correct one or, more precisely, that are on the right of a singularity point that may arise when the denominator term in the recursive form, which is a function of β and A, becomes zero. Despite this little limitation (which is promptly over-passed by the very good estimate provided by the above-commented good fit), convergence turns out quite fast and, most important, over all the range of the admissible values of A. Indeed, the recursive proposal reported in Table 4 allows for a fast convergence on both extremes of the values of A, being actually slower in the usually considered range of half-angles of embrace that are lower than π/2. Other proposals might work better in this range but would display a slower convergence on the opposite side or might have no singular points but achieve a much slower convergence for the different values of A. So, the presented option constitutes a reasonable compromise to handle all the possible cases.

Briefly, the adopted recursive formula that has been reported on top of Table 1 has been obtained as follows. Take Eq. (44)_a, isolate a pivoting β_p term and solve with respect to that. Different possibilities arise, which can be checked right away for a possible recursive convergence. The chosen one has originated from the following choice (for compactness, again, S = sin β, C = cos β):

$$ \beta \,\,\frac{{2\beta_{p} \,C/S + C^{2} + 1}}{{2\beta_{p} + S\,C - S}} = A\;\;\; $$

(60)

Table 1 Recursive determination of inner-hinge position for classical Heyman solution

Also, to further improve the convergence rate, additional splits of the pivoting term have been attempted and optimized through numerical trials by the following proposal:

$$ \beta \,\,\frac{{12\,\beta_{p} C/S - 10\beta \,C/S + C^{2} + 1}}{{6\beta_{p} - 4\beta + S\,C - S}} = A $$

(61)

which, by solving with respect to β_p, leads to the following expression, finally useful towards a recursive evaluation of root β in Heyman solution:

$$ \beta_{iter}^{H} = \frac{{\beta \,\left(10\beta \,C - S\,\text{(}1 + C^{2} )\right) - S\,\left(4\beta + S\text{(}1 - C)\right)\,A}}{{6\text{(}2\beta \,C - S\,A\text{)}}} $$

(62)

This expression has been used to generate the recursive estimations reported in Table 1, for five given values of A (A = α = π/2 and two values on the right and on the left of that). Starting from simple fitting guess [7]:

$$ \beta_{0}^{H} \, = \sqrt {3} \,\,\,\left( {2 - A} \right)^{1/2} \,\,\left(1 - 1\,/6\,\text{(}2 - A\text{)}^{2/3} \right) $$

(63)

a few iterations turn out enough to recover the correct recursive estimate of root β.

4.1.2 CCR Solution

Expressions (47) on A could numerically be used for a recursive determination of the value of A at given pre-peak 0 ≤ β ≤ β_sβ. Specifically, the first two relations in the first line of Eq. (47) would lead to root values A⁺ ≥ A_sβ and the second two to root values A⁻ ≤ A_sβ. Indeed, as it has been illustrated, according to Eq. (51)_a, function A(β) is double-valued, at a given value of 0 ≤ β ≤ β_sβ. To let successfully run the recursive iteration in all four cases reported in Eqs. (47), initial roots A₀ should be chosen as defect estimates for the A⁺ root and as excess estimates for the A⁻ root.

Similarly, Eqs. (48) and (49) could be used as well for a recursive determination of η and h, with above comments applying in the same way. However, recall that the role of η is inverted, as opposed to that of A and h, since the triplets have to be sorted out in orders (A⁺,η⁻, h⁺) and (A⁻,η⁺, h⁻). The estimate of $\hat{h}$ = ηh might either go through similar resolutions of Eqs. (50) or directly by the product of the found η and h estimates.

A root recursive estimate of angular inner-hinge position β that works quite well on all sides of the solution branches is provided in Table 2, similarly to what reported in Table 1 for Heyman solution. It is based on the first expression in the second line of Eq. (47), which can lead to:

$$ \beta_{iter}^{CCR} \,\, = \,\,\frac{{2 \, fg\,A - S\text{(}2g - S\text{)}A^{2} }}{g} \, - \, S\,C\,\,\, = \,\,\,2fA - 2\,S\,A^{2} + \frac{{S^{2} A^{2} }}{g} - \, S\,C $$

(64)

Table 2 Recursive determination of inner-hinge position for CCR solution

Starting from easy-to-remember CCR fitting guess [7]

$$ \beta_{0}^{CCR} \, = 3\,/2\,\,\left( {A - 2\,/3} \right)^{1/4} \,\left( {2 - A} \right)^{1/2} $$

(65)

which one could derive for CCR solution based on the trends experienced for α small (A near 2) and α near α^CCR_l (A near 2/3), see Fig. 4, a good refinement is achieved in not many iterations (in general terms, a bit more than for the previous Heyman recursive evaluation).

4.1.3 Milankovitch Solution

As above, Table 3 presents a recursive strategy for refining angular hinge position β in Milankovitch solution, according to expressions

$$ \begin{aligned} \beta_{iter}^{M\,} &= \,\frac{{ - 2g^{3} + 3fg^{2} A + 3fg\,S\,A^{2} + S^{2} \text{(}2S - 3g\text{)}A^{3} }}{{3g\,S\,A^{2} }}\, \, - \,\,S\,C \hfill \\ & = f - SA + \frac{{2S^{2} A}}{3g\,\,} + \frac{\,f\,g\,}{S\,A} - \frac{{2g^{2} }}{{3\,S\,A^{2} }}\,\, - \,\,S\,C \hfill \\ \end{aligned} $$

(66)

starting from an appropriate fitting guess, such as [7]

$$ \beta_{0}^{M} \, = 1.53\,\,\,\,\left( {A - ({\sqrt 3 - 1})} \right)^{1/4} \,\left( {2 - A} \right)^{1/2} $$

(67)

leading to a reasonable convergence on all branches of the solution characteristics.

Table 3 Recursive determination of inner-hinge position for Milankovitch solution

Overall, the number of iterations may become the highest, among the three solution instances, due to the increasing complexity in, respectively, dealing with Heyman, CCR and Milankovitch solutions (appreciate the global increasing height of represented Tables 1, 2 and 3).

4.2 DDA Least-Thickness Results

A least-thickness self-standing evaluation of the masonry arch may be elaborated by using Discrete Element Method (DEM) quasi-static simulations of discrete voussoir arches [36,37,38]. To provide an independent numerical validation of the achieved analytical results, an available Discontinuous Deformation Analysis (DDA) tool was adopted in [8, 9], to deliver the estimates of the critical thickness and the appearance of the corresponding collapse mode (notice that the five-hinge purely rotational collapse mode is assumed from scratch, in the analytical analysis, while in such a case is numerically evaluated, out of the analysis). Here, further complementary and completing results are reported. On the description of the employed methodology, and the framing in the competent literature, see [9]. The adopted DDA computational tool was freely taken from the web (sourcefoce.net, “DDA for Windows”, Limerick version 1.6), as developed from researchers at the University of Berkeley (see [15,16,17]).

Symmetric discretized arches with four blocks (with radial joints), at variable half-opening angles α between 60° and 140°, each 10° (thus encompassing under-complete, semicircular complete and over-complete arches) are DDA analysed, at given inner-joint position, for which the minimum thickness is numerically estimated (alias, the minimum thickness still preventing, or the maximum thickness still inducing, arch collapse). At each given value of α, a target value of angular inner-joint position β is assumed, on the basis of the derived analytical solution (CCR solution is here taken for reference). For each arch, three analyses have been performed, for that value of β and for two fork, inferior/superior values of β, as ± 0.5° that value of β, as to reveal and confirm a possible maximum trend of η at variable β, as it should be in denouncing the critical least-thickness condition [9].

A summary of such DDA results, with a direct comparison to CCR and Milankovitch solutions is provided in Table 4 and compactly illustrated in Figs. 18 and 19.

Table 4 Summary of numerical DDA result on critical thickness, with comparison to CCR and Milankovitch solutions

Fig. 18

Numerical DDA results for critical η and relevant rotational collapse mode

Fig. 19

Numerical DDA results and comparison to CCR and Milankovitch solutions

Table 4 reports the recorded values of critical η = t/r, at the given values of α and β, determined with the procedure explained in [9]. On the side, the target discrete η values for CCR and Milankovitch solutions [7] are as well reported, for comparison purposes. As it may be appreciated, differences look really minimal and results show that the features of the analytical treatment are correctly reproduced and that the numerical DDA tool is able to provide consistent estimates of the critical thickness (having reasonably guessed the position of the inner-joint). Clearly, little dislocations of the inner-joint position do not alter much the numerically recorded values of critical thickness.

Similar comments may be deciphered from the reading of Figs. 18 and 19. Figure 18 globally shows in grid-view the array of the various masonry-arch openings, with inner-joint set from CCR solution, with the corresponding recorded value of critical thickness-to-radius ratio and the attached apparent rotational collapse mode. The sequence of plots clearly illustrates the monotonic increasing trend of thickness necessary for the arch to withstand, at increasing arch opening, over-complete arches becoming really much thicker to stand up. The purely rotational collapse mode is anyway correctly reproduced (a high value of friction coefficient is set in the simulations, to avoid possible manifestations of any form of sliding [9], se also next subsection). This looks thus consistent with source Heyman hypothesis of no sliding failure, at the basis of the present theoretical treatment (finite friction effects are separately analysed in [10,11,12,13]).

Figure 19 further gathers, on the top plot, the imposed values of angular inner-joint position β at variable half-opening values α of the arch, with fork values around CCR solution (Milankovitch solution is never much dissimilar, if not, a bit, for the last case of α = 140°) and going through a maximum of β for an α at around 120°–130° [14]. Just below, in the centre plot, the recorded critical values of η are plotted, as a function of imposed β, with almost undistinguishable η results for the considered β positions. This scored plot provides a sort of discretized functional representation of the continuous plots earlier depicted in Fig. 16_b (comparison to Figs. 10_b shall also confirm an increasing divergence to Heyman solution, also for the DDA results). The bottom plot finally represents the DDA trends of η as a function of α, where clearly the differences from the underlying values in β position really disappear in terms of η. The plot truly confirms a continuous, increasing trend of critical thickness, at increasing arch opening, as made apparent in looking, all together, at the array of arches gathered in grid-view in previous Fig. 18.

4.3 Least-Thickness Optimization by a Complementarity Problem/Mathematical Programming Formulation

A Complementarity Problem/Mathematical Programming (CM/MP) numerical formulation and self-made implementation, recently developed by the authors [12, 13] within a MATLAB environment, with the target to specifically enquire finite-friction effects on masonry arches [22, 53,54,55,56,57,58,59,60,61,62], is here adapted to the numerical analysis of symmetric circular masonry arches relying on infinite (say high) friction and employed for a further independent validation and interpretation of the arch characteristics in the least-thickness condition, as by the solution of a numerical optimization problem.

The general formulation is first briefly introduced, in its main traits, and described in the needs of adaptation of the computational implementation, towards the present numerical optimization analysis. Then, it is run, and salient results are selected and presented, in view of complementing and confronting the previous analytical and numerical outcomes.

The formulation sets in within similar modelization attempts in the limit analysis of masonry constructions [63,64,65,66,67,68,69,70,71,72,73,74] and relies on a series of works in stating the optimization problem as a Complementarity/Mathematical Programming Problem [75,76,77,78,79,80], as ruled by typical theoretical literature references [81,82,83,84,85].

The masonry arch is subdivided into n rigid blocks, by (n + 1) radial joints (Fig. 20) and the relative voussoir kinematics at each joint is described by 4 non-negative kinematic variables (ϑ⁺, ϑ⁻, s⁺, s⁻; see Fig. 21) [13].

Fig. 20

Geometric representation of a general circular masonry arch (with thickness-to-radius ratio η = t/r), with redundant reactions at the right shoulder, internal actions within the circular masonry arch and block weight, in the CP/MP formulation

Fig. 21

Kinematic description and variables of possible relative (rotational and sliding) movements among adjacent rigid masonry blocks, in the CP/MP formulation

The 3n velocities ruling the rigid-body movements of the n rigid blocks are collected in vector ${\dot{\mathbf{U}}}_{G}$:

$$ {\dot{\mathbf{U}}}_{G} = \left[ {\dot{u}_{{G_{1} }} ,\dot{v}_{{G_{1} }} ,\dot{\psi }_{1} , \ldots ,\dot{u}_{{G_{k} }} ,\dot{v}_{{G_{k} }} ,\dot{\psi }_{k} , \ldots ,\dot{u}_{{G_{n} }} ,\dot{v}_{{G_{n} }} ,\dot{\psi }_{n} } \right]^{T} $$

(68)

where $\dot{u}_{{G_{i} }} ,\dot{v}_{{G_{i} }} ,\dot{\psi }_{i} ,\,\,\,\,i = 1 \ldots n$, are the linear and angular velocities of each i-th block and the 4(n + 1) relative velocities relevant to the (n + 1) joints are gathered in vector ${\dot{\varvec{\lambda }}}$:

$$ {\dot{\varvec{\lambda }}} = \left[ {\dot{\vartheta }_{1}^{ + } ,\dot{\vartheta }_{1}^{ - } ,\dot{s}_{1}^{ + } ,\dot{s}_{1}^{ - } , \ldots ,\dot{\vartheta }_{k}^{ + } ,\dot{\vartheta }_{k}^{ - } ,\dot{s}_{k}^{ + } ,\dot{s}_{k}^{ - } , \ldots ,\dot{\vartheta }_{n + 1}^{ + } ,\dot{\vartheta }_{n + 1}^{ - } ,\dot{s}_{n + 1}^{ + } ,\dot{s}_{n + 1}^{ - } } \right]^{T} $$

(69)

The two sets of kinematic variables are linearly related by compatibility matrix ${\mathbf{ B}}_{G}$:

$$ {\dot{\mathbf{U}}}_{G} {\mathbf{ = B}}_{G} \,{\dot{\varvec{\lambda }}}\text{,}\;\;\;\;\;\;\;\;{\dot{\varvec{\lambda }}} \ge {\mathbf{0}} $$

(70)

The physical constraints on the masonry arch require the fulfilment of the following kinematic relations:

$$ {\mathbf{B}}\,\dot{\varvec{\lambda }}\, =\, {\mathbf{0}}\text{,}\;\;\;\;\;\;\;\;{\dot{\varvec{\lambda }}} \ge {\mathbf{0}} $$

(71)

where B is a (kinematic) constraint matrix [13] while, at each joint, internal actions must fulfil the following static inequalities:

$$ {{\varvec{\varphi}} = \left[ {\begin{array}{*{20}c} {{\varvec{\varphi }}_{1} } \\ \cdots \\ {{\varvec{\varphi }}_{k} } \\ \cdots \\ {{\varvec{\varphi }}_{n + 1} } \\ \end{array} } \right] \le {\mathbf{0}},}\;\;\;\;\;\;\;\;\;\;\;\;\;{\varvec{\varphi }}_{k} = \left[ {\begin{array}{*{20}c} {\varphi_{{\vartheta^{ + } }}^{k} }\\ {\varphi_{{\vartheta^{ - } }}^{k} } \\ {\varphi_{{s^{ + } }}^{k} } \\ {\varphi_{{s^{ - } }}^{k} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {\;\;\,M_{k} - N_{k} {\displaystyle{\frac{t}{2}}}} \\ { - M_{k} - N_{k} {\displaystyle{\frac{t}{2}}}} \\ {\;\;\;\;T_{k} - \mu N_{k} } \\ {\; - T_{k} - \mu N_{k} } \\ \end{array} } \right] \le {\mathbf{0}} $$

(72)

Once internal actions are computed by equilibrium, relationships (72), in terms of activation functions φ^k, define an “admissible static configuration” within the masonry arch, the independent variables being three redundant reactions H, V and W at the right built-in shoulder of the arch (Fig. 20):

$$ {\varvec{\varphi}} = {\mathbf A}\,{\mathbf{H}} + {\mathbf T}^{w} \le {\mathbf{0}},\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathbf{H}} = \left[ {\begin{array}{*{20}c} H \\ V \\ W \\ \end{array} } \right] $$

(73)

where A and T^w are, respectively, the matrix and the vector governing equilibrium and joint activation [13].

Each relative velocity shall result orthogonal to the corresponding static activation function:

$$ \varphi_{{\vartheta^{ + } }} \,\dot{\vartheta }^{ + }=0\;\;\;\;\;\;\varphi_{{\vartheta^{ - } }} \,\dot{\vartheta }^{ - }=0\;\;\;\;\;\;\varphi_{{s^{ + } }} \,\dot{s}^{ + }=0\;\;\;\;\;\;\varphi_{{s^{ - } }} \,\dot{s}^{ - }=0 $$

(74)

Such relations also include a complete detachment (see [13]), for which, in order to eliminate a numerical multiplicity (not corresponding to a physical one), a (weak) “orthogonality” condition among non-negative variables $\dot{s}^{ + }$ and $\dot{s}^{ - }$ shall be added to the problem statement:

$$ \dot{s}^{ + } \;\dot{s}^{ - } {\mathbf{ = }}\text{0} $$

(75)

The external power is then readily computed as:

$$ \dot{L}_{e} \text{ = }{\mathbf{w}}^{T} {\dot{\mathbf{U}}}_{G} \ge \text{0,}\;\;\;\;\;\;\;\;\;\;{\mathbf{w}}\text{ = }\left[ {\begin{array}{*{20}c} {{\mathbf{w}}_{1} } \\ \cdots \\ {{\mathbf{w}}_{k} } \\ \cdots \\ {{\mathbf{w}}_{n} } \\ \end{array} } \right],\;\;\;\;\;\;\;\;\;\;\;{\mathbf{w}}_{k} = \left[ {\begin{array}{*{20}c} 0 \\ { - \left( {\beta_{k + 1}^{*} - \beta_{k}^{*} } \right)\gamma \,t\,d\;r} \\ 0 \\ \end{array} } \right] $$

(76)

Finally, the velocity field can arbitrarily be normalized by setting:

$$ {\mathbf{p}}^{T} {\dot{\varvec{\lambda }}} = \text{1,}\;\;\;\;\;\;\;\;\;\;p\left( k \right) = 1,\;\;\;\;\;\;\;\;\;k = 1, \ldots ,n + 1 $$

(77)

In general terms, in a limit equilibrium configuration at incipient collapse, kinematic (${\dot{\varvec{\lambda }}}$) and static (H) variables shall fulfil the following linear Complementarity Problem (see, e.g. [77]):

$$ \left\{ {\begin{array}{*{20}c} {f = - {\varvec{\varphi }}^{T} {\dot{\varvec{\lambda }}} = 0} \\ {{\varvec{\varphi }} = {\mathbf{A}}\,{\mathbf{H}} + {\mathbf{T}}^{W} \le {\mathbf{0}}} \\ {{\mathbf{B}}\;{\dot{\varvec{\lambda }}} = {\mathbf{0}}} \\ {\dot{L}_{e} \text{ = }{\mathbf{w}}^{T} {\mathbf{B}}_{G} \,{\dot{\varvec{\lambda }}} \ge \text{0}} \\ {{\mathbf{p}}^{T} {\dot{\varvec{\lambda }}} = \text{1}} \\ {{\dot{\varvec{\lambda }}} \ge {\mathbf{0}}} \\ {{\dot{\mathbf{s}}}^{ + T} \,{\dot{\mathbf{s}}}^{ - } = {\dot{\varvec{\lambda }}}^{T} \,{\mathbf{B}}_{s}^{ + T} \,{\mathbf{B}}_{s}^{ - } \,{\dot{\varvec{\lambda }}} = 0} \\ \end{array} } \right. $$

(78)

Additionally, a convenient solution to complementarity system (78) may be obtained by the following non-linear Mathematical Programming problem, in which the (nonlinear) orthogonality condition on variables ${\varvec{\varphi }}$ and ${\dot{\varvec{\lambda }}}$ is used as an objective function (see, e.g. [78]), to be led to a zero value:

$$ f_{\hbox{min} } = \mathop {min}\limits_{{{\dot{\varvec{\lambda }}},{\mathbf{H}}}} \left\{ { - {\varvec{\varphi }}^{T} {\dot{\varvec{\lambda }}}\left| {\begin{array}{*{20}c} {{\varvec{\varphi }} = {\mathbf{A}}\,{\mathbf{H}} + {\mathbf{T}}^{W} \le {\mathbf{0}}} \\ {{\mathbf{B}}\;{\dot{\varvec{\lambda }}} = {\mathbf{0}}} \\ {\dot{L}_{e} \text{ = }{\mathbf{w}}^{T} {\mathbf{B}}_{G} \,{\dot{\varvec{\lambda }}} \ge \text{0}} \\ {{\mathbf{p}}^{T} {\dot{\varvec{\lambda }}} = \text{1}} \\ {{\dot{\varvec{\lambda }}} \ge {\mathbf{0}}} \\ {{\dot{\mathbf{s}}}^{ + T} \,{\dot{\mathbf{s}}}^{ - } = {\dot{\varvec{\lambda }}}^{T} \,{\mathbf{B}}_{s}^{ + T} \,{\mathbf{B}}_{s}^{ - } \,{\dot{\varvec{\lambda }}} = 0} \\ \end{array} } \right.} \right\} = 0 $$

(79)

Further, also the orthogonality condition on variables ${\dot{\mathbf{s}}}^{ + }$ and ${\dot{\mathbf{s}}}^{ - }$ may conveniently be transferred to the objective function, by adding (non-negative) scalar product ${\dot{\mathbf{s}}}^{ + T} \,{\dot{\mathbf{s}}}^{ - } = 0$ to (non-negative) term $- {\varvec{\varphi }}^{T} {\dot{\varvec{\lambda }}}$, in this way leading to a non-linear Mathematical Programming problem with linear constraints only [13].

Notice that, in Mathematical Programming problem (79), thickness t and friction coefficient μ are assumed as free parameters to be arbitrarily changed, to get to a zero value as an optimum value for objective function $f = - {\varvec{\varphi }}^{T} {\dot{\varvec{\lambda }}}$.

Multiple solutions to (non-convex) problem (78), and thus to programming problem (79), may generally be expected (see, e.g. [79] and [82,83,84,85]), as linked to the context of non-standard Limit Analysis, in the realm of masonry arches with finite friction effects [53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73].

Additional details on the source CP/MP formulation are delivered in [13]. After an appropriate normalization of variables, the solution of the optimization problem has been obtained within MATLAB^®, by either the “interior point” or the “active set” minimization algorithm in (built-in) optimization function “fmincon”, with tolerances kept at 10⁻¹⁰. The following dedicated strategies have here been considered and implemented, to handle the specific problem at hand, and leading to the outcomes then gathered in Figs. 22 and 23:

Fig. 22

Grid-view of numerical CP/MP results (referred to CCR solution)

Fig. 23

Numerical CP/MP arch optimization (referred to CCR solution): maximum horizontal thrust (top); widest inner-hinge position from the crown (bottom)

• In a way to let all arch characteristics to be extracted in output from the numerical implementation, a procedure has been devised to record out the value of intrinsic non-dimensional horizontal thrust ĥ, out of the normalization of static primary-variable redundant reaction H (Fig. 20), at the achieved solution instance. The accuracy on the determination of such a static variable has been monitored and found as almost comparable to that recorded for the source geometrical variables in the optimization problem. This appears as a rather consistent feature of the present adaptation implementation, since it comes up able to correctly reproduce all the characteristic features of the least-thickness arch (geometrical, kinematical and statical). Specifically, the quest was posed in terms of evaluating the maximum horizontal thrust that the arch is able to transfer in the minimum-thickness condition (Fig. 23), as explicitly analytically addressed in [14]. The recorded matching, on the estimation of such a maximum value of thrust, turns out astonishingly good.

• Issue of symmetry. As briefly introduced, the present formulation for circular masonry arches is rather general and allows for the analysis of non-symmetric arches (Fig. 20). Moreover, the collapse mechanism to be plotted is obtained out of the numerical calculation on the kinematic variables and, for symmetric arches, is by no means pre-imposed to be symmetric [13]. In other words, what becomes uniquely know is the geometrical position of the failure joints (geometrical variables), while the mechanism undertakes an arbitrariness linked to an exuberant number of dofs. By instead implicitly imposing symmetry, for the resulting mechanism, this becomes a single-dof mode, and fully respecting symmetry, with a transverse zero displacement at the crown of the symmetric masonry arch. This has been set in the present adaptation, through the insertion of a specific control flag, within the implementation, up to pre-set the symmetry condition in the treatment of the geometrical reference configuration and of the kinematic variables of the formulation and, for the definition of the symmetric collapse mode, some equations are added to set the relations for the centres of gravity of the i-th and j-th homologous chunks (by the axis of vertical symmetry at the crown), namely u_i + u_j = 0, v_i − v_j = 0, φ_i + φ_j = 0. If the number of chunks is odd, for central chunk r including the crown, the following symmetry constraints are considered, for a consistent symmetric arch kinematics: u_r = 0, φ_r = 0.

• Considering, still, the representation of the (symmetric) collapse mode (see Figs. 22 and 23), the following dedicated strategies have been implemented, in view of achieving a convincing and realistic representation out of the act of motion within linear kinematics (displacements for velocities). The positions of the vertices of the masonry block are determined through linear kinematics, thus conserving compatibility among the chunks and with external constraints. Moreover, the relative angle among the joints of the voussoir is the same, as the rigid-body act of motion lefts it unaltered, still within an arbitrarily amplified, linear kinematics (i.e. it represents a conformal mapping). Intrados and extrados vertices are, respectively, joined by circular segments, with a slightly varied radius, given by the prolongation of the facing joints of each chunk.

• Since finite-friction effects are here not of a target, the information from the “line of friction” [12, 13] has been eliminated, and a high value of friction coefficient µ among the blocks has been set, up to prevent any form of sliding, and then warrant a purely rotational collapse mode, as indeed recorded in all the analysed cases with 60° ≤ α ≤ 140° (Figs. 22 and 23), where a high value of µ = 10 (friction angle φ = 84.3°) has been imposed. Indeed, any form of sliding, for α ≤ α_lm = 2.48716 = 142.504°, shall be prevented by μ ≥ μ_lm= 1.41527 (φ_lm= 54.7558°) [11], as also confirmed by the numerical analyses in [13].

• Given the fact that the amount of friction is here immaterial, with respect to what advanced in [13], a single plot with the line of thrust, turning out strictly contained within the arch, in the critical, optimality configuration of least-thickness, having superimposed the resulting purely rotational collapse mode, has been set, also in the perspective to facilitate the unitary grid-view representation in Fig. 22. Thus, the algorithm has been “specialized” to conveniently represent the features of symmetry-constrained arches, in the hypothesis of “infinite” friction, for inquiring and representing purely rotational collapse, as consistent with the mathematical derivation in the first part of the paper, for the continuous masonry arch.

Numerical results for symmetric arches with four blocks are condensed in Figs. 22 and 23 (here, comparison reference should be made to CCR solution, since uniform self-weight distribution is taken along geometrical centreline). Figure 22 shows, in compact grid-view, as earlier done in Fig. 18 for the DDA treatment, the results for the same half-angles of embrace of the arch, from 60° to 140°, each 10°. Results are perfectly matching, also with the numerical values that shall be compared with the analytical formulation (CCR solution). Figure 23 further shows two specific cases of a stationarity interest, namely the arch with “maximum horizontal thrust” and with “widest angular inner-hinge position”, made the subject of a specific analytical investigation in [14], with wholly consistent results.

5 Conclusions

In the first part of the paper, different, equivalent analytical derivation approaches of the ruling equations for the statics of circular masonry arches in the least-thickness, optimality condition have been provided, specifically those based on the so-called upper and lower horizontal thrusts (“Coulomb’s static approach”) and on the balance of virtual work or power (“Mascheroni’s kinematic approach”), both leading to same optimization outcomes:

• The proposed formulations reconcile previous treatments and bring in a unifying and additional interpretation and understanding of the conditions that rule the least-thickness solution.

• Specifically, beside established equilibrium at incipient collapse, the tangency condition of the line of thrust (locus of pressure points) at the intrados of the arch can be handled, as a stationary condition on the horizontal thrust itself, without the need of going through the definition of the line of thrust, namely of the equation of its eccentricity with respect to the geometrical centreline of the arch.

• The generalization to Milankovitch solution for the real self-weight distribution is also reconsidered, and the features of the three arising solutions further described, with an additional illustration. Thus, analytical features and results herein presented complement those earlier provided in previous source companion work [7].

In the second part of the manuscript, an independent numerical validation investigation has then been developed, by three separate approaches, with truly consistent results:

• First, a recursive determination of angular inner-hinge position of a continuous arch is devised, to avoid the direct numerical solution of the governing system of equations. This could conveniently be handled even by hand computations, or by a simple programming calculator or a spreadsheet, leading to the determination of the angular inner-hinge position, after which the main features of the arch characteristics, in the least-thickness critical condition, could directly flow down, by substitution into the provided explicit analytical representations.

• Second, a further numerical Discrete Element Method (DEM), Discontinuous Deformation Analysis (DDA) investigation of discrete voussoir arches has been developed, to deliver a confirmation of the assumed purely rotational collapse mechanism and relevant estimation of the critical thickness. The achieved results complement and complete those earlier presented in [9].

• Third, an innovative Complementarity Problem/Mathematical Programming (CP/MP) formulation has been adapted and operated on symmetric circular masonry arches at high friction, for a consistent, final confirmation of all the characteristics of the circular masonry arch in the least-thickness condition. The attached presented results complement those delivered in [12, 13].

The current attempt has considered the so-called Couplet–Heyman problem in the statement of a classical form-optimization problem in the Mechanics (statics) of (symmetric) circular masonry arches, as pertinent to the present editorial collocation. The optimization process has been set as follows:

• Initially, by analytical procedures, through the appropriate statement of the underlying stationary conditions, as those corresponding to Heyman literal description, and has allowed to derive closed-form explicit analytical representations, for the problem at hand.

• Additionally, by numerical strategies, where the optimization problem has been scheduled as a direct numerical target (objective function), which opens up the way to further perspectives towards a general interpretation and practical application, for cases that shall be going beyond those codified ones that may conveniently be handled by analytical derivations.

Indeed, further subsequent investigations may consider different shapes of the masonry arch and possible true implications of the presence of finite friction, in terms of non-standard Limit Analysis, about the bearing capacity of the masonry arch and the resulting collapse mode, which may include sliding and other potential effects linked to the arising non-normality kinematics of the arch collapse, as initially investigated, analytically in [10, 11], and numerically in [12, 13].