1 Introduction

The Mechanics of masonry arches constitutes a fascinating theme of research, even in present times, towards the understanding of general and specific functional principles, and resulting features, for this peculiar type of natural and classical form of bearing structures, both from the theoretical and the practical points of view.

Specifically, codified geometrical shapes of the masonry arch may be considered, as attached to typical historical architectural forms. Among those, the shape of the (symmetric) circular arch is that which not only has been typically employed, since ancient constructions but also allows for more convenient inspection and analytical analysis.

In that realm, several rational studies in “modern” times have considered the investigation of specific critical (or optimum) conditions for the masonry arch, as that leading to the least-thickness form optimization under self-weight, as “recently” codified in the pioneering mechanical works by Jacques Heyman [1,2,3,4,5,6]. According to his three classical hypotheses of arch behaviour (namely: no tensile strength; infinite compressive strength; no sliding failure), solutions are sought for the typical characteristics at incipient (purely rotational) collapse, namely critical thickness t (to-radius r ratio, \(\eta =t/r\)), angular position of the inner hinge of the collapse mechanism (from the crown) \(\beta \) and (non-dimensional) horizontal thrust H acting within the masonry arch (\(h=H/(w r)\), \(w=\gamma t d\), \(\gamma \) specific weight, d out-of-plane depth of the arch, or intrinsic \(\hat{h}={H/(\gamma d r^2)=\eta h}\), as here further defined in Sect. 3).

Within a previous mainstream of work by the present authors [7,8,9,10,11,12,13,14,15], such a classical optimization problem in the Mechanics (statics) of symmetric circular masonry arches has been revisited, and framed within the relevant, now updated, literature [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34], by providing new, explicit analytical derivations and representations of the solution characteristics. Specifically, classical Heyman solution [3] has been shown to constitute a sort of approximation of the true solution (here labelled as “CCR” [7]), in Heyman assumption of uniform self-weight distribution (location of the centres of gravity along the geometrical centreline of the circular arch), while Milankovitch solution [35,36,37] may as well be consistently derived, in the consideration of the true self-weight distribution, though at the price of an increasing complexity in the explicit analytical handling of the governing equations (here newly resolved to a very end).

Indeed, Heyman, CCR and Milankovitch solutions have been shown to arise from the solution of a first-order (linear), a second-order (quadratic) and a third-order (cubic) algebraic problem, respectively, with resulting outcomes that are all rather similar, in terms of “external” arch characteristics (thickness-to-radius ratio \(\eta \) and non-dimensional horizontal thrust h), though differing, in Heyman’s case, about the prediction of “internal” arch feature angular inner-hinge position \(\beta \). In other words, CCR and Milankovitch solutions turn out to result very close to each other, for all three basic arch characteristics \(\beta , \eta , h\); Heyman solution is still able to locate the main features of the self-standing arch (and in a much simpler manner), though fails in predicting the correct position of the inner intrados hinge in the associated purely rotational collapse mode, with an increasing divergence at growing opening angle of the arch. This becomes much apparent for over-complete (horseshoe) arches, which constitute one rather emblematic arch shape in historic constructions, typically in Moorish architecture, with different forms, like rounded, pointed or lobed (the present work considers just the circular arch shape, with radial joints, in the whole analysis).

Within such a framework of research, the present study investigates further specific features of the solution cases, through a complete analytical derivation, up to deliver closed-form explicit mathematical representations of the sought arch characteristics, as developed all together, for the above three solution instances. Moreover, particular interesting extremal conditions are explicitly investigated, as earlier revealed by the CCR analysis in [7], namely those corresponding to the maximum intrinsic non-dimensional horizontal thrust \(\hat{h}\) that the circular masonry arch is able to withstand and transfer in the critical least-thickness condition under self-weight, at given material (\(\gamma \)) and geometrical (d, r) properties, and to the foremost wide angular inner-hinge position \(\beta \) from the crown appearing in the purely rotational mode developing at incipient collapse.

Specifically, the present analytical note newly:

  • Considers a complete description of the governing equations of the “Couplet-Heyman problem” through appropriate on/off control \(\delta \) flags, allowing to automatically shift among the three solution instances;

  • Develops a unitary representation of the up-to-cubic polynomials governing the resulting algebraic problem;

  • Provides explicit analytical solutions with closed-form expressions, markedly also for the cubic algebraic problem attached to the true Milankovitch self-weight distribution;

  • Examines the stationary (maximum) condition of the intrinsic non-dimensional horizontal thrust of the arch, leading to define the “arch of maximum thrust” in the critical least-thickness condition;

  • Derives the stationary (maximum) condition of foremost wide angular inner-hinge position from the crown, by explicitly codifying analytical expressions linked to that.

  • Illustrates the arising results and characteristic features linked to the above steps, in terms of mechanical understanding and practical value.

As per the treatment adopted in the paper, namely a full analytical one, the present study is framed between similar related investigations that have enquired the problem by analytical or semi-analytical approaches, between which specific important developments that shall be mentioned in sharing the vision and perspectives here pursued are those concisely outlined as follows: pioneering and fundamental Heyman contributions [1,2,3,4,5,6]; attempts by the present authors [7,8,9,10,11,12,13,14,15], specifically: analytical [7], analytical-numerical, accompanied by a Discrete Element Method (DEM) investigation (Discontinuous Deformation Analysis, (DDA) [8, 9, 15], and including reducing friction effects and resulting mixed collapse modes [10, 11]; analytical-numerical, by a Complementarity Problem/Mathematical Programming formulation truly accounting for finite friction implications [12,13,14,15]; other studies on mainstream arch derivation and masonry arch mechanics after Heyman, specifically considering the issue of least thickness [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34]; Milankovitch formulation, as a cornerstone of thrust-line analysis [35,36,37]; modern and recent thrust-line-like analyses, particularly in view of form optimization [38,39,40,41,42,43,44].

A further wider literature framing within the vast subject of the analysis of masonry arches has already been provided in [7], and therein quoted references, with the integration of the above-quoted updates. Thus, the here-mentioned literature setting does not seek an intent of completeness, just that of focusing on self-coherent formulations that may consider further related, though different and complementary features, such as: different geometrical shapes of the masonry arch [22, 25, 30,31,32, 34, 39, 41], issues of “stereotomy” (i.e., the shape of the cut of the masonry blocks; here, a continuous circular arch with just potential rupture radial joints is considered) [6, 28, 33], dedicated use of “thrust-line analysis” or graphical-analytical methods, and so on, though all sharing the common characteristics of seeking and setting an analytical method of analysis, or being linked to the optimization target of acquiring the least-thickness condition. Numerical methods, which are massively employed in the analysis of masonry structures, including for masonry arches (vaults and domes), shall also be pursued and mentioned (see, for instance, about the adoption of the Discrete Element Method for the evaluation of the least-thickness condition [8, 9, 15, 30, 31], and therein quoted references), but these are here intentionally taken out from an explicit discussion. The study also focusses on the statics of masonry arches, not on the dynamics, which as well constitutes an important, though separate, topic, typically tied to numerical methods of analysis, and just for the self-supporting static condition under self-weight, though other external conditions or loadings may apply and be considered (e.g., lateral, as it may be linked to dynamical, seismic excitation conditions, for instance) [21, 24, 26, 32].

The present manuscript attempts to synthetically outline the main steps of the analytical derivations, to deliver the sought closed-form explicit representations and to present the illustration of the resulting outcomes, which are pursued by analytical plots. Section 2 gathers and develops the source governing equations and the relevant explicit expressions. Section 3 derives the stationary condition for the arch of maximum (intrinsic non-dimensional) horizontal thrust. Section 4 deepens the analysis of the foremost wide angular inner-hinge position from the crown, in the arising purely rotational collapse mode. Conclusions finally attempt to shortly recapitulate the main achievements of the study, as already above described and outlined, in itemized form, and draw down a few perspectives of possible further research investigation.

2 Main structure and analytical derivation

2.1 Governing equations

The classical definition of the main non-dimensional circular masonry arch quantities may be deciphered from the illustrations depicted in Figs. 1 and 2.

Specifically, for thickness-to-radius ratio \(\eta \) and non-dimensional horizontal thrust h, the following definitions may classically be stated:

$$\begin{aligned} \eta =\displaystyle \frac{t}{r}; \quad h=\displaystyle \frac{H}{w r}\quad (w=\gamma \,t\,d); \end{aligned}$$
(1)

where r, t and d are the radius, thickness and out-of-plane depth of the circular masonry arch, H is the horizontal thrust and \(\gamma \) and w are the specific weights per unit volume and per unit length of geometrical centreline of the arch. Together with angular inner-hinge position \(\beta \) from the crown (\(0\le \beta \le \alpha \), \(\alpha \) half-opening angle of the symmetric circular arch), they form characteristic solution triplet (\(\beta \), \(\eta \), h) of least-thickness “Couplet-Heyman problem”, and the attached determination of the associated purely rotational collapse mode (warranted by high friction, up to prevent any sliding), at variable \(\alpha \).

Fig. 1
figure 1

Symmetric circular masonry arch under self-weight (\(\gamma \) specific weight per unit volume) with all involved characteristic parameters (d out-of-plane depth)

Full size image
Fig. 2
figure 2

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

Full size image

The underlying governing equations linked to the analysis of the “Couplet-Heyman problem” of least-thickness optimization, in sought solution characteristic triplet (\(\beta \), \(\eta \), h), may be expressed, in terms of variable h, with appropriate on/off control \(\delta \) flags, up to shift among possible different treatments, and attached solutions [7, 9, 11, 15]:

$$\begin{aligned} \left\{ \begin{aligned} h&=h_1=h_L=\displaystyle \frac{(2-\eta ) \beta S-2 (1-C) \, {fac}_M}{2+\eta -(2-\eta )C };\\ h&=h_2=A-\displaystyle \frac{2}{2+\eta } \, {fac}_M;\\ h&=h_e=\!\!\!\!\!\underbrace{\displaystyle \frac{f-S}{S}}_{\displaystyle h_H=\beta C/S}\!\!\!\!\!+\delta _{{CCR}}\left( 1-\displaystyle \frac{2}{2-\eta } \, {fac}_M\right) ; \end{aligned} \right. \end{aligned}$$
(2)

where Milankovitch factor \({fac}_M\) of correct self-weight distribution (true location of the centres of gravity of each theoretical chunk of the continuous arch) [35,36,37] is

$$\begin{aligned} {fac}_{M}=\left( 1+\delta _{M}\eta ^2/12\right) ; \end{aligned}$$
(3)

as ruled by on/off control flag \(\delta _M\) (off, \(\delta _M=0\), then \({fac}_M=1\), approximate Heyman-like self-weight distribution along geometrical centreline; on, \(\delta _M=1\), then \({fac}_M=1+\eta ^2/12\), true Milankovitch-like self-weight distribution). Also, the following functions are introduced and adopted, for the convenience of the sought analytical representation [7]:

$$\begin{aligned} \begin{aligned} f&=(\beta S)'=S+\beta C; \quad g=f C+\beta S^2=\beta +S C;\\ S&=\sin \beta ; \quad C=\cos \beta ; \quad \text {with link}\quad \beta S^2+C f-g=0; \end{aligned} \end{aligned}$$
(4)

where symbol \('\) denotes a first-order differentiation (with respect to \(\beta \)) and, moreover, variable

$$\begin{aligned} A=\alpha \,\cot (\alpha /2) \end{aligned}$$
(5)

turns out the one that truly brings in the dependence on half-opening angle \(\alpha \) of the symmetric circular arch. Functions f and g, with mutual relation \(\beta S^2 + C f - g=0\), turn out useful to express the solution in compact form. Notice that \(g \ge 0\) for \(\beta \ge 0\); also, \(g - f=(1 - C)(\beta - S) \ge 0\) and \(f + g=(1 + C)(\beta + S)\ge 0\) for \(\beta \ge 0\). Functions f and g both vanish at \(\beta =0\) (see illustration in [15] and limit implications in [7, 15]).

The meaning of the stated governing equations goes as follows. Static equations (2)\(_{\mathrm{{a}}}\) and (2)\(_{\mathrm{{b}}}\) represent two equilibrium relations, respectively the rotational equilibrium of any upper portion AB of the half-arch (symmetry conditions apply), with respect to inner intrados hinge at haunch B (represented through thrust \(h_1\); namely, “lower thrust” \(h_L\) [15, 16]), and of whole half-arch AC, with respect to extrados hinge at shoulder C (represented through thrust \(h_2\)); the latter (only) sets in the dependence on opening angle \(\alpha \), through variable \(A=\alpha \,\cot (\alpha /2)\). Optimality equation (2)\(_{\mathrm{{c}}}\) states the tangency condition of the line of thrust (locus of pressure points) at haunch intrados B in the least-thickness condition and may also be derived from stationary condition \(h_1'=dh_1/d\beta =0\), leading to \(h=num_{h_1}'/den_{h_1}'=h_e\), namely \(h_e\) is obtained by the ratio between the derivatives of the numerator and the denominator in the expression of \(h_1\) in Eq. (2)\(_{\mathrm{{a}}}\) [7, 15]).

Indeed, the equation of the eccentricity of the line of thrust (locus of pressure points) from the geometrical centreline of the arch (or of non-dimensional eccentricity \({\hat{e}}=2e/t\)) may be obtained from a more general version of equilibrium relation (2)\(_{\mathrm{{a}}}\), of any upper portion of the half arch, by its rotational equilibrium with respect to the centre of pressure located at eccentricity e, at a general position \(\beta \) along the arch:

$$\begin{aligned} {\hat{e}}=\displaystyle \frac{e}{t/2}= \displaystyle \frac{2\beta S-2(1-C) {fac}_M-h \big ((2+\eta ) -2C\big )}{\eta (\beta S+h C)}; \end{aligned}$$
(6)

while, with \(h = h_2\) in Eq. (2)\(_{\mathrm{{c}}}\), as set to impose the passage of the line of thrust also at the shoulder extrados of the arch at \(\beta =\alpha \), one finally has:

$$\begin{aligned} {\hat{e}}=\displaystyle \frac{e}{t/2}= \displaystyle \frac{2 (2+\eta ) \big (\beta S-(1-C)A\big )-\eta \big ((2+\eta )A-2 C {fac}_M\big )}{\eta (2+\eta ) \big (\beta S-(1-C)A\big )+\eta \big ((2+\eta )A-2 C {fac}_M\big )}. \end{aligned}$$
(7)

Then, above third governing equation (2)\(_{\mathrm{{c}}}\), \(h=h_e\), actually comes from the stationary condition of eccentricity e (or \({{\hat{e}}}'(\beta )=0\)), at intrados inner-hinge position B, where, at the same time, \({\hat{e}}=1\) (\(e=t/2\)) [7], leading to equate the derivatives of the numerator and the denominator in Eq. (6), \(1=num_{{\hat{e}}}\,'/den_{{\hat{e}}}\,'\), which, solved with respect to h, precisely leads to condition \(h=h_e\).

Thus, the correct representation of such a tangency condition brings in a main difference to classical Heyman solution, which instead stated the tangency condition on the thrust force itself (\(h = h_H=\beta \cot \beta \) in Eq. (2)\(_{\mathrm{{c}}}\), and may be considered as leading to a (reasonable) solution approximation, at least until \(\eta \) (as \(\alpha \)) keeps small [see the expression of \({fac}_M\), quadratic in \(\eta \), in Eq. (3)]. This is ruled by on/off control flag \(\delta _{CCR}\) (off, \(\delta _{CCR}=0\), classical Heyman solution; on, \(\delta _{CCR}=1\), correct non-Heyman, CCR and Milankovitch solutions, considering the true tangency condition of the line of thrust at the intrados).

Then, the meaningful couples of the on/off control \(\delta \) flags, for the three underlying solutions that are considered in the equations and in the sequel are:

  • (\(\delta _{CCR}=0\), \(\delta _{M}=0\)) for Heyman solution;

  • (\(\delta _{CCR}=1\), \(\delta _{M}=0\)) for CCR solution;

  • (\(\delta _{CCR}=1\), \(\delta _{M}=1\)) for Milankovitch solution.

Least-thickness governing system (2) shall be solved to deliver characteristic solution triplet (\(\beta \), \(\eta \), h) at given half-opening angle \(\alpha \) of the arch (thus, at given \(A=\alpha \,\cot (\alpha /2\))), namely as (\(\beta (\alpha )\), \(\eta (\alpha )\), \(h(\alpha )\)), or (\(\beta (A)\), \(\eta (A)\), h(A)). However, ruling system (2) shows that it would actually be easier to conversely determine solution triplet (\(A(\beta )\), \(\eta (\beta )\), \(h(\beta )\)), corresponding to a least-thickness arch condition displaying the inner hinge at a given angular inner-hinge position \(\beta \). This constitutes a convenient way to turn around the analytical problem. Also, the use of variable A, which plays a linear role in Eq. (2)\(_{\mathrm{{b}}}\), instead as of \(\alpha \), turns out much convenient, for a full analytical treatment, leaving to a later stage the issue of then solving transcendental relation \(A=\alpha \,\cot (\alpha /2)\), to go back to the sought underlying dependence of three arch solution characteristics (\(\beta \), \(\eta \), h) at a given half-angle of embrace \(\alpha \).

2.2 Solution instances

Analytical solutions of governing system (2), in the just-defined above spirit, are now made explicit, for the three solution cases of Heyman [3], CCR [7] and Milankovitch [35,36,37].

2.2.1 Heyman solution (\(\delta _{{CCR}}=0\), \(\delta _M=0\)): linear algebraic problem

Classical Heyman solution [3] is readily obtained from governing system (2) once \(\delta _{CCR}=0\), \(\delta _M=0\), in terms of triplet \(A(\beta )\), \(\eta (\beta )\), \(h(\beta )\), since variables (A, \(\eta \), h) linearly enter the three equations (linear algebraic problem):

$$\begin{aligned} \left\{ \begin{aligned} A_H&=h_H+\displaystyle \frac{2}{2+\eta _H}=\beta \cot \beta \, \displaystyle \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 };\\ \eta _H&=2\,\displaystyle \frac{(\beta -\sin \beta )(1-\cos \beta )}{\beta (1+\cos \beta )};\\ h_H&=\beta \cot \beta . \end{aligned} \right. \end{aligned}$$
(8)

Indeed, for \(\delta _{CCR}=0\), \(h_H\) is directly obtained from Eq. (2)\(_{\mathrm{{c}}}\); then, substituted into Eq. (2)\(_{\mathrm{{a}}}\), which is solved with respect to \(\eta \), to get \(\eta _H\), and finally, by solving Eq. (2)\(_{\mathrm{{b}}}\) with respect to A, with \(\eta =\eta _H\), \(h=h_H\) inserted into it, to get \(A_H\). This “linear” feature of Heyman solution, extensively described in [7], even with different arrangements of the resulting equations, constitutes a beautiful characteristic property apt to simplify the analysis and still leads to reasonably approximate solutions until \(\beta \) keeps small (then \(\alpha \), or anyway for \(\alpha \) not nearing or even exceeding 90\(^{\circ }\), i.e., for under-complete arches). Resulting trends \(A_H(\beta ),~\eta _H(\beta ),~h_H(\beta )\) turn out to be single-valued and monotonic, which leads to an increasing divergence of the solution prediction at increasing \(\beta \), and eventually to the uncorrect prediction of angular inner-hinge position \(\beta \) in the collapse mechanism, as it will later be shown, in the confrontation of all the solution cases. End predictions of “external” features \(\eta _H\) and \(h_H\) may keep reasonable, at variable A or \(\alpha \), despite missing a correct evaluation of “internal” geometrical feature \(\beta _H\).

Next, for non-Heyman, CCR and Milankovitch solutions (\(\delta _{ CCR}=1\)), setting the correct tangency condition of the line of thrust at the inner hinge at the haunch intrados, an algebraic elimination from source governing system (2) of two out of three variables \(A,\,\eta ,\,h\) leads to an up-to-cubic algebraic problem in \(A(\beta ),\,\eta (\beta ),\,h(\beta )\), described by the following three polynomial equations with \(\delta _{M}\) flag:

$$\begin{aligned} \left\{ \begin{aligned}&{pol}_A=S^2\left( 3 g-\displaystyle \frac{3+\delta _M}{2} S\right) A^3 + 3 g S \left( g-f-\displaystyle \frac{1-\delta _M}{2} S\right) A^2-3 g^2\left( f-\displaystyle \frac{1-\delta _M}{2}S\right) A+\displaystyle \frac{3+\delta _M}{2} g^3=0;\\&{pol}_\eta =\delta _MS \,\eta ^3 + 3 (f+g)\,\eta ^2-12(g-S)\,\eta +12 (g-f)=0;\\&{pol}_h=6 S^2 \,h^3-3 S (3 f-g-2 S)\,h^2 -3(g-f)(f-2 S)\,h +\displaystyle \frac{3+\delta _M}{2} (g-f)^2=0. \end{aligned} \right. \end{aligned}$$
(9)

Notice that a full cubic degree of the governing polynomials is clearly ruled by a nonzero \(\delta _M\) flag, as specifically apparent from Eq. (9)\(_{\mathrm{{b}}}\), which, for \(\delta _M=0\), immediately degenerates into a quadratic equation, as it actually is, after simplification, for the other two equations, by eliminating redundant solutions (as coherent with a direct setting of a vanishing \(\delta _M\) control flag right on top of the process, namely directly on source governing equations (2)).

2.2.2 CCR solution (\(\delta _{CCR}=1\), \(\delta _M=0\)): quadratic algebraic problem

By setting \(\delta _M=0\) in above polynomial equations (9), after reduction, or directly from source equations (2), with also \(\delta _{CCR}=1\), by an algebraic elimination of two out of three variables \(A,\,\eta ,\,h\), the following quadratic algebraic problem in triplet \(A(\beta ),\,\eta (\beta ),\,h(\beta )\) is recovered:

$$\begin{aligned} \left\{ \begin{aligned} {pol}_A^{CCR}&=S(2g-S)\,A^2-2 f g \,A+g^2=0;\\ {pol}_\eta ^{CCR}&= (f+g)\,\eta ^2-4(g-S)\,\eta +4 (g-f)=0;\\ {pol}_h^{CCR}&=2 S \,h^2 -2(f-S)\,h+g-f=0. \end{aligned} \right. \end{aligned}$$
(10)

The arising solution can readily be expressed by the analytical solution of the second-order equations, as follows [7]:

$$\begin{aligned} \begin{aligned} A_{CCR}=g\;\displaystyle \frac{f \pm \sqrt{f^2-2 g S+ S^2}}{S(2g-S)};\\ \eta _{CCR}=2\;\displaystyle \frac{g-S{\mp }\sqrt{f^2-2 g S+S^2}}{f+g};\\ h_{CCR}=\displaystyle \frac{f-S\pm \sqrt{f^2-2 g S+S^2}}{2 S};\\ \end{aligned} \end{aligned}$$
(11)

where term \(f^2-2 g S+S^2\) under square root is proportional to the discriminant of the obtained quadratic equations and is positive in the range of \(\beta \) of interest, until it vanishes at \(\beta _{s\beta }^{CCR}=1.129085087576187\) (as further discussed in later Section 4), assuring real-valued branches in solution triplet \(A(\beta )\), \(\eta (\beta )\), \(h(\beta )\). Then, notice that, with respect to monotonic Heyman solution, each instance of the triplet becomes a two-valued function of angular inner-hinge position \(\beta \). The signs in front of the square root terms are sorted out either in triplet (+,−,+) or in triplet (−,+,−), i.e., with sorting \((A^+,\eta ^-, h^+)\) in the first branch and \((A^-,\eta ^+, h^-)\) in the second branch (different colours, for the two branches, will be adopted, in the illustrations that follow).

2.2.3 Milankovitch solution (\(\delta _{{CCR}}=1\), \(\delta _M=1\)): cubic algebraic problem

By setting \(\delta _M=1\) in polynomial equations (9), one directly gets the following cubic algebraic problem in triplet \(A(\beta ),\,\eta (\beta ),\,h(\beta )\):

$$\begin{aligned} \left\{ \begin{aligned} {pol}_A^M&=S^2(3 g-2 S) \,A^3 + 3 gS(g-f)\,A^2-3 f g^2 \,A+2 g^3=0;\\ {pol}_\eta ^M&= S \,\eta ^3 + 3 (f+g)\,\eta ^2-12(g-S)\,\eta +12 (g-f)=0;\\ {pol}_h^M&=6 S^2\,h^3-3 S (3 f-g-2 S)\,h^2 -3 (g-f) (f-2 S)\,h +2 (g-f)^2=0.\\ \end{aligned} \right. \end{aligned}$$
(12)

The solution of the cubic equations can be analytically developed in a classical (though little involved) way and represented as follows (two meaningful real-valued branches; same sign, either \(+\) or − in the two terms where ± appears, i.e., \(+\), \(+\) or −, −, or \({\mp }\) appears, i.e., −, − or \(+\), \(+\), in the given order):

$$\begin{aligned} \begin{aligned} A_M&=\displaystyle \frac{g}{S(3g-2S)}\Bigg (-(g-f)\\&\quad +\bigg (\displaystyle \frac{1\pm \mathrm {i} \sqrt{3}}{2^{2/3}}\bigg )\displaystyle \frac{ (f+g)^2-f (g +2\text { }S)}{\Big ((g-f)\big (2(f+g)^2+f (g-6 S)\big )+2S (3 g-2 S)^2 +(3g-2 S)\; {{SQRT}}\Big )^{1/3}}\\&\quad +\bigg (\displaystyle \frac{2^{2/3}}{1\pm \mathrm {i} \sqrt{3}}\bigg )\Big ((g-f)\big (2(f+g)^2+f (g-6 S)\big )+2S (3 g-2 S)^2 +(3 g-2 S) \; {{SQRT}}\Big )^{1/3}\Bigg ); \end{aligned}\nonumber \\ \end{aligned}$$
(13)
$$\begin{aligned} \begin{aligned} \eta _M&=\displaystyle \frac{1}{S}\Bigg (-(f+g)\\&\quad +\bigg (\displaystyle \frac{1{\mp } \mathrm {i} \sqrt{3}}{2}\bigg )\displaystyle \frac{(f+g)^2+4 S(g-S)}{\Big ((f+g)^3+6 gS (f+g)-12 f S^2+2S \; {{SQRT}}\Big )^{1/3}}\\&\quad +\bigg (\displaystyle \frac{2}{1{\mp } \mathrm {i} \sqrt{3}}\bigg )\Big ((f+g)^3+6 gS (f+g)-12 f S^2+2S \; {{SQRT}}\Big )^{1/3}\Bigg ); \end{aligned}\nonumber \\ \end{aligned}$$
(14)
$$\begin{aligned} \begin{aligned} h_M&=\displaystyle \frac{1}{6 S}\Bigg (3 f-g-2S\\&\quad +\bigg (\displaystyle \frac{1\pm \mathrm {i} \sqrt{3}}{2}\bigg )\displaystyle \frac{3 f^2+g^2-4 S (2 g-S)}{\Big (-g\left( 9 f^2-g^2\right) -2S\big (9 f (g-f)-12 g (g-S) -4 S^2\big )+3 (g-f)\; {{SQRT}}\Big )^{1/3}}\\&\quad +\bigg (\displaystyle \frac{2}{1\pm \mathrm {i} \sqrt{3}}\bigg )\Big (-g\left( 9 f^2-g^2\right) -2S\big (9 f (g-f)-12 g (g-S) -4 S^2\big )+3 (g-f)\; {{SQRT}}\Big )^{1/3}\Bigg ); \end{aligned}\nonumber \\ \end{aligned}$$
(15)

with

$$\begin{aligned} {SQRT}=\sqrt{9 S^2(g-f)^2 -2S (g-S) \left( 9 f^2-g^2-16 g S+8 S^2\right) -3 (f+g)^2 \left( f^2-2 g S+S^2\right) }; \end{aligned}$$
(16)

where i is the imaginary unit and the term under square root in Eq. (16) is proportional to minus the discriminant of the cubic equations and is negative in the range of \(\beta \) of interest, until it vanishes at \(\beta _{s\beta }^M=1.119864231443752\), assuring real-valued solution branches (distinguished by colours, in the subsequent illustrations). Again, solution triplet \(A(\beta )\), \(\eta (\beta )\), \(h(\beta )\) is represented by two-valued functions of angular inner-hinge position \(\beta \).

2.3 Graphical representations of the derived solutions

The three derived solutions can be analytically plotted from the solution branches represented in the above equations as a function of \(\beta \), for solution triplet \(A(\beta )\), \(\eta (\beta )\), \(h(\beta )\), as depicted in Fig. 3. Additionally, the least-thickness solution for arch triplet \(\beta , \eta , h\) is also analytically plotted in Fig. 4, from the above equations, as an analytical parametric plot as a function of A, thus for triplet \((A(\beta ),\beta )\), \((A(\beta ),\eta (\beta ))\), \((A(\beta ),h(\beta ))\), where, since \(A=\alpha \,\cot (\alpha /2)\), a half-angle of embrace \(\alpha \) increasing from 0 leads to an A decreasing from 2. Notice that no numerical evaluation is directly needed to produce the analytical plots displayed in Figs. 3 and 4, through the achieved closed-form explicit functions. Use of differentiating colours is made, for the underlying analytical branches of each characteristic function.

Fig. 3
figure 3

Analytical plot of triplet \(A(\beta )\), \(\eta (\beta )\), \(h(\beta )\), for the three solutions

Full size image
Fig. 4
figure 4

Analytical parametric plot of triplet \(\beta (A)\), \(\eta (A)\), h(A), for the three solutions

Full size image
Fig. 5
figure 5

Analytical-numerical plot of triplet \(\beta (\alpha )\), \(\eta (\alpha )\), \(h(\alpha )\), for the three solutions

Full size image

Figure 3 shows that the \(\beta \)-dependencies diverge for Heyman solution from some stage on, at increasing \(\beta \), since the described trends keep monotonic, showing that the prediction of this solution is not truly correct, in terms of “internal (geometrical) ingredient” \(\beta \). In fact, CCR and Milankovitch solutions, correctly accounting for the tangency condition of the line of thrust at the haunch intrados, reverse down after a maximum value of angular inner-hinge position \(\beta \) is reached. Conversely, CCR and Milankovitch solutions stay quite close to each other, with CCR a bit overshooting (CCR being in between Heyman and Milankovitch), in the post-stationary reversing paths, the predictions on non-dimensional horizontal thrust h keeping the closest.

Figure 4 confirms that a main divergence of Heyman solution with respect to CCR/Milankovitch solutions is recorded on the \(\beta \) dependencies on variable \(A=\alpha \,\cot (\alpha /2)\) (top plot), with similar overshooting trends. Amazingly, despite the prediction on \(\beta \) being increasingly wrong, all trends for \(\eta \) and h as a function of A keep rather close to each other, in resulting monotonic (with \(\eta \) and h, respectively, increasing and decreasing, at decreasing A, namely at increasing \(\alpha \)). Surprisingly, Heyman prediction on h turns out the closest to Milankovitch response.

Finally, given the source quest of inspecting the arch characteristics, in the least-thickness condition, at variable half-angle of embrace \(\alpha \), the further plot in Fig. 5 depicts the final dependencies of solution triplet \(\beta (\alpha ), \eta (\alpha ), h(\alpha )\). This step eventually requires a numerical evaluation (or a representation in implicit, rather than explicit form), since A and \(\alpha \) are linked by transcendental equation \(A=\alpha \,\cot (\alpha /2)\), and the values of \(\beta (\alpha ), \eta (\alpha ), h(\alpha )\) can be numerically evaluated at given values of \(\alpha \) (scored points), either from a direct numerical solution of governing system (2) or by numerical evaluations of the analytical solutions expressed for the cases of Heyman, in Eq. (8), of CCR, in Eq. (11), and of Milankovitch, in Eqs. (13)–(16).

3 Arch of maximum thrust: stationary condition on intrinsic non-dimensional horizontal thrust

3.1 Evaluation of intrinsic non-dimensional horizontal thrust

Now, main specific features of the least-thickness solutions are investigated. Particularly, the representation of the solutions in terms of “intrinsic non-dimensional horizontal thrust”, as introduced below, is investigated, with an analytical derivation of peculiar stationary properties. Indeed, referring just to given material (\(\gamma \)) and geometrical (d, r) parameters of the circular masonry arch, the following alternative definition of non-dimensional horizontal thrust may be stated:

$$\begin{aligned} \hat{h}={\displaystyle \frac{H}{\gamma d r^2}=\eta h}; \end{aligned}$$
(17)

meaning that the horizontal thrust that is present within the circular masonry arch is read in terms on known characteristic features of the arch, before setting the optimization problem in terms of the least-thickness condition on thickness t (or thickness-to-radius ratio \(\eta \)). Thus, despite appearing as \(\hat{h}=\eta h\), position (17) versus (1)\(_{\mathrm{{b}}}\), sets a definition explicitly independent on thickness t (or \(\eta \)), which constitutes one main subject of the sought least-thickness optimization process. Such a definition has been considered in [7] and shown to lead to interesting bell-shaped trends of \(\hat{h}\) (as opposed to monotonic trends of h), given relation \(\hat{h}=\eta h\), as it is going to be fully analytically investigated in the following.

Given the above definition of intrinsic non-dimensional horizontal thrust \(\hat{h}={\eta h}\), and considering the derived ruling equations and solution representations, concerning \(\eta \) and h, the following expressions of \(\hat{h}\) in the critical least-thickness condition can be derived, for the three solution instances.

For Heyman solution, linear algebraic problem

$$\begin{aligned} \hat{h}_H=\eta _H h_H=2\cot \beta \;\displaystyle \frac{(\beta -\sin \beta )(1-\cos \beta )}{1+\cos \beta }; \end{aligned}$$
(18)

directly obtained from Eq. (8), and complementing the relations within it.

For non-Heyman, CCR and Milankovitch solutions (\(\delta _{CCR}=1\)), tracing the above path, one may derive the following up-to-cubic polynomial equation in terms of on/off control flag \(\delta _M\), complementing those earlier given in Eq. (9):

$$\begin{aligned} {pol}_{\hat{h}}= & {} 2 \delta _MS^3\,\hat{h}^3 -3S (g-f) (f+g+2 \delta _MS)\,\hat{h}^2 \nonumber \\&+ 6 (g-f)^2 \big (f+(1+\delta _M) S\big )\,\hat{h}-2 (3+\delta _M) (g-f)^3=0; \end{aligned}$$
(19)

leading to the following specific representations of \(\hat{h}\).

For CCR solution (\(\varvec{\delta }_{\mathbf{M}}\varvec{=0}\) ), quadratic algebraic problem

$$\begin{aligned} {pol}_{\hat{h}}^{CCR}= S(f+g) \,\hat{h}^2-2 (g-f) (f+ S)\,\hat{h}+2 (g-f)^2; \end{aligned}$$
(20)

complementing those earlier given in Eq. (10), with solution out of the quadratic equation (plus/minus signs as in the above solution for \(h_{CCR}\)):

$$\begin{aligned} \hat{h}_{ CCR}=\eta _{CCR} h_{CCR}=(g-f)\;\displaystyle \frac{f+S\pm \sqrt{f^2-2 g S+S^2}}{S(f+g)}; \end{aligned}$$
(21)

complementing the representations provided in Eq. (11). It can be checked that expression (21) for \(\hat{h}_{CCR}\) may be obtained by the product of \(\eta _{CCR}\) and \(h_{CCR}\) in Eqs. (11)\(_{\mathrm{{b}}}\) and (11)\(_{\mathrm{{c}}}\).

For Milankovitch solution (\(\varvec{\delta }_{\mathbf{M}}\varvec{=1}\) ), cubic algebraic problem

$$\begin{aligned} {pol}_{\hat{h}}^M=2S^3\,\hat{h}^3 -3 S(g-f) (f+g+2 S) \,\hat{h}^2 + 6 (g-f)^2 (f+2 S)\,\hat{h}-8 (g-f)^3=0; \end{aligned}$$
(22)

complementing those earlier given in Eq. (12), with solution out of the cubic equation (minus/plus signs as in the above solution for \(\eta _M\)):

$$\begin{aligned} \begin{aligned} \hat{h}_M=\eta _{M} h_{M}&=\displaystyle \frac{g-f}{ 2 S^2} \Bigg (f+g+2 S\\&\quad -\bigg (\displaystyle \frac{1{\mp } \mathrm {i} \sqrt{3}}{2}\bigg ) \displaystyle \frac{(f+g)^2+4S (g-S)}{\Big ((f+g)^3+6 gS (f+g)-12 fS^2 +2 S \; {SQRT}\Big )^{1/3}}\\&\quad -\bigg (\displaystyle \frac{2}{1{\mp } \mathrm {i} \sqrt{3}}\bigg )\Big ((f+g)^3+6 gS (f+g)-12 fS^2 +2S \; {SQRT}\Big )^{1/3}\Bigg ); \end{aligned} \end{aligned}$$
(23)

complementing the representations provided in Eqs. (13)–(16). It can be checked that expression (23) for \(\hat{h}_{M}\) may be obtained by the product of \(\eta _{M}\) and \(h_{M}\) in Eqs. (14) and (15).

Now that the expressions of \(\hat{h}\) in the critical condition of least-thickness of the symmetric circular masonry arch have been derived in Eqs. (18), (21), (23), for the three solutions instances, and which can be coupled to Eqs. (8), (11) and (13)–(16), for the concomitant appreciation of the other arch characteristics, and thus, at variable half-opening angle of the arch, these can be analytically plotted and specific stationary conditions of interest that appear from the plots can then be analytically investigated.

Figure 6 indeed shows the dependencies of intrinsic non-dimensional horizontal thrust \(\hat{h}=\eta h\) on angular inner-hinge position \(\beta \), for the three solution instances. Besides for the already commented single-valued trend of \(\hat{h}\) at a monotonically increasing \(\beta \) for Heyman solution, as opposed to the double-valued \(\hat{h}(\beta )\) returning-to-zero branches of \(\hat{h}\) and \(\beta \) for CCR and Milankovitch solutions, the analytically plotted functions show that a characteristic stationary (maximum) value \(\hat{h}_{s\hat{h}}\) appears, which actually keeps rather similar for the three solutions (and strictly similar for CCR and Milankovitch solutions). In that sense, despite Heyman is failing in reproducing the correct \(\beta \) developments, still correctly estimates the values of \(\hat{h}\), and also peak value \(\hat{h}_{s\hat{h}}\). Notice that the horizontal tangency of curve \(\hat{h}(\beta )\) at the peak value of \(\hat{h}(\beta )\) for CCR and Milankovitch solutions is recorded along the first branch of solutions (21), \(+\) sign, and (23), − sign. Also, the two ±/\({\mp }\) branches of the solutions (distinguished by colours) are separated by the maximum reachable value of \(\beta \), where a vertical tangency of curve \(\hat{h}(\beta )\) is attained (which will be the subject of further analytical analysis in next Section 4). Thus, the maximum of \(\hat{h}\) occurs earlier, at increasing \(\beta \), than the reaching of the maximum value of \(\beta \), before that the reversing path of \(\beta \) (and \(\hat{h}\)) going back to zero has begun. Plots depicted in Fig. 6 for \(\hat{h}(\beta )\) complement and complete those illustrated in Fig. 3 and specifically provide a dual vision to those already traced for \(h(\beta )\) in Fig. 3c.

Even nearer trends of \(\hat{h}(A)\) are recorded for intrinsic non-dimensional horizontal thrust \(\hat{h}\) at variable \(A=\alpha \,\cot (\alpha /2)\), with eliminated angular inner-hinge position \(\beta \) as an “internal (geometrical) ingredient” of the formulation. As a rather interesting feature of newly defined non-dimensional horizontal thrust \(\hat{h}\), recorded trends \(\hat{h}(A)\) are “bell-shaped” at decreasing A, i.e., at increasing \(\alpha \) (instead of the monotonic, decreasing trends of h(A), in Fig. 4c). Obviously, this is induced by the two null values of \(\eta \) and h on the two sides of the interval of A, as depicted in the plots in Figs. 4b, c, given the fact that \(\hat{h}=\eta h\) corresponds to the product of \(\eta \) and h. Surprisingly, CCR solution is a bit overshooting, with respect to Milankovitch solution, and Heyman solution keeps rather close to Milankovitch solution in the post-peak branch at decreasing A (all solutions staying very close to each other in the pre-peak branch), with some little discrepancies in the evaluation of the maximum value of the horizontal thrust. Plots in Fig. 7 for \(\hat{h}(A)\) complement and complete those in Fig. 4, and specifically provide a dual vision to those already traced for h(A) in Fig. 4c.

Fig. 6
figure 6

Analytical plot \(\hat{h}(\beta )\) of intrinsic horizontal non-dimensional thrust \(\hat{h}=\eta h\), for the three solutions

Full size image
Fig. 7
figure 7

Analytical parametric plot \(\hat{h}(A)\) of intrinsic horizontal non-dimensional thrust \(\hat{h}=\eta h\), for the three solutions

Full size image
Fig. 8
figure 8

Analytical-numerical plot \(\hat{h}(\alpha )\) of intrinsic horizontal non-dimensional thrust \(\hat{h}=\eta h\), for the three solutions

Full size image

The same trends are confirmed by a further representation of the analytical-numerical dependencies of \(\hat{h}(\alpha )\), which again require a scored numerical evaluation linked just to the transcendental solution of A = \(\alpha \,\cot (\alpha /2)\). Here, the pre-peak branches, at increasing \(\alpha \), stay much close to each other, for sure for under-complete arches, for the complete semicircular arch and for not much over-complete arches with \(\alpha \) up to around 100\(^\circ \). Then, some undershooting for Heyman and overshooting for CCR solutions is recorded, with respect to reference Milankovitch solution. Again, for the further increasing \(\alpha \) values approaching those leading to a vanishing \(\hat{h}\) (and h), see extensive discussion in [7], curiously Heyman solution gets the closest to Milankovitch solution. The figure also clearly shows that the peak value of \(\hat{h}\) is reached for half-opening angles \(\alpha \) in a range at around 120\(^\circ \)–125\(^\circ \), with maximum values of \(\hat{h}\) that are in the order of 0.11 (detailed values will later be provided in table form at the end of the paper). Plots in Fig. 8 for \(\hat{h}(\alpha )\) complement and complete those in Fig. 5 and specifically provide a dual vision to those already traced for \(h(\alpha )\) in Fig. 5c.

3.2 Analysis of the stationary condition for the maximum thrust

Referring now to the plots in Fig. 6, where the extremum (maximum) condition of horizontal non-dimensional thrust \(\hat{h}\) is marked by the horizontal tangency condition at the peak of the \(\hat{h}=\hat{h}(\beta )\) curve, stationary condition \(\hat{h}'(\beta )=d\hat{h}(\beta )/d\beta =0\) can now be explicitly evaluated, given above analytical representations (18), (21), (23). This is straightforward, for Heyman solution. Alternatively, for non-Heyman, CCR and Milankovitch solutions, it may conveniently be derived as reported below, with closed-form symbolic representations.

For Heyman solution, directly from the first-order derivative of \(\hat{h}\) with respect to \(\beta \) in Eq. (18), setting its vanishing, for stating the stationary condition of function \(\hat{h}(\beta )\):

$$\begin{aligned} \hat{h}_H'=0 \quad \Rightarrow \quad (1-2 \cos \beta )\beta -\sin \beta \big (\sin ^2\beta -\cos \beta \big )=0; \end{aligned}$$
(24)

with numerical root, and corresponding intrinsic non-dimensional horizontal thrust:

$$\begin{aligned} \beta _{s\hat{h}}^H=1.3225601781836527 \quad \Rightarrow \quad \hat{h}_{s\hat{h}}^H=0.10842222971698252. \end{aligned}$$
(25)

For non-Heyman, CCR and Milankovitch solutions, all together, instead of directly working out the stationary condition on the first-order derivatives of the expressions (first branches) of \(\hat{h}(\beta )\) in Eqs. (21) and (23), the stationary analysis may conveniently be carried out as follows.

One may notice that polynomial function \({pol}_{\hat{h}}\) in horizontal non-dimensional thrust \(\hat{h}\) in Eq. (19), in solution instance \(\hat{h}(\beta )\) such that \({pol}_{\hat{h}}=0\), ultimately becomes a function of angular inner-hinge position \(\beta \) in identity \({pol}_{\hat{h}}\big (\beta ,\hat{h}(\beta )\big )=0\). Then, total differentiation of such an identity, with respect to \(\beta \), leads to:

$$\begin{aligned} \displaystyle \frac{d {pol}_{\hat{h}}}{d \beta }=\frac{\partial {pol}_{\hat{h}}}{\partial \beta }+\frac{\partial {pol}_{\hat{h}}}{\partial \hat{h}} \; \frac{\partial \hat{h}}{\partial \beta }=0. \end{aligned}$$
(26)

Since at the horizontal tangency of the peak of \(\hat{h}\), stationary condition \(\hat{h}(\beta )'\,{=}\,d\hat{h}/d\beta \,{=}\,\partial \hat{h}/\partial \beta =0\) has to be set, from above Eq. (26) the second additive term shall vanish, and one is then lead to simultaneously set the following system of up-to-cubic polynomial equations in \(\hat{h}\):

$$\begin{aligned} \left\{ \begin{aligned} {pol}_{\hat{h}}&=0;\\ \frac{\partial {pol}_{\hat{h}}}{\partial \beta }&=0. \end{aligned} \right. \end{aligned}$$
(27)

Then, to work out such a system of two polynomial equations in \(\hat{h}\), the up-to-cubic polynomial from Eq. (19) may compactly be re-written as follows:

$$\begin{aligned} {pol}_{\hat{h}}=a\,\hat{h}^3+ b\,\hat{h}^2+c\,\hat{h}+d=0; \end{aligned}$$
(28)

and its partial derivative with respect to \(\beta \) then derived as:

$$\begin{aligned} \displaystyle \frac{\partial {pol}_{\hat{h}}}{\partial \beta }=a'\,\hat{h}^3+b'\,\hat{h}^2+c'\,\hat{h}+d'=0; \end{aligned}$$
(29)

where, from above Eq. (19), polynomial coefficients a, b, c, d are functions of \(\beta \):

$$\begin{aligned} \left\{ \begin{aligned} a&=2 \delta _MS^3;\\ b&=-3 S(g-f) (f+g+2 \delta _MS);\\ c&=6 (g-f)^2 \big (f+(1+\delta _M) S\big );\\ d&=-2 (3+\delta _M) (g-f)^3; \end{aligned} \right. \end{aligned}$$
(30)

and \(a',b',c',d'\) are their first-order derivatives with respect to \(\beta \), whose expressions can readily by evaluated.

Now, by eliminating variable \(\hat{h}\) from two Eqs. (28) and (29), of system (27), the following implicit relation in \(\beta \) for the stationary condition of \(\hat{h}\) can be obtained:

$$\begin{aligned} \begin{aligned}&(a d'-a' d)^3\\&\quad +(b d'-b' d)^2 (a b'-a' b) +(a c'-a' c)^2 (c d'-c' d)-(a d' -a' d)^2(b c' -b' c)\\&\quad -(a b'-a' b) (c d'-c' d)\big ((b c' -b' c)+3(a d' -a' d)\big )=0; \end{aligned} \end{aligned}$$
(31)

leading to:

$$\begin{aligned} \begin{aligned}&3(3+4 \cos \beta +\cos 2\beta ) \Bigg (\beta ^3(1-2 \cos \beta )\displaystyle \frac{\cot \beta /2}{\sin ^2\beta /2}+ 3 \beta ^2-4\beta \sin \beta -\sin ^2\beta \Bigg )\\&+\quad \delta _M\;\Bigg (2 \beta ^3 (6 -11\cos \beta + \cos 2\beta )\displaystyle \frac{\cot \beta /2}{\sin ^2\beta /2}+3 \beta ^2 (13+10 \cos \beta +9 \cos 2\beta )\\&\quad \qquad \qquad \quad -6\beta \sin \beta (15-4 \cos \beta +\cos 2\beta )+2 \sin ^2\beta (1-\cos \beta )(26+\cos \beta )\Bigg )\\&+\quad \delta _M^{2} \; 18(1-\cos \beta ) (\beta -\sin \beta )^2=0. \end{aligned} \end{aligned}$$
(32)

For CCR solution (\(\varvec{\delta }_\mathbf{M}\varvec{=0}\) ), Eq. (32) reduces to:

$$\begin{aligned} \beta ^3(1-2 \cos \beta )\displaystyle \frac{\cot \beta /2}{\sin ^2\beta /2}+ 3 \beta ^2-4 \beta \sin \beta -\sin ^2\beta =0; \end{aligned}$$
(33)

with numerical root, and corresponding intrinsic non-dimensional horizontal thrust:

$$\begin{aligned} \beta _{s\hat{h}}^{CCR}=1.123938439075214 \quad \Rightarrow \quad \hat{h}_{s\hat{h}}^{CCR}=0.11479135769565468. \end{aligned}$$
(34)

For Milankovitch solution (\(\varvec{\delta }_\mathbf{M}\varvec{=1}\) ), Eq. (32) comes to:

$$\begin{aligned} \begin{aligned}&\beta ^3 (19+10 \cos \beta +3 \cos 2 \beta )(1-2 \cos \beta ) \displaystyle \frac{\cot \beta /2}{\sin ^2\beta /2}+ 12 \beta ^2 (7+4 \cos \beta +3 \cos 2 \beta )\\&\quad -6\beta \sin \beta (27-2 \cos \beta +3 \cos 2 \beta )+4 \sin ^2\beta (15-20 \cos \beta -\cos 2 \beta )=0; \end{aligned} \end{aligned}$$
(35)

with numerical root, and corresponding intrinsic non-dimensional horizontal thrust:

$$\begin{aligned} \beta _{s\hat{h}}^M=1.1154928768047523\quad \Rightarrow \quad \hat{h}_{s\hat{h}}^M=0.11230191519065938. \end{aligned}$$
(36)

Note

Notice that, for CCR solution, a simpler quadratic problem, in a similar way, may directly lead to:

$$\begin{aligned} \left\{ \begin{aligned} {pol}_{\hat{h}2}&=a\,\hat{h}^2+b\,\hat{h}+c=0;\\ \displaystyle \frac{\partial {pol}_{\hat{h}2}}{\partial \beta }&=a'\,\hat{h}^2+b'\,\hat{h}+c'=0; \end{aligned} \right. \end{aligned}$$
(37)

where

$$\begin{aligned} \left\{ \begin{aligned} a&=S(f+g);\\ b&=-2(g-f)(f+S);\\ c&=2(g-f)^2; \end{aligned} \right. \end{aligned}$$
(38)

leading to:

$$\begin{aligned} (a c'-a' c)^2-(a b'-a' b)(b c'-b' c)=0; \end{aligned}$$
(39)

namely:

$$\begin{aligned} \beta ^3\left( 3-\cot ^2\beta /2 \right) \cot \beta /2 + 3 \beta ^2 - 4\beta \sin \beta -\sin ^2\beta =0; \end{aligned}$$
(40)

which is equivalent to above Eq. (33), with same numerical solution (34).

4 Foremost wide angular inner-hinge position: stationary condition on \(\beta \)

Referring again to Fig. 6, just for non-Heyman, CCR and Milankovitch solutions, setting the condition of foremost wide angular inner-hinge position from the crown may be handled by imposing the vertical tangency condition of returning curve \(\beta =\beta (\hat{h})\), namely of its stationary condition with respect to \(\hat{h}\). This can be worked out by the following treatment, which correspondingly leads to the vanishing of the discriminant of all the common polynomial source equations (namely to the null values of the terms under square root in Eqs. (11) and (21), and in Eqs. (13)–(16) and (23), i.e., that setting to zero term SQRT in Eq. (16)).

In this case, one may interpret that polynomial function \({pol}_{\hat{h}}\) in horizontal non-dimensional thrust \(\hat{h}\) in Eq. (19), in solution instance \(\beta (\hat{h})\) such that \({pol}_{\hat{h}}=0\), ultimately becomes a function of \(\hat{h}\) in identity \({pol}_{\hat{h}}\big (\beta (\hat{h}),\hat{h}\big )=0\). Then, total differentiation of such an identity, with respect to \(\hat{h}\), leads to:

$$\begin{aligned} \displaystyle \frac{d{pol}_{\hat{h}}}{d \hat{h}}=\frac{\partial {pol}_{\hat{h}}}{\partial \hat{h}}+\frac{\partial {pol}_{\hat{h}}}{\partial \beta } \; \frac{\partial \beta }{\partial \hat{h}}=0. \end{aligned}$$
(41)

Since at the vertical tangency of the maximum of \(\beta \), stationary condition \(d\beta /d\hat{h}=\partial \beta /\partial \hat{h}=0\) has to be set, from above Eq. (41) the second additive term shall vanish and one is then lead to simultaneously set the following system of one up-to-cubic polynomial equation and one up-to-quadratic polynomial equation in \(\hat{h}\):

$$\begin{aligned} \left\{ \begin{aligned} {pol}_{\hat{h}}&=0;\\ \frac{\partial {pol}_{\hat{h}}}{\partial \hat{h}}&=0. \end{aligned} \right. \end{aligned}$$
(42)

Then, since \({pol}_{\hat{h}}\) can be taken from up-to-cubic representation (28), one easily gets, from the partial derivative with respect to \(\hat{h}\):

$$\begin{aligned} \left\{ \begin{aligned} {pol}_{\hat{h}}&=a\,\hat{h}^3+b\,\hat{h}^2+c\,\hat{h}+d=0;\\ \displaystyle \frac{\partial {pol}_{\hat{h}}}{\partial \hat{h}}&=3a\,\hat{h}^2+2b\,\hat{h}+c=0. \end{aligned} \right. \end{aligned}$$
(43)

By eliminating \(\hat{h}\) from the two polynomial equations one has:

$$\begin{aligned} c^2 \left( b^2-4 a c\right) -d (4 b^3+27 a^2 d-18 a b c)=0; \end{aligned}$$
(44)

which coincides with the vanishing condition of cubic discriminant

$$\begin{aligned} \Delta =b^2 c^2-4 a c^3-4 b^3 d-27 a^2 d^2+18 a b c d=0; \end{aligned}$$
(45)

leading to:

$$\begin{aligned} -9S^2(g-f)^2 \, \delta _M^2 +2 S(g-S) \left( 9 f^2-g^2-16 g S+8 S^2\right) \, \delta _M+3 (f+g)^2 \left( f^2-2g S +S^2\right) =0; \end{aligned}$$
(46)

and, explicitly, in \(\beta \):

$$\begin{aligned} \begin{aligned}&-9 \delta _M^2 (1-\cos \beta )^2 (\beta -\sin \beta )^2 \sin ^2\beta \\&\quad +2 \delta _M\sin \beta \Big (\beta -(1-\cos \beta ) \sin \beta \Big ) \Big (8 \sin ^2\beta +9 (\sin \beta +\beta \cos \beta )^2-(\beta +\sin \beta \cos \beta )^2-8 \sin \beta (2 \beta +\sin 2\beta )\Big )\\&\quad +3 (1+\cos \beta )^2 (\beta +\sin \beta )^2 \Big (\sin ^2\beta +(\sin \beta +\beta \cos \beta )^2-2 \sin \beta (\beta +\sin \beta \cos \beta )\Big )=0. \end{aligned} \end{aligned}$$
(47)

For CCR solution (\(\varvec{\delta }_\mathbf{M}\varvec{=0}\))

Polynomial equation (46) comes to correspond to set the vanishing of term \(f^2-2 g S+S^2\) under square root in Eqs. (11) and (21), namely:

$$\begin{aligned} f^2-2 g S+S^2=(\sin \beta +\beta \cos \beta )^2-2(\beta + \sin \beta \cos \beta ) \sin \beta +\sin ^2\beta =0; \end{aligned}$$
(48)

which indeed corresponds to the vanishing condition of the discriminant of earlier quadratic polynomial (37)\(_{\mathrm{{a}}}\) (a, b, c from Eq. (38)):

$$\begin{aligned} b^2-4 a c=0 \quad \Rightarrow \quad f^2-2 g S+S^2=0; \end{aligned}$$
(49)

readily coming from system

$$\begin{aligned} \left\{ \begin{aligned} {pol}_{\hat{h}}&=a \, \hat{h}^2+b\, \hat{h}+c=0;\\ \displaystyle \frac{\partial {pol}_{\hat{h}}}{\partial \hat{h}}&=2a\, \hat{h}+b=0; \end{aligned} \right. \end{aligned}$$
(50)

by solving for \(\hat{h}=-b/(2a)\) from the second equation and inserting it into the first equation, leading to \(b^2-4 a c=0\). All this eventually leads to the following numerical solution of foremost wide angular inner-hinge position:

$$\begin{aligned} \beta _{s\beta }^{CCR}=1.129085087576187\;\text{ rad }= 64.69181022927447^\circ . \end{aligned}$$
(51)

For Milankovitch solution (\(\varvec{\delta }_\mathbf{M}\varvec{=1}\) )

Polynomial equation (46) comes to correspond to set the vanishing of the (opposite of the) term under square root in SQRT, Eq. (16), namely

$$\begin{aligned} 9 S^2 (g-f)^2 -2S (g-S) \left( 9 f^2-g^2-16 g S+8 S^2\right) -3 (f+g)^2 \left( f^2-2 g S+S^2\right) =0. \end{aligned}$$
(52)

This explicitly leads to:

$$\begin{aligned} \begin{aligned}&3 \beta ^4\cos ^2\beta (1+\cos \beta )^2\\&\quad -2 \beta ^3 \bigg (4+3 \cos \beta \Big (1-\cos \beta \big (5+\cos \beta (3+\cos \beta )\big )\Big )\bigg )\sin \beta \\&\quad -3 \beta ^2 \bigg (15-\cos \beta \Big (14-\cos \beta \big (6-\cos \beta (10+\cos \beta )\big )\Big )\bigg ) \sin ^2\beta \\&\quad +6 \beta (1-\cos \beta ) \big (15-(6-\cos \beta ) \cos \beta \big ) \sin ^3\beta \\&\quad -(1-\cos \beta ) \big (37-\cos \beta (53+8 \cos \beta )\big ) \sin ^4\beta =0; \end{aligned} \end{aligned}$$
(53)

or (equivalent form):

$$\begin{aligned} \begin{aligned}&3 \beta ^4\cos ^2\beta (1+\cos \beta )^2\\&\quad -2 \beta ^3 \left( 4+3 \cos \beta -15 \cos ^2\beta -9 \cos ^3\beta -3 \cos ^4\beta \right) \sin \beta \\&\quad -3 \beta ^2 \left( 15-14 \cos \beta +6 \cos ^2\beta -10 \cos ^3\beta -\cos ^4\beta \right) \sin ^2\beta \\&\quad +6 \beta \left( 15-21 \cos \beta +7 \cos ^2\beta -\cos ^3\beta \right) \sin ^3\beta \\&\quad -(37-90 \cos \beta +45 \cos ^2\beta +8 \cos ^3\beta ) \sin ^4\beta =0; \end{aligned} \end{aligned}$$
(54)

with numerical solution of foremost wide angular inner-hinge position:

$$\begin{aligned} \beta _{s\beta }^{M}=1.119864231443752\;\text{ rad }= 64.16349408938858^\circ . \end{aligned}$$
(55)

The characteristic extremal arch features in the least-thickness condition out of the present analytical analysis are finally gathered in Table 1, for the three solution instances of Heyman, CCR and Milankovitch. Notice that the following summarizing peculiar features may be itemized:

  • the \(\alpha \) for the stationary value of \(\hat{h}\), i.e., \(\alpha _{s\hat{h}}\), is in the order of 120\(^\circ \) \(-\)123\(^\circ \), while the \(\alpha \) for the stationary value of \(\beta \), i.e., \(\alpha _{s\beta }\), is in the order of 128\(^\circ \) \(-\)125\(^\circ \). First, stationarity on \(\hat{h}\) is reached, at increasing \(\alpha \), then stationarity on \(\beta \) is achieved (see also next Fig. 9).

  • Though being rather close to each other, the two stationary conditions differ, while one may anyway conclude that these extremal features are recorded for over-complete symmetric circular masonry arches with half-opening angles in the order of 120\(^\circ \)–128\(^\circ \).

  • The widest angular inner-hinge position from the crown turns out to be in the order of 64\(^\circ \) \(-\)65\(^\circ \).

  • The maximum value of \(\hat{h}\) is in the order of 0.11, meaning that the maximum value of horizontal thrust that a symmetric circular masonry arch with assigned specific weight per unit volume \(\gamma \), out-of-plane depth d and radius r, as given material and geometrical properties, may be able to transfer to the abutments, in the critical least-thickness condition, is in the order of \(H^{max}=0.11 \gamma d r^2\).

  • A rather similar value of horizontal thrust is also recorded for the opening angles leading to the foremost wide angular inner-hinge position.

  • In conclusion, both stationary conditions basically lead to a circular masonry arch pushing at around \(H^{max}=0.11\gamma d r^2\), with an angular inner-hinge position from the crown at about 64\(^\circ \).

Table 1 Characteristic extremal features (\(\alpha \) half-opening angle, \(A=\alpha \,\cot \)(\(\alpha /2\)), \(\beta \) angular inner-hinge position from the crown, \(\eta =t/r\) thickness-to-radius ratio, \(h=H/\)[(\(\gamma t d\))\(\,r\)] non-dimensional horizontal thrust, \(\hat{h}=H/\)(\(\gamma d r^2)=\eta h\) intrinsic non-dimensional horizontal thrust), for a symmetric circular masonry arch of maximum intrinsic non-dimensional horizontal thrust (stationarity of \(\hat{h}\), label “\(s\hat{h}\)”, first block lines) and of foremost wide angular inner-hinge position (stationarity of \(\beta \), label “\(s\beta \)”, second block lines), for the three considered solutions. Values are academically reported with 6 significative digits, to help in inspecting the differences between the three derived solutions
Full size table
Fig. 9
figure 9

Analytical sketch of peculiar least-thickness half-arches with line of thrust, for CCR solution

Full size image

Further, Fig. 9 consistently shows the analytical geometrical plots of the arch in the critical least-thickness condition (for CCR solution), which strictly include the line of thrust (being tangent to the intrados at the inner-hinge position), from the eccentricity expressed in Eq. (7), with values of half-opening angle \(\alpha \) allowing to make a fork around the two stationary values on \(\hat{h}\) and on \(\beta \), which are, respectively, represented on the top-right and bottom-left analytical sketches.

Independent numerical validation results (and further graphical illustration) by a Complementarity Problem/Mathematical Programming formulation [13, 14] are separately reported in [15], with an achieved full consistency in the various outcomes, on all the characteristic features of the least-thickness condition, namely kinematical (angular inner-hinge position \(\beta \), defining the purely rotational collapse mode), geometrical (thickness-to-radius ratio \(\eta \)), statical (non-dimensional horizontal thrust h or \(\hat{h}\)).

5 Conclusions

The present study has allowed to derive explicit closed-form analytical representations of the solution characteristics for (symmetric) circular masonry arches in the least-thickness condition (relying on infinite friction), developed all together, for the free solution instances of Heyman, CCR and Milankovitch. For the latter, explicit mathematical expressions have newly been provided, by the analytical solution of the underlying cubic equations.

Then, interesting aspects of the features of the derived solutions have been inspected, specifically in view of explicitly characterizing two extremal properties, namely the condition of maximum intrinsic non-dimensional horizontal thrust and of foremost wide angular inner-hinge position from the crown. The two instances occur for half-opening angles at around 120\(^\circ \)–128\(^\circ \), thus in the realm of over-complete (horseshoe) arches, as peculiar in Moorish architecture, though slightly differ from each other. Both extremal conditions reveal a non-monotonic trend for the intrinsic non-dimensional horizontal thrust (while thickness-to-radius ratio and regular non-dimensional horizontal thrust keep monotonic, one increasing, the other decreasing, at increasing arch opening) and for the angular inner-hinge position. The former shows that the horizontal thrust within the arch in the critical least-thickness condition goes to a maximum at increasing opening angle; the latter also that, at some stage, the inner hinge pulls back towards the crown, after an initial increase, at increasing arch opening. These mechanical features, appearing for little-opened horseshoe arches constitute intriguing characteristics of the studied mechanical system, concerning the mechanics (statics) of circular masonry arches, for sure from an academical interpretation point of view but also in terms of practical interest, in the engineering understanding of adopted (historical) architectural forms.

One general remark may be that, on one side, even a rather simple mechanical system shall require a quite considerable effort to develop full closed-form explicit solutions by a complete analytical approach, thus constituting, still, a real challenge towards achieving a consistent Mechanics modelization; on the other side, it displays such rich and non-obvious engineering features, in exploring a full range of arch openings, then going to entirely appreciate the mechanical matter, also for over-complete arches, with analytical results that anyway constitute a fundamental reference for the understanding of the masonry arch behaviour and in terms of safe dimensioning.

Indeed, concerning the statement of achieved academical and practical value of the present investigation, the following resuming remarks may be outlined. The analysis first seeks a solution (by a complete analytical treatment) in all possible ranges of the variables ruling the mechanical problem. Thus, in principle, half-angles of embrace from 0\(^\circ \) to 180\(^\circ \) are taken into account (then, coming to include horseshoe arches). This is just from a pure Mechanics stand point, to avoid losing any possible solution and aiming at contributing to the general understanding of the Mechanics of masonry arches. The value of the manuscript shall mainly lay on an academic ground but the attached practical value shall not be limited. It is just meant that “optimum” arches in the studied least-thickness condition are obviously far to be adopted in real human constructions, due to the need of achieving proper safety margins from the studied collapse states. Thus, supercritical arch thicknesses are expected and indeed adopted, in practice (also due to additional loading that is usually present). The analytical representations that are derived within the paper shall be unprecedented (as far as the writers know). They constitute the main contribution and can be employed to enquire and derive basic conditions of masonry arch Mechanics, through explicit formulas, which could be adopted for first dimensioning or for benchmark reference, towards comparison of built shapes. In the end, they may also be used for teaching purposes.

Finally, to truly appreciate the type of scientific content and technical derivations herein presented, and aim and value of the present endeavour, the interested reader may attempt trying to independently replicate the mathematical steps that have allowed to achieve the present derivations and aspired closed-form explicit representations. All necessary elements shall have been given in the paper, in a self-contained manner, to make this challenge possible, still as an unknown way to be traced or, at least, as an “a posteriori” verification of the mathematical forms that have been pursued, and now wrapped up, to be shared with the community by the enclosed description.

Perspectives in terms of the mechanical study may further consider different geometrical arrangements, for instance for non-symmetric arches, which shall set in an additional geometric parameter into the picture, with implications that may even go much beyond that, due to expected changes in the collapse mechanism to be recorded, and linked to possible different considered form shapes of the masonry arch. Moreover, the role of finite friction shall systematically be enquired, as a possible weakening argument, to handle also collapse instances that may include sliding [10, 11, 13, 14].