Abstract
The paper deals with the variational setting of the optimal archgrid construction. The archgrids, discovered by William Prager and George Rozvany in 1970s, are viewed here as tensionfree and bendingfree, uniformly stressed gridshells forming vaults unevenly supported along the closed contour of the basis domain. The optimal archgrids are characterized by the least volume. The optimization problem of volume minimization is reduced to the pair of two auxiliary mutually dual problems, having mathematical structure similar to that known from the theory of optimal layout: the integrand of the auxiliary minimization problem is of linear growth, while the auxiliary maximization problem involves test functions subjected to meansquare slope conditions. The noted features of the variational setting governs the main properties of the archgrid shapes: they are vaults over a subregion of the basis domain being the effective domain of the minimizer of the auxiliary problem. Thus, the method is capable of cutting out the material domain from the design domain; this process is built in within the theory. Moreover, the present paper puts forward new methods of numerical construction of optimal archgrids and discusses their applicability ranges.
Introduction
Funiculars are planar frameworks which, albeit subject to a transverse and transmissible load, do not undergo flexure, i.e. they remain bendingfree and, consequently, are not subject to transverse shear forces, while the axial stress resultants are characterized by a fixed sign: all the bars are in tension (or in compression). The load is called transmissible, if it follows the design, yet keeping the given (usually vertical) direction and keeping the value of its intensity q(x) measured per unit length of the line orthogonal to the direction of the load [q] = [N/m], cf. Fuchs and Moses (2000). Of the funiculars being fully and uniformly stressed (up to a given limit (equal, say, − σ_{C}, σ_{C} being the permissible stress in compression)) one can find one of the least volume. Its rise cannot be too high but also cannot be too small, since big values of thrusts (i.e. horizontal reactions at the pin supports) lead to an inevitable increase of areas of the crosssections, resulting in the increase of the volume. Thus, the rise should be appropriately chosen. Let z = z(x) represent the shape of the funicular pin supported at A and B at the levels z(0) = h_{A}, z(l) = h_{B}.
The optimal rise of the leastvolume funicular is determined by the condition
discovered by Rozvany and Prager (1979) and Rozvany and Wang (1983) and called there: the mean square slope condition. An illustrative example of the optimal funicular corresponding to a particular transverse load is shown in Fig. 1. The level function of a funicular follows the diagram of the bending moment in a simply supported beam with supports of the same ordinates as points A and B and subject to the same load q(x). Thus, for a given vertical transmissible load of intensity q(x), referred to the unit area of the basis plane, one can construct a family of funiculars and select one of the least volume, see Section 6.1 in Lewiński et al. (2019a). The main conclusion is that in a given plane one can construct frameworks, usually having forms of arches, which are simultaneously tensionfree and do not undergo bending. Admitting tension leads to Michelllike designs of great complexity, see Darwich et al. (2010). A similar theoretical construction in space is much less obvious. A natural counterpart of a planar bendingfree arch is a membrane shell. Of three local equilibrium equations, one is algebraic:
where N_{1}, N_{2} represent the membrane stress resultants corresponding to the principal curvature parameterization, while the principal curvatures are denoted by R_{1}, R_{2}; q represents intensity of normal load, [q] = [N/m^{2}]. Thus, if the normal load is absent and the curvatures are positive, the stress resultants N_{1}, N_{2} must have opposite signs. At least in this case, a membrane shell cannot be subject only to compression and work in a bendingfree state. On the other hand, we remember the solutions of statics of membrane shells of revolution under selfweight: at the pole, both the membrane stress resultants are negative and equal but along the meridian, the diagrams of these stress resultants have different distribution: the meridional force (i.e. its absolute value) increases, while the initially negative circumferential stress resultant increases and, crossing a zero value at a certain circumference, attains positive values along the support, see Flügge (1960). These examples show that, in general, one cannot construct a membrane shell which under a transmissible vertical load would work in a tensionfree state. That is why instead of membrane shells, the tensionfree gridshells are preferred in many applications. The gridshells are tensionfree spatial frameworks formed on a single or on several surfaces, their form being selected such that the bending is vastly eliminated, see Richardson et al. (2013) and Jiang et al. (2018).
Gridshells are usually constructed in two steps. First, the geometry of the deformed configuration is found, e.g. by dynamic relaxation, see Day (1965). Then the optimal position of nodes is found to minimize the total weight and impose additional design conditions, see Richardson et al. (2013). The recent paper by Jiang et al. (2018) discusses three more efficient methods of constructing gridshells: the force density method (FDM) by Linkwitz and Schek (1971) and Schek (1974); the potential energy method (PEM); and the ground structure method (GSM). In the FDM, the values of the member forces are predefined, while the nodal position is determined by the equilibrium conditions in the undeformed state. Thus, the designer has a full control over the values of axial stresses and it is possible to eliminate tension (or compression) a priori. The GSM is also capable of eliminating tension (or compression) but the solution does not constitute a single surface structure; in general, a multistorey structure appears, cf. Lewiński et al. (2019b).
An alternative method of construction of gridshells can be inferred from the theory of archgrids, proposed by Rozvany and Prager (1979) and extended in Rozvany et al. (1980, 1982), Rozvany and Wang (1983), and Wang and Rozvany (1983). This theory has been recently revisited in Lewiński et al. (2019a, Sec. 6.3.6; 2019b), Czubacki and Lewiński (2019), and Czubacki et al. (2018). Archgrids are composed of infinitely thin arches formed in the mutually orthogonal planes: x = const., y = const., orthogonal to the basis Ω corresponding to z = 0. The given load of intensity q is transmissible along the z axis and is decomposed into the loads q_{x}, q_{y} acting on the arches in y = const. and x = const. planes. The shape of the arches is chosen such that the bending is fully eliminated. Consequently, also the transverse shear forces vanish. Thus, the state of stress in the arches is axial: the only stress components are normal to the crosssections. One can design the areas of crosssections of the arches such that the state of stress is fully uniform within the arches and equals − σ_{C}. Of such uniformly stressed structures, one can find the one of the smallest volume. Rozvany and Prager (1979) discovered that this least volume structure is formed on a single surface. More precisely, if z = z_{x}(x,y) represents the surface composed of arches formed in y = const. planes, and if z = z_{y}(x,y) represents the surface composed of arches formed in x = const. planes, then the conditions of minimum volume imply the identity: z_{x}(x,y) = z_{y}(x,y) = z(x,y) there, where both the arches are subject to nonzero loads. On the other hand, arches in one direction can also appear and even some part of the basis Ω can be uncovered. The point loads are also allowed and lead to arches of finite crosssections. Nevertheless, the optimal structure is a grid shell in all cases; the problem of its construction belongs to the class of topology optimization problems.
The numerical solutions shown in the mentioned papers by Rozvany and coauthors have been constructed iteratively, the optimization problem itself had not been explicitly formulated. Czubacki and Lewiński (2019) noted that the archgrid problem reduces to the two mutually dual problems, to some extent similar to those constituting the scalar Beckmann problem, see Santambrogio (2015), forming the Michell’s problem, see Lewiński et al. (2019a) as well as the free material design problems, see Bouchitté and Buttazzo (2001), Czubacki and Lewiński (2015), and Czarnecki and Lewiński (2017). The minimization problem is characterized by the functional of linear growth; its dual version is a maximization problem of the virtual work with locking conditions. Just this decomposition governs the topology optimization features of the setting and of the final form of solutions. Moreover, the numerical and mathematical methods to attack the problem should be appropriately chosen, as for functionals with integrands of linear growth.
The mathematical theory of archgrids evenly supported along the closed contour of the basis domain has been put forward in Czubacki and Lewiński (2019). The present paper is aimed at extensions of this theory towards archgrids unevenly supported. The level of the edge support is determined by the values of the ruled surface given by the equation: z_{0}(x,y) = a_{1}xy + a_{2}x + a_{3}y + a_{4}. In the first step, the surfaces z_{x}(x,y), z_{y}(x,y) are generated by the family of planar arches in the sections y = const., x = const., supported along the contour Γ of the domain Ω on the level z_{0} along Γ. These surfaces determine the grid shells carrying the loads q_{x}(x,y), q_{y}(x,y) such that q(x,y) = q_{x}(x,y) + q_{y}(x,y). The decomposition of the given load q into q_{x},q_{y} is one of the unknowns of the archgrid theory. Upon introducing the areas of crosssections for assuring a uniform state of stress in the whole structure, we state the problem of the volume minimization. This leads to various optimality conditions, one of them being the identity z_{x}(x,y) = z_{y}(x,y) = z(x,y) valid in these subdomains, where two families of arches carry the load. All the optimality criteria will be discussed and the mathematical formulation of the optimum design problem will be stated; it turns out that the level function z(x,y) of the vault is governed by a wellposed problem composing of two auxiliary mutually dual problems. Upon finding the static unknowns, one can compute the optimal areas of the crosssections to assure the uniform distribution of stress. The optimum archgrid is the lightest among all statically admissible and fully stressed.
The least volume bendingfree arch under vertical transmissible load
Construction of the bendingfree arches
The present section shows that for any transverse and transmissible load, one can construct a threehinge arch which is bendingfree.
Consider a threehinge planar arch subject to a vertical load of intensity \(q(\overline {x})\), \(0\leq \overline {x}\leq l\), referred to its projection on the horizontal axis x. The hinges at the support lie on the levels: z = h_{A} and z = h_{B}, see Fig. 2.
Along with the arch problem, an auxiliary beam supported at points lying below the supports A and B of the arch, of ordinates x = 0 and x = l, is considered. The beam is subject to the same load \(q(\overline {x})\), \(0\leq \overline {x}\leq l\).
The vertical reaction at the left support A equals
The bending moment in this beam at the section of ordinate x is expressed by
The condition of the total equilibrium of the arch reads:
where H represents the horizontal reaction at the support A. We consider only such loads q for which H > 0.
Note that the vertical reactions of the arch and the auxiliary beam are linked by
where \(\tan (\alpha ) = (h_{B}h_{A})/l\). The bending moment in the arch at the section of the ordinate x can be expressed by
The shape of the arch can now be adjusted to the given load by requiring that M(x) = 0 for 0 ≤ x ≤ l. This condition determines the shape of the bendingfree arch:
Thus,
where \(\widehat {T}\) represents the transverse shear force in the auxiliary beam of Fig. 2b. This transverse shear force satisfies the differential equilibrium equation
Its variational counterpart has the form
where \(\mathcal {V}_{0}\) represents the set of trial functions vanishing at x = 0 and x = l; the regularity conditions involved in the definition of \(\mathcal {V}_{0}\) will not be discussed in the present paper.
The equalities (9) and (10) imply
called the funicular equation.
One can prove that the arch of the form given by the elevation function (8) is not only bendingfree but is also not subject to transverse shear. It is subject only to compression. The axial force N(x) in this arch is linked with the horizontal reaction H by the formula
where φ(x) is the angle of inclination of the tangent at point (x,z(x)) to the axis x, hence \(\tan (\varphi (x)) = dz/dx,\) cf. Fig. 3. Thus, N(x) < 0 along the whole arch.
The areas A(x) of the crosssections can be chosen such that the normal stresses are made uniform within the whole arch and equal − σ_{C}, σ_{C} being the permissible stress in compression. To this end, we assume N(x)/A(x) = σ_{C}, or
The bendingfree arch of least volume
Among the arches which are bendingfree and uniformly stressed, one can find an arch of the least volume. The volume of the arch of areas of crosssections given by (14) equals
By using the formulas
and considering that the integral of \(\widehat {T}(x)\) along the span of the auxiliary beam vanishes, we obtain
where \(\tilde {l} = (1+\tan ^{2}(\alpha )) l\) and ∥⋅∥ stands for the norm^{Footnote 1}
The optimum design problem (\(\mathcal {P}_{0}\))
Among bendingfree threehinge arches corresponding to a given transmissible load q(x) find the leastvolume arch.
The solution.
All such arches are characterized by the level functions z(x) given by (8), involving the design variable H. By minimizing the volume V over H, we obtain the optimal value of the horizontal reaction \(H = \check {H}\), or
while the volume of the lightest arch equals
Its elevation function \(z=\check {z}\) is given by (8) and (19), or
see Fig. 1. Let us differentiate both sides of (21):
The right hand side has the unit norm; consequently, the norm of the left hand side is equal to 1 and this condition leads to, cf. (1)
This condition means that the optimal arch cannot be too shallow and also cannot be too high. If it were very shallow, the load would cause great values of the thrust H; consequently, the areas of the crosssections would be appropriately big to assure a constant normal stress. On the other hand, if an arch is higher and higher than its volume increases. The condition (23) indicates that the rise of the arch should be perfectly chosen. It is worth noting that this is an integral condition and not a pointwise one.^{Footnote 2} Moreover, note that the bound increases along with the increase of α. Thus, if the difference between the levels of the supports is bigger and bigger, the mean rise increases accordingly. This fact is not obvious and can serve as an hint for a designer.
The arch of minimal volume has the shape given by the elevation function (21), while the areas of its crosssections are given by (14) and (19). The optimal areas of crosssections are expressed as below
Note that \(\check {A}\) cannot degenerate along the arch and assumes finite values.
Remark 1
It is worth noting that the optimal volume \(\check {V}\) (given by (20)) is proportional to the potential energy of the load, namely
The proof runs as follows. Since \(\tilde {z}\) vanishes at the ends of the beam, we can choose \(v=\tilde {z}\) in (11), which implies
Note that
Substitution of (27) into (26) gives
which ends the proof.
Remark 2
Consider the following two auxiliary problems of the variational calculus.
Problem (\(\mathcal {P}_{0}^{\prime }\))
Find the minimizer \(Q = \check {Q}\) of the problem
Problem (\(\mathcal {P}_{0}^{\prime \prime }\))
Find the maximizer \(v = \check {v} \in \mathcal {V}_{0}\) of the problem
If regularity assumptions are appropriately chosen, one can prove that Z_{1} = Z_{2}, cf. Czubacki and Lewiński (2019). Moreover, \(\check {v}\) and \(\check {Q}\) are linked by
and \(\check {Q}\) represents just the transverse shear force appearing in the auxiliary beam of Fig. 2b, or \(\check {Q} =\widehat {T}\). Thus,
where \(\widehat {M}\) is bending moment in the auxiliary beam of Fig. 2b. Moreover, according to (21), the elevation function of the optimal arch can be expressed by
which implies \(\tilde {z}(x) = \sqrt {\tilde {l}} \check {v}(x)\). It is seen that the field \(\check {v}(x)\) attains the bound in the condition involved in (30): \(\ \mathrm {d}\check {v} /\mathrm {d}x\ = 1\).
In the problem considered, the variational problems (\(\mathcal {P}_{0}^{\prime }\)), (\(\mathcal {P}_{0}^{\prime \prime }\)) are not helpful in construction of the solution, since the solution (21) has been directly constructed. Nevertheless, the formulations of problems (\(\mathcal {P}_{0}^{\prime }\)), (\(\mathcal {P}_{0}^{\prime \prime }\)) are inspiration for the construction of the optimal archgrids which are the main subject of the present study.
Archgrids unevenly supported along a closed contour
The aim of this section is to put forward a variationally consistent theory of optimal archgrids.
The stressbased construction of archgrids
The term archgrid is indissolubly bonded with its leastvolume property. It will be shown that just the leastvolume condition leads to shaping a single surface on which the arches (in the planes x = const., y = const.) are formed. The basis of the archgrids to be designed will be a plane domain Ω parameterized by a Cartesian coordinates system (x,y). Assume that the straight lines y = const. cut the contour Γ of the domain Ω at two points: (x_{1}(y),y), (x_{2}(y),y); c ≤ y ≤ d. Similarly, the lines x = const. cut the contour Γ at two points (x,y_{1}(x)), (x,y_{2}(x)); a ≤ x ≤ b, cf. Fig. 4. Assume that the domain Ω is included in a rectangle of dimensions a_{0} and b_{0} and let the point (0,0) be one of its vertices; The sides of the rectangle are given by y = 0, y = b_{0}, x = 0, x = a_{0}, see Fig. 4. Let us form over a given rectangle a ruled surface whose elevation function is expressed by
Note that
The angles of inclination of the lines x = const., y = const. are expressed by:
The vault will be designed over the basis Ω, along the contour Γ. Its elevation along the contour Γ will be given: z = z_{0}(x,y) if (x,y) ∈Γ. The elevations of the vault at the points of the contour Γ: (x_{1}(y),y), (x_{2}(y),y), (x,y_{1}(x)), (x,y_{2}(x)) are denoted by h_{1}(y), h_{2}(y), g_{1}(x), g_{2}(x), respectively. Thus, we have
The following equalities hold
for the subsequent crosssections y = const., x = const., see Fig. 5.
The vault to be designed is assumed as an archgrid composed of two families of arches in the directions y = const., x = const., hence formed in the mutually orthogonal planes. The given load q, referred to a unit area of a basis plane, [q] = [N/m^{2}], being transmissible in the z direction (orthogonal to Ω) is decomposed as follows
The arches in the sections y = const., subjected to the vertical load of intensity q_{x}(x,y), are formed such that they are bendingfree with using the method of Section 2. Consider a simply supported beam of length l_{x}(y). The load q_{x}(x,y) causes the transverse shear force \(\widehat {T}_{x}(x,y)\) and the bending moment \(\widehat {M}_{x}(x,y)\). These fields satisfy the differential equations of equilibrium:
The level function of the arch is assumed such that the arch is bendingfree. Thus, according to (8), the ordinates of the arch formed for y = const. are given by
In the crosssections x = const., the optimal arches, subjected to the vertical load of intensity q_{y}(x,y), are constructed in a similar manner. Consider now simply supported beams of lengths l_{y}(x). The load of intensity q_{y}(x,y) causes the fields \(\widehat {T}_{y}(x,y)\), \(\widehat {M}_{y}(x,y)\) satisfying
The elevation function of the arch is assumed by the rule (8), or
The arch thus constructed is bendingfree. The horizontal reactions H_{x}(y), H_{y}(x) are assumed to be nonnegative.
Note that the fields \(\widehat {T}_{x}\), \(\widehat {T}_{y}\) satisfy the equilibrium equation
At this stage of the study, the designed vault is formed on two surfaces of elevation functions: z_{x}(x,y), z_{y}(x,y). Their values coincide along the contour of Γ.
The values of z_{0} along the contour are determined by the ruled surface (34). In the sequel, our aim is to construct a vault of the minimal volume. It will occur that this condition will result in the identity z_{x} = z_{y} within Ω. Just this identity will be interpreted as a junction of both the families of arches to form a one shelllike structure.
The computation of volumes of the both families of arches requires analysis of their slopes. By differentiating (40) and (42), we get
Let φ_{x}, φ_{y} represent the inclination of the tangents to the arches in the planes y = const., x = const., respectively, see Fig. 6. Thus,
To make the distribution of the axial stress uniform and equal − σ_{C}, we shall choose the densities of areas of crosssections in an appropriate manner. By using the rule (14), one obtains
The elementary volumes of the arches are expressed by
where
Substitution of (47), (49) into (48) gives
The total volume of the arches parallel to the axis x is expressed by
while the total volume of the arches parallel to the axis y equals
Substitution of (46), (45) into (51) results in
which reduces to
where
In a similar manner, we obtain
where
Now we are ready to formulate the main problem of the present paper, namely the problem of constructing the leastvolume archgrid:
The minimizing pairs \((\check {H}_{x}, \check {H}_{y})\), \((\check {Q}_{x}, \check {Q}_{y})\) are stationary points of the Lagrangian \({\mathscr{L}}= {\mathscr{L}}(H_{x},H_{y},\widehat {T}_{x},\widehat {T}_{y}, \lambda )\) given by
where λ(x,y) is a Lagrangian multiplier corresponding to the equilibrium condition (43) and vanishing along the contour Γ. By integrating by parts, we rearrange the Lagrangian to the form
where
Let us note that the condition
is equivalent to (43). Thus, the fields \(\widehat {T}_{x}\), \(\widehat {T}_{y}\) involved in (59) will be treated as independent. Moreover, v = λ/2 plays in (59) the role of a Lagrangian multiplier.
We require \(\delta {\mathscr{L}} = 0 \) with respect to independent variations δv, δH_{x}, δH_{y}, \(\delta \widehat {T}_{x}\), \(\delta \widehat {T}_{y}\). In this way, the stationary conditions of the functional (59) are obtained: the condition (61) and
where the norms ρ_{x}(⋅ ; ⋅, ⋅), ρ_{y}(⋅ ; ⋅, ⋅) are defined by
for a(x) < b(x).
The stationarity conditions (62) and (63) refer to the case of H_{x} > 0, H_{y} > 0.
Substitution of the formulas
into (61) leads to:
Problem (\(\mathcal {P}_{1}\))
Find v vanishing on the contour Γ such that
where the bilinear form at the lefthand side reads
If H_{x}, H_{y} are given and H_{x} > 0, H_{y} > 0, then the solution v of problem (\(\mathcal {P}_{1}\)) is unique. This function v determines uniquely the fields \(\widehat {T}_{x}\), \(\widehat {T}_{y}\) by (65). The formulae (62) may be rewritten as
where now v is the solution to the problem (\(\mathcal {P}_{1}\)).
By combining (68) with (45), one gets
hence,
The functions s_{1}(y), s_{2}(x) can be chosen such that z_{x} = z_{y} in the whole basis domain including the boundary. Indeed by taking
one obtains the identity
cf. (36). Since v vanishes on Γ, we note that z = z_{x} = z_{y} assumes the values z_{0}(x,y) for (x,y) lying on Γ. Thus, if H_{x}(y), H_{y}(x) are given, then the function
determines the shape of the vault. To assure that this shape refers to the structure of the minimal volume, the functions H_{x}(y), H_{y}(x) should be appropriately selected. The relations (63) make it possible to express the total volume V = V_{x} + V_{y} in terms of the fields \(\widehat {T}_{x}\), \(\widehat {T}_{y}\). The direct substitution of (63) into (54), (56) results in
where for w = (u,v) we define
The vector field \(\widehat {\boldsymbol {T}} = (\widehat {T}_{x},\widehat {T}_{y})\) is subject only to the condition (43) (or equivalently: (61)). Now the problem \((\mathcal {P})\) of minimization of the volume of the archgrid can be explicitly formulated as (cf. the theoretical scheme given in Fig. 7):
The problem above is well posed since the set of vector fields Q of given divergence is nonempty. Moreover, let us note that the functional ℘(w) possesses the properties of a norm of a vector function of argument w = (u(x,y),v(x,y)).
Since the functional in (76) is convex and the set of admissible vector fields Q is nonempty, the problem (\(\mathcal {P}^{\prime }\)) is wellposed. The problem of existence depends on regularity of the data.
Let \(\pmb {\check {Q}} = (\check {Q}_{x},\check {Q}_{y})\) be the minimizer of (76). By substituting \(\widehat {\boldsymbol {T}} = \pmb {\check {Q}}\) into (63), one obtains the optimal distributions \(\check {H}_{x}(y)\), \(\check {H}_{y}(x)\) of the boundary forces. Substitution of these results into (67) fixes the coefficients in the bilinear form \(a_{\check {H}}(\cdot , \cdot )\) governing the optimal solution. If \(\check {H}_{x} >0\), \(\check {H}_{y} >0\), then this bilinear form is positive definite and the problem (\(\mathcal {P}_{1}\)) possesses a unique solution \(\check {v}\) which determines the surface of the archgrid by (73).
The displacementbased construction of the archgrid
It turns out that the stressbased problem (\(\mathcal {P}^{\prime }\)) has its displacementbased counterpart and just the latter problem will be, usually, easier to solve. The aim of this section is to put forward this dual formulation.
Let us start by noting that the gradient of \(\check {v}\) can be computed directly from (62) by putting
Let us substitute: \(H_{x} = \check {H}_{x}\), \(H_{y} = \check {H}_{y}\), and
into (63). If \(\check {H}_{x}>0\), \(\check {H}_{y}>0\), then the relations (63) imply
The field \(\check {v}\) determined by the solution of the problem (\(\mathcal {P}_{1}\)) with the optimal values of the coefficients \(H_{x} = \check {H}_{x}\), \(H_{y} = \check {H}_{y}\) of the quadratic form (67) is the maximizer of
Indeed, the problems (\(\mathcal {P}^{\prime }\)), (\(\mathcal {P}^{\prime \prime }\)) are mutually dual; Z_{1} = Z_{2}, or the duality gap vanishes. For the proof, the reader is referred to Czubacki and Lewiński (2019). The proof applies here since the formal structure of problem (\(\mathcal {P}^{\prime \prime }\)) and problem (78) in Czubacki and Lewiński (2019) is the same. It is sufficient to replace l_{x}, l_{y} with \(\tilde {l}_{x}\), \(\tilde {l}_{y}\) and if \(\pmb {\check {Q}} = (\check {Q}_{x},\check {Q}_{y})\) is the minimizer of (76) and if \(\check {v}\) is the maximizer of (80), then these fields are linked by the optimality conditions
Having \(\check {Q}_{x}\), \(\check {Q}_{y}\), one can find \(\nabla \check {v}\) but not vice versa; the relations (81) are noninvertible. Therefore, the construction of the optimal archgrid necessitates solving the two independent problems: (\(\mathcal {P}^{\prime }\)), (\(\mathcal {P}^{\prime \prime }\)), while (81) can play the role of a check. Upon finding \(\check {v}\), one can construct the level function \(\check {z}\) of the optimal archgrid by the rule
where z_{0}(x,y) is given by (34). One can prove that the function \(\check {z}(x,y)\) satisfies the equalities:
These equations are counterparts of the condition (23) known from the theory of planar funiculars. Let us prove the first of the conditions (83). To this end, we insert (cf. (69))
into the first of relations (79). One finds
While computing the l.h.s of (85), one can make use of the formulas:
and thus rearrange (85) to the form of the first of equalities (83).
One can say that the conditions (83) introduce the upper bounds on the slope of the arches in both directions. These bounds are increasing functions of arguments α(y), β(x). The theoretical construction of the optimal archgrids governed by problem (\(\mathcal {P}\)) is reduced to solving problems (\(\mathcal {P^{\prime }}\)), (\(\mathcal {P^{\prime \prime }}\)), their juxtaposition being presented in the diagram in Fig. 7.
Having solved the problems (\(\mathcal {P^{\prime }}\)), (\(\mathcal {P^{\prime \prime }}\)), one can determine the areas of arches’ crosssections, which will be discussed in the sequel.
Determination of the densities of the areas of arches’ crosssections
The solution to the problem (80) determines the shape of the optimal archgrid according to (82) but, in general, does not allow for a direct computation of the densities of areas of crosssections \(\check {A}_{x}\), \(\check {A}_{y}\) in the optimal structure. These fields are determined by the solution \((\check {Q}_{x},\check {Q}_{y})\) to the problem (76). To find the explicit formulas for the densities of the area crosssections \(\check {A}_{x}\), \(\check {A}_{y}\), we shall make use of (47). The forces H_{x}, H_{y} involved in these formulas are found by putting \(\widehat {T}_{x} =\check {Q}_{x}\), \(\widehat {T}_{y} =\check {Q}_{y}\) in (63).
By using (46), (84), and (68), one finds
In this manner, one gets
Now the forces \(\check {H}_{x}\), \(\check {H}_{y}\) can be expressed in terms of \(\check {Q}_{x}\), \(\check {Q}_{y}\) by putting \(\widehat {T}_{x} = \check {Q}_{x}\), \(\widehat {T}_{y} = \check {Q}_{y}\) into (63). The final formulae read
In the above formulas, the arguments in the functions
are omitted. The fields \(\check {A}_{x}\), \(\check {A}_{y}\) are positive everywhere. Having found the minimizer \(\pmb {\check {Q}} = (\check {Q}_{x}, \check {Q}_{y})\) of problem (\(\mathcal {P}^{\prime }\)) and the maximizer \(\check {v}\) of problem (\(\mathcal {P}^{\prime \prime }\)), one can compute the volume of the lightest archgrid either by
or by
which shows that the optimal volume is proportional to the potential energy of the loading.
Optimal archgrids over rectangular domains—the complete numerical approach with using Legendre polynomials
The present section is aimed at numerical constructions of optimal gridworks over rectangular bases, corresponding to selected transmissible loads including the loads concentrated along curves. The construction is based on tackling simultaneously the problems \((\mathcal {P^{\prime }})\) and \((\mathcal {P^{\prime \prime }})\) along with checking the optimality conditions (81). Both the problems \((\mathcal {P^{\prime }})\), \((\mathcal {P^{\prime \prime }})\) are solved numerically by using Legendre polynomials.
Construction of the solution to the problem (\(\mathcal {P^{\prime \prime }}\))
We confine consideration to the case of the bases being rectangular domains [−L_{x},L_{x}] × [−L_{y},L_{y}], hence x_{1}(y) = −L_{x}, x_{2}(y) = L_{x}, y_{1}(x) = −L_{y}, y_{2}(x) = L_{y}.
The n th Legendre polynomial is given by
where \(\left [\cdot \right ]\) is a function which rounds down its argument to the nearest integer. Let us introduce polynomials T_{n}(t), called antiderivatives of P_{n}(t), cf. Abedian and Düster (2019), such that:
and T_{n}(− 1) = T_{n}(1) = 0 for n ≥ 1. The explicit form of T_{n} reads
while T_{0}(t) = t does not satisfy the above boundary conditions.
The test functions in (80) will be expressed by the truncated series involving the T_{n} functions, as below
with \(\tilde {x} = x/L_{x}\), \(\tilde {y} = y/L_{y}\); v_{ij} are unknown coefficients, the numbers m_{x}, m_{y} control accuracy of the method. Due to (94) the partial derivatives of v assume the form:
hence,
or
and the expression \(\rho _{x}\left (\frac {\partial {v}}{\partial {y}} ; y_{1}(x), y_{2}(x) \right )\) can be written in a similar manner. Note that a significant number of components at the r.h.s. of (99) cancel by orthogonality of Legendre polynomials:
where δ_{ij} denotes the Kronecker delta. This simplifies the formulation (80) essentially and paves the way for an efficient programming of the maximization problem, see the r.h.s. of Fig. 8. The boundary condition v = 0 on Γ is satisfied identically. The coefficients v_{ij} are the only unknowns of the problem.
The constraints in (80) are satisfied along the lines: x = x_{i}, \(i=1,\dots ,s_{x}\) and along the lines y = y_{j}, \(j=1,\dots ,s_{y}\), where the numbers s_{x} and s_{y} are fixed and chosen simultaneously with numbers m_{x}, m_{y}. For each x_{i}, \(i=1,\dots ,s_{x}\), the corresponding slopes β(x_{i}) and the auxiliary lengths \(\tilde {l}_{y}(x_{i})\) are computed by (37)_{2} and (57) respectively. For each y_{j}, \(j=1,\dots ,s_{y}\), the corresponding slopes α(y_{j}) and the auxiliary lengths \(\tilde {l}_{x}(y_{j})\) are computed by (37)_{1} and (55) respectively. In the problems considered in the present paper, the choice: m_{x} = m_{y} = m occurred to be efficient.
The numerical approximations to the definite integrals are performed with using the NIntegrate function available in Wolfram Mathematica software.
Construction of the solution of the problem (\(\mathcal {P^{\prime }})\)
The unknowns of the problem (\(\mathcal {P^{\prime }}\)) are expressed by using the Legendre polynomials P_{n} and their antiderivatives T_{n} as follows
and the vertical load is represented by
where
The differential condition nested in (76) makes it possible to eliminate the unknowns \(Q^{y}_{ij}\) by the algebraic equation
thus making the unknowns \(Q^{x}_{ij}\) the only design variables of the problem. The subsequent steps of the method are schematically shown in the flowchart in Fig. 8. Upon finding the minimizer \(Q^{x}_{ij}\), \(Q^{y}_{ij}\), one can either compute the volume by (91) or multiplying by the factor 2 the objective function F, see l.h.s. of Fig. 8. This method delivers directly the distribution of the horizontal forces along the edges, see (63). However the method does not provide a direct formula for the optimal level function \(\check {z}(x,y)\). To find it, one should integrate the equalities (62) and augment them with the boundary conditions. Let us emphasize that the optimal shape is directly delivered by solving problem (\(\mathcal {P^{\prime \prime }}\)).
Archgrids over a rectangular domain for a uniform load
Consider a rectangular domain Ω of side lengths 2L_{x} = 3L and 2L_{y} = 2L. The aim is to construct the optimal archgrids over this domain, of the level function z assuming the values z_{0} given by (34), now denoted by z_{0(λ)}, with (c_{1},c_{2},c_{3},c_{4}) = (λL, 0, λL, 0), λ being a parameter, designed for the transmissible uniformly distributed load q = const.; [q] = [N/m^{2}]. Due to symmetries of the problem v_{kl} = 0 for k, l being even numbers and the integral constraints in (80) are applied only for x_{i} = \((i1)\frac {3L}{2s}\) and y_{i} = \((i1)\frac {L}{s}\) for \(i = 1,\dots ,s\), where s = 2m. Both methods of archgrid constructions are applied according to the flowchart in Fig. 8. The computations based on solving problem (\(\mathcal {P^{\prime \prime }}\)) are performed for m = 15,45,75; the values of the optimal volumes V_{(λ)} (for m = 75) increase along with the increase of λ, cf. Table 1. The increase is nonlinear, cf. the last column of this table. The optimal volume for λ = 0 and L_{x} = L_{y} compares favourably with the result found in Rozvany and Prager (1979), see Fig. 7 therein, and coincides with the result found in Czubacki and Lewiński (2019) by using: (a) the Fourier series method applied to the problem (\(\mathcal {P}^{\prime }\)), (b) the Fourier series method applied to the problem (\(\mathcal {P}^{\prime \prime }\)).
For the same case of λ = 0, the results obtained by the method (\(\mathcal {P^{\prime }}\)), for selected values of n_{x} = n_{y} = n, are listed in Table 2. The values of volumes obtained by the method (\(\mathcal {P^{\prime \prime }}\)) corresponding to the increasing number m_{x} = m_{y} = m form an increasing sequence (a lower converging sequence), see the first row in Table 1, while the sequence of volumes obtained by the method (\(\mathcal {P^{\prime }}\)) for increasing number n_{x} = n_{y} = n is a decreasing sequence (an upper converging sequence), see Table 2. Both sequences tend to the same optimal value, thus confirming the duality gap between problems (\(\mathcal {P^{\prime }}\)) and (\(\mathcal {P^{\prime \prime }}\)) being zero.
Of two methods shown in Fig. 8, the method based on (\(\mathcal {P^{\prime }}\)) produces the sequence of results converging faster. That is why the results in Table 2 could be confined to n = 30, which corresponds to the accuracy of the method (\(\mathcal {P^{\prime \prime }}\)) for m = 75.
The method \((\mathcal {P^{\prime \prime }})\) gives the optimal shapes of optimal archgrids directly. They are not concave, cf. Fig. 9. Let us emphasize that \(\check {z}_{(\lambda \neq 0)} \neq z_{0(\lambda )} + \check {v}_{(\lambda = 0)}\). On the other hand, the method \((\mathcal {P^{\prime }})\) delivers directly the distribution of the horizontal reactions H_{x}, H_{y} and the areas densities A_{x}, A_{y}. The method \((\mathcal {P^{\prime \prime }})\) for m = 75 gives the results of the gradient of v of better accuracy than the accuracy with which the method \((\mathcal {P^{\prime }})\) for n = 30 delivers the quotients Q_{x}/H_{x}, Q_{y}/H_{y}, see Fig. 10 corresponding to the case λ = 0. The same figure shows the precision with which the optimality conditions (81) are fulfilled. To measure the relevant error, the function:
is introduced and the function error_{y}(x) is defined similarly. These errors, corresponding to the sections x = const., y = const., depend on the distance to the edges. The errors are small far from the edges and increase in the boundary zone, see Table 3. This phenomenon appears because in the problem (\(\mathcal {P^{\prime }}\)), the boundary conditions are not explicitly imposed.
The optimal archgrids—other illustrative designs
The optimal archgrids to be constructed in this section will be found by solving the problem (\(\mathcal {P^{\prime \prime }}\)) in order to have directly the optimal level function \(\check {z}(x,y)\). Two kinds of vertical loads will be discussed: the load concentrated along a contour of a circle and the load uniformly distributed within a circle.
The optimal archgrids over a square domain for the load uniformly distributed along a contour of a circle
The optimal archgrids will be constructed over a square domain Ω, of side of length 2L, loaded uniformly along a contour of a circle of origin in the middle of domain and a radius R, see Fig. 11. The load acts along the z axis orthogonal to Ω, the intensity of the load being p = const., [p] = [N/m]. The two cases of the support will be considered:

(i)
z = 0 along Γ

(ii)
z = z_{0} along Γ where z_{0} is defined by (34) with (c_{1},c_{2},c_{3},c_{4}) = (L, 0, L, 0), see Fig. 11.
Due to symmetries of the problem v_{kl} = 0 for k, l being even numbers and v_{kl} = v_{lk} due to symmetry of a square domain. Thus, integral constraints in (80) are applied only for x_{i} = \((i1)\frac {L}{s}\) for \(i = 1,\dots ,s\) and s = 2m.
The accuracy of the result of the optimal volume is controlled by the number m. For increasing m, the volume tends to
cf. Table 4.
The shape of the optimal archgrid supported as in case (i) is surprising, see Fig. 12a. The vault is formed not over the whole basis domain; the four squares of the basis: O_{1}A_{1}B_{1}A_{8}, \(\dots \), O_{4}A_{7}B_{4}A_{6} are not covered, cf. Fig. 12d. The load acting along the circle is transmitted to the sides of the basis domain by straight arches forming four ruled surfaces over the domains: A_{1}A_{2}B_{2}C_{1}B_{1}, \(\dots \), A_{7}A_{8}B_{1}C_{4}B_{4}, cf. Fig. 13. The net of bars over the interior square B_{1}B_{2}B_{3}B_{4} forms a surface of nonzero Gauss curvature, see Fig. 14. This surface is continuously connected with the ruled surfaces mentioned above by the ruled surfaces over the domains: B_{1}C_{1}B_{2}, \(\dots \), B_{4}C_{4}B_{1}. Note that although the archgrid is not concave, all the arches are in compression and the stress distribution is uniform, as assured by appropriate design of areas densities (89). Let us stress that the optimization process has determined the position of the load p over the contour C_{1}C_{2}C_{3}C_{4}. The points \(B_{1},\dots , B_{4}\) are at the highest levels (see Fig. 14), which may suggest that the load is concentrated at these vertices, which is obviously not the case. Let us stress here that the exact analytical solution of this result is still pending.
The optimal archgrid for the case (ii) of the support has two symmetry planes, see Figs. 15 and 16, where the crosssections y = const. are shown. We note that the arches touching the contour of the basis form the ruled surfaces as in case (i). These ruled surfaces and the surface over the interior region B_{1}B_{2}B_{3}B_{4} are connected by other four ruled surfaces, see Fig. 17. The optimal volume equals
and this value is bigger than that of case (i) by 1.26 %. Thus, the construction is heavier, as expected.
The optimal archgrid over a square domain for the load uniformly distributed within a circle
The aim is to construct an optimal archgrid over a square domain of the side of length 2L for the load of intensity q uniformly distributed within a circle of origin at (0,0) and radius R = 0.75L, see Fig. 18a. The level function \(\check {z}\) is constructed by using Legendre polynomials by the method described above.
The optimal shape of the archgrid (Fig. 18b) is similar to that of Fig. 12a, but now the central part is concave and not convex as before. This structure is composed of four ruled surfaces A_{1}A_{2}B_{2}C_{1}B_{1}, \(\dots \), A_{7}A_{8}B_{1}C_{4}B_{4} extended towards the interior of the basis domain by the surfaces B_{1}C_{1}B_{2}, \(\dots \), B_{4}C_{4}B_{1}—which here are not ruled surfaces—and completed by the surface of nonzero Gauss curvature over the B_{1}B_{2}B_{3}B_{4} square (Fig. 19).
For increasing m, the volume tends to
cf. Table 5.
Final remarks
The theory of optimum design of trusses teaches us that the condition of minimum compliance with a prescribed value of the volume leads to designs in which all the bars are stressed up to the same given limit (cf. Lewiński et al. (2019a), Section 2.1). This suggests that the design optimization of archgrids based on the uniform stress condition should lead to the designs of least compliance. Yet this conjecture needs deeper studies due to issues concerning deformation behaviour of archgrids; note that a single funicular is geometrically variable, which makes the kinematic analysis nontrivial.
The identity of both the surfaces z_{x}(x,y) = z_{y}(x,y) of two families of arches shown in Section 3 refers obviously to the undeformed state. Yet the arches of both the families are not connected at nodes on this surface. Consequently, they may form two different surfaces upon deformation. Yet the choice of the areas of the crosssections of both the families of arches assures uniform stress distribution. We conjecture that this property singles out the rational design which assures the required deformation compatibility. In the forthcoming papers, this question will be addressed by a rationalization of archgrids by replacing them by frames of finite number of arches connected by stiff joints. Like Michell structures (see He and Gilbert (2015)), the archgrids may be effectively rationalized, thus showing their lightness, stiffness, and applicability in the civil engineering practice.
The theory put forward concerns the archgrids composed of two families of arches in x = const and y = const planes. Such constructions can be improved, noting that the load applied close to the corners is better transmitted by arches in the planes going in ± 45^{∘} directions, cf. Sec.6.3 in Lewiński et al. (2019a) and examples in Lewiński et al. (2019b) constructed by a discretized version of the archgrid model admitting arches in four directions at the point: 0^{∘}, 90^{∘}, ± 45^{∘}. A continuum model of archgrids composed of arches going in selected directions, or even in all possible planes, is the subject of the current research.
Notes
 1.
The unit of ∥u∥ is different than the unit of u, namely \([\u\] = \sqrt {\mathrm {m}}[u]\), where m is the unit of length.
 2.
References
Abedian A, Düster A (2019) Equivalent Legendre polynomials: numerical integration of discontinuous functions in the finite element methods. Comput Methods Appl Mech Engrg 343(5):690–720
Bendsøe MP (1983) On obtaining a solution to optimization problems for solid, elastic plates by restriction of the design space. J Struct Mech 11:501–521
Bouchitté G, Buttazzo G (2001) Characterization of optimal shapes and masses through MongeKantorovich equation. J Eur Math Soc 3:139–168
Czarnecki S, Lewiński T (2017) On material design by the optimal choice of Young’s modulus distribution. Int J Solids Struct 110–111:315–331
Czubacki R, Lewiński T (2015) Topology optimization of spatial continuum structures made of nonhomogeneous material of cubic symmetry. J Mech Mater Struct 10(4):519–535
Czubacki R, Lewiński T (2019) On optimal archgrids. In: Altenbach H, Chróścielewski J, Eremeyev VA, Wiśniewski K (eds) Recent developments in the theory of shells. Springer Nature, Cham, pp 203–225
Czubacki R, Dzierżanowski G, Lewiński T (2018) On funiculars and archgrids of minimal weight. In: 41st Solid Mechanics Conference book of abstracts (SOLMECH2018); 398–399, solmech2018.ippt.pan.pl/BookOfAbstracts.pdf (access: 31.03.2019)
Darwich W, Gilbert M, Tyas A (2010) Optimum structure to carry a uniform load between pinned supports. Struct Multidiscip Optim 42:33–42
Day A (1965) An introduction to dynamic relaxation. The Engineer 219(5688):218–221
Flügge W (1960) Stresses in shells. Springer, Berlin
Fuchs M B, Moses E (2000) Optimal structural topologies with transmissible loads. Struct Multidiscip Optim 19:263–273
He L, Gilbert M (2015) Rationalization of trusses generated via layout optimization. Struct Multidiscip Optim 52:677–694
Jiang Y, Zegard T, Baker W F, Paulino G H (2018) Form finding of gridshells using the ground structure and potential energy methods: a comparative study and assessment. Struct Multidiscip Optim 57:1187–1211
Lewiński T, Sokół T, Graczykowski C (2019a) Michell structures. Springer, Cham
Lewiński T, Czubacki R, Dzierżanowski G, Sokół T (2019b) Optimal archgrids revisited: variational approach and numerical methods. In: Guo X, Huang H (eds) Advances in structural and multidisciplinary optimization, Proceedings of the 13th World Congress of Structural and Multidisciplinary Optimization (WCSMO 13). ISBN: 9787894372079. Dalian University of Technology Electronic & Audiovisual Press, Dalian, pp 196–201
Linkwitz K, Schek H (1971) Einige Bemerkungen zur Berechnung von vorgespannten Seilnetzkonstruktionen. IngenieurArchiv 40(3):145–158
Niordson F (1983) Optimal design of plates with a constraint on the slope of the thickness function. Int J Solids and Struct 19:141–151
Richardson J N, Adriaenssens A, Coelho R F, Bouillard P h (2013) Coupled formfinding and grid optimization approach for single layer grid shells. Eng Struct 52:230–239
Rozvany G I N, Nakamura H, Kuhnell B T (1980) Optimal archgrids: allowance for selfweight. Comput Meth Appl Mech Eng 24:287–304
Rozvany G I N, Prager W (1979) A new class of structural optimization problems: Optimal archgrids. Comput Meth Appl Mech Eng 19:127–150
Rozvany G I N, Wang C M, Dow M (1982) Pragerstructures: archgrids and cable networks of optimal layout. Comput Meth Appl Mech Eng 31:91–113
Rozvany G I N, Wang C M (1983) On plane Pragerstructures I. Int J Mech Sci 25:519–527
Santambrogio F (2015) Optimal transport for applied mathematicians. Calculus of variations, PDEs, and modeling. Springer Int. Publ., Switzerland
Schek H J (1974) The force density method for form finding and computation of general networks. Comput Meth Appl Mech Eng 3:115–134
Wang C M, Rozvany G I N (1983) On plane Pragerstructures II. Nonparallel external loads and allowances for selfweight. Int J Mech Sci 25:529–541
Funding
The paper was prepared within the Research Grant no. 2019/33/B/ST8/00325 financed by the National Science Centre (Poland), entitled: Merging the optimum design problems of structural topology and of the optimal choice of material characteristics. The theoretical foundations and numerical methods
Author information
Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Replication of results
The maximization problem (\(\mathcal {P}^{\prime \prime }\)) and the minimization problem (\(\mathcal {P}^{\prime }\)) have been solved with using FindMaximum and FindMinimum functions respectively, available in Wolfram Mathematica 11.3.0.0 software, on a computer equipped with the Intel Xeon CPU E52697 v3 @ 2.60GHz (2 processors), 48 GB RAM, 64bit Windows 7. The computations have been partly performed symbolically and partly numerically, see the diagram in Fig. 8. Along with the increase of the order of polynomials represented symbolically, one should increase the WorkingPrecision of the internal computations within the Mathematica software. It has turned out that the choice of WorkingPrecision of 100 digits assures a sufficiently high accuracy of the results. This number of digits was chosen accordingly to the order of polynomials involved in the process of computations. In Sections 5.1 and 5.2, the objective function (80) has been computed with the help of function NIntegrate, while in Section 4.3, the definite integrals have been computed by the function Integrate.
Responsible Editor: Mehmet Polat Saka
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Czubacki, R., Lewiński, T. Optimal archgrids: a variational setting. Struct Multidisc Optim 62, 1371–1393 (2020). https://doi.org/10.1007/s0015802002562y
Received:
Revised:
Accepted:
Published:
Issue Date:
Keywords
 PragerRozvany archgrids
 Michell structures
 Leastvolume design