Simply supported Michell trusses generated by a lateral point load

The paper deals with the optimum design problem posed by George Rozvany: find the lightest fully stressed truss transmitting a given concentrated force to two supports forming a line parallel to the force. One of the supports is a hinge while the second one is a roller. The feasible domain is a square domain over the line linking the supports. The problem thus formulated belongs to the class of the three force problem, till now unsolved in general. In the problem stated here two of the three forces are mutually orthogonal. The family of solutions to this problem is parameterized by the coordinates of the point of the force application, hence is a two-parameter family. This seemingly simple problem generates a vast family of extremely interesting solutions, some of them being known, some being only partly resolved, while others turn out to be surprising and not resolved till now. The present paper delivers exact solutions to the optimum designs corresponding to the force position being a sufficiently big distance to the line linking the supports. The kinematic and static approaches are used, both leading to the same exact results. Other solutions are constructed numerically by the adaptive ground structure method.


Introduction
The present paper refers to the three force problem in the following meaning: find the lightest plane structure of a bounded stress level capable of bearing three selfequilibrated forces acting in the plane at fixed application points. This problem belongs to the class of Michell's (1904) problems. It involves many parameters and its general solution has not been found till now. A selected class of solutions to this problem has been discovered by H.S.Y. Chan (1966), see Rozvany (1998) for a review. The problem of transmitting an arbitrarily directed point load applied between two point supports also belongs to this class; then the points of application of the three forces (two of them being reactions) lie along one line. New and surprisingly complicated highly accurate numerical solutions to this problem have been recently reported by He and Gilbert (2015).
The present paper is aimed at delivering the solutions to a subclass of the three force problem posed to the present authors by George Rozvany: find the Michell truss transmitting a given horizontal force to the two point supports, one of them being a roller, see Fig. 1.
The solution depends on the feasible domain and the position of the point of application of the force, see Fig. 2, where the exemplary family of numerical solutions is displayed, all of them referring to the feasible domain being a square domain over the line linking the supports. This family is driven by the point load, applied at many possible places within the feasible square domain. The subdomain where the point load can be applied will be denoted Fig. 1 Posing George Rozvany's problem by P . Some of these layouts are well known from the literature, some of them have been only announced, reported without a derivation. Yet some of the layouts shown have never been reported.
For the sake of enabling a direct reference to the previous papers, especially to Sokół and Lewiński (2010), the problem will be posed as in Fig. 3; the line of supports is a vertical line, while the point load of magnitude P is applied upwards. Note that the axes x, y of the Cartesian system are different than those used in the paper cited above, which should not lead to misunderstandings. The unknown structures are trusses of finite or infinite number of members. The axial stress σ in the members is assumed to be bounded: |σ | ≤ σ p ; σ p being the stress limit common for tension and compression George Rozvany's problem to be discussed reads: find the least volume truss within a given feasible domain which, under the condition |σ | ≤ σ p at each member, Fig. 3 Scheme of notation of the problem to be solved coordinate transmits the given force P to the hinge support R and to the roller at H, see Fig. 3.
In Rozvany's problem as originally posed the feasible domain is a square: 0 ≤ x ≤ h, 0 ≤ y ≤ h denoted by 1 . In the present paper we shall discuss Rozvany's problem in this domain or we shall extend this domain to make the problem simpler. The solution depends on the coordinates (x N , y N ) of point N of the application of the force, hence it depends on two non-dimensional parameters: Two possible cases: a) ξ 2 + η 2 ≤ 1, b) ξ 2 + η 2 ≥ 1 will be discussed separately. The preliminary numerical results of Rozvany's problem were reported at the CMM2011 conference (cf. Sokółand Rozvany 2011). For the special case of the point load applied at the middle point of the feasible square domain, or for ξ = 1/2, η = 1/2, the  optimal layout as well as the numerical layout have been announced (without derivations) in Fig. 4 in Lewiński and Sokół (2014). This is probably the first Michell truss in which a curvilinear fan appears. The problem of minimizing the weight of the truss under the condition |σ | ≤ σ p reduces to two auxiliary mutually dual problems. The primal problem is formulated as a minimization problem of an integral of a norm of the statically admissible stress within the two-dimensional domains and along the reinforcing bars, see (2.6)-(2.7) in Graczykowski and Lewiński (2006a). This formulation will be discussed in Sec. 6. The dual problem reduces to maximizing the virtual work L(w) over the virtual displacement fields w = (u, v) being continuous and satisfying the locking-like condition pointwise in the feasible domain: (1.2) Here B is the unit ball: the set of second rank tensors ε with bounded eigenvalues: |ε I | ≤ 1, |ε II | ≤ 1. The volume of the optimal truss is expressed by V = 1 σ p max {L(w) |w being kinematically admissible and The equivalence of both the problems mentioned: primal and dual has been proved in Bouchitté et al. (2008) with using subtle tools of the variational calculus. In Rozvany's problem considered the condition of kinematic admissibility of a field w means its continuity in , vanishing at point R of both its components and vanishing of its x-component at point H. In the present paper we shall make use of the construction of this field based on the two pointwise conditions at each point x of the domain m occupied by the optimal structure. The above conditions mean that the pair (ε I , ε II ) lies at the vertex (1, −1) of the locking locus B. The maximization operation in (1.3) imposes the strain field to lie as far as possible from the origin (0, 0) being the middle point of B. In the problems in which shear deformation prevails, and Rozvany's problem is of this type, the conditions (1.4) are satisfied in the whole material domain m identically. The field w should be continuously extended beyond this region such that the condition (1.2) is fulfilled in the whole feasible domain . Note, however, that such an extension is not unique. Such extensions will not be constructed in the present paper, since the question of correctness has been resolved by repeating the computation of the volumes with using the stress-based formula mentioned above. In all cases both the results occur to be the same up to arbitrary accuracy. The duality gap is thus zero, which ends the proof of the correctness of the constructions put forward. The conditions (1.4) are highly restrictive: they determine very specific shapes of trajectories of the virtual principal strains. These trajectories will be parameterized by the curvilinear system forming Hencky nets (α, β), known from the theory of plasticity. In this paper we shall make use of the known parameterizations of these nets. However, the construction of the virtual displacement field w satisfying (1.4) for Rozvany's problem is only partly known. We shall perform one construction of the virtual field w to make the paper complete. To this end we make use of the analytical methods by Carathéodory and Schmidt (1923), Chan (1966, 1967), Hemp (1973 and by the techniques developed in Lewiński (2006a, b, 2007) and Sokół and Lewiński (2010).
The present paper delivers both analytical and numerical solutions to a wide subclass of Rozvany's problem.
We depict the domains of point load application for which the solutions constructed are fully complete, justified by the criterion of the duality gap being zero and confirmed numerically by the newest version of the ground structure method. It turns out that convexity is one of the properties which determines the range of applicability of the analytical layouts discussed.
To make the denotation of references to frequently cited papers as short as possible we shall adopt the following convention: Eq. (n.m) from the paper by Sokół and Lewiński (2010) will be referred to as (SL n.m). Similarly, equation Eq. (n) from the paper by Lewiński et al. (1994aLewiński et al. ( ) (or 1994b will be referred to as Eq. (LZRa n) (or Eq. (LZRb n)). Moreover Eq. (GLa n.m), Eq. (GLb n.m) refer to Eq. (n.m) in Graczykowski and Lewiński (2006a, b).
The Hencky net is formed by curvilinear parametric lines (α, β). The curvilinear coordinates of a specific point are usually denoted by (λ, μ). The virtual displacements along the (α, β) coordinate lines are denoted by (u, v). The displacements along the axes x or y are denoted by w x , w y . Points will be denoted by straight fonts, like A, B, R, N, H, while functions will be denoted by italic fonts, like F , G. In particular, Lamé coefficients are functions and will be denoted by A, B. Frequently used two-argument Chan's functions F n , G n are defined as in Eqs (LZRa 1-2).
2 Parameterization of the feasible domain. Case of the domain P being ξ 2 + η 2 ≤ 1 The structure to be constructed should be supported at the fixed hinge R and should have a roller support at H, with the slide along the vertical direction. The structure should lie to the right of the line RH (the line x = 0) and admit members tangential to this line. Let us note that the point N of the application of the force P lies in the quarter of the circle of origin at R and of radius h: P = {(x N , y N ) | x 2 N + y 2 N ≤ h 2 , x N ≥ 0}. It turns out that a broad class of solutions is encompassed by the layout shown in Fig. 4, composed of: a) the circular fan RAB b) the circular fan NAC c) the domain ABDC generated by two circular arcs d) the domain BDH tangent to RH, called Chan's domain The numerical solutions will show that the solutions to George Rozvany's problem considered assume the form of Fig. 4 if the point N lies sufficiently close to the axis x (y = 0). A detailed description of this condition will be given in Section 7.
In the circular fan domains RAB and NAC the polar coordinate systems are introduced: (r, β) and (r, α) where β ∈ [0, θ 2 ] and α ∈ [0, θ 1 ]. In these domains the virtual displacement field w = (u, v) is constructed: u -along the radii and v -along the circles of r = const.
The parametric lines (α, β) in the domain ABDC are extensions of the radial lines running from within the fans RAB, NAC. The geometry of the parametric lines (α, β) in the domain ABDC has been for the first time mathematically described by Carathéodory and Schmidt (1923). This construction can be found in Hemp (1973), cf. Lewiński et al. (1994a). The parametric lines are defined with respect to the Cartesian system (x 0 , y 0 ) of origin at A, see Fig. 4. These formulae are given by Eqs (SL 2.5-2.7). The Lamé coefficients in the domain ABDC at a point of coordinates (λ, μ) are denoted by A(λ, μ), B(λ, μ) and expressed by Eqs (SL 2.4). The extension of the parameterization (α, β) into the domain BDH has been elaborated by Chan (1967). We repeat here the Eqs (GLa 7.30) expressing the parameterization referred to the Cartesian coordinate system (x 0 , y 0 ) These equations determine Lamé coefficients A(α, β),B(α, β), cf. Eqs (GLa 7.31) Let us note that the Lamé coefficient B vanishes along the line BH, since there β = α + θ 2 and one can easy check that B(α, α + θ 2 ) = 0. Let us add that the points (x 0 , y 0 ) lie along the line RBH if The parameters r 1 , r 2 , θ 1 , θ 2 determine the position of point N by the equations x N = r 1 cos θ 2 + r 2 sin θ 2 y N = −r 1 sin θ 2 + r 2 cos θ 2 (2.5) The equilibrium condition Qh − P x N = 0 yields The coordinates (α, β) of point H are (θ 1 , θ 1 + θ 2 ). The equations (2.1), (2.2) determine the coordinates (x H 0 , y H 0 ) of the point H. One can prove that the substitution λ = θ 1 , μ = θ 1 + θ 2 into these formulae results in the expressions for x H 0 , y H 0 which satisfy (2.4) identically. This proof will be omitted. In particular we find By inverting (2.5) and substitution the results into (2.8) one finds a new form of the equation (2.9)

Construction of the virtual, kinematically admissible displacement field. Case of the domain
The point N of application of the force is assumed to lie within the quarter of the circle of the origin at R and of the radius h or ξ 2 + η 2 ≤ 1. The field w = (u, v) is constructed consecutively: starting from the domains RAN, RAB, NAC, through the domain ABDC and then is found in the domain BDH by subsequent extensions of the solutions through the boundaries AB, AC, BD. The field w should vanish at point R and its horizontal component should vanish at H. To simplify the analysis we shall retain only the first condition and determine w up to an angle of rigid rotation around R. Thus the expression for the virtual work involves now two terms and the final result is obviously independent of the angle of rigid rotation around R. Just the expression (3.1) will be used to compute the optimal volume of the structure by the formula (1.3). The operation max in (1.3) will be performed over the parameters r 1 , r 2 , θ 1 , θ 2 of the design predicted in the form as in Fig. 4. These parameters are linked by the constraints determined by the position of the force and the position of the support at H. To compute the optimal volume one should construct the fields (u, v) in the whole feasible domain and compute w N y , w H x . The fields (u, v) have already been built in the triangle RAN and in the domain RBDCNA, see Sokół and Lewiński (2010), but its extension along BDH has never been constructed. This extension will be performed in the present paper. The displacement fields (u, v) in the domain NAC are given by Eqs (SL 2.55). The parameter ψ involved there represents an angle of rigid rotation of the whole structure around R. The value of ψ does not affect the value of the virtual work, if computed by (3.1). The value of ψ will be chosen to be equal −1, as in the papers by Chan (1966) and Sokół and Lewiński (2010). Under this choice the virtual displacements of point N, referred to the Cartesian system Thus the virtual work (3.1) is expressed by where v H is the displacement of point H measured along β coordinate line. Computation of the virtual work is thus reduced to computing v H . To this end one should construct the virtual field (u, v) first in ABDC and then in BDH. The fields (u, v) in ABDC are given by Eqs (SL 2.70, 2.72); these results correspond to ψ = −1.
To extend these fields towards the domain BDH we make use of the integral representation, see Eqs (LZRa 74-75) involving the angle of infinitesimal in-plane rotation ω The angle of rotation of the whole structure around R is equal ψ = −1, hence ω(0, 0) = −1. We insert this value into the formula for the field ω, see Eq.
(LZRa 73) to find This result will be used in (3.4). The Lamé coefficient B in the domains ABDC and BDH can be expressed as below, cf. Eq. (GLb 6.21) Note that the first of the above formula coincides with the formula for B in the domain ABDC, see Eq. (SL 2.4). Substitution of (3.6)- Note that the expression (3.10) can be written as follows and because B * (λ, β) is equal to B(λ, β) determined within the domain ABDC, hence the function v * (λ, μ) is expressed by Eq. (SL 2.72), although now this expression applies to a different range of the arguments (λ, μ). Thus Now, to find v 1 (λ, μ), we insert (3.6) and (3.8) into (3.11) and perform integration over β. Thus we arrive at the result being a sum of four terms: and the function ω(λ, β) is given by (3.6). The above integrals can be explicitly expressed in terms of functions F n , G n . According to Eqs (LZRa 16) one finds The remaining integrals are expressed by while the functions g 5 , g 6 are defined by Eqs (LZRa A.11, A.12). The above results make it possible to express the function v 1 (λ, μ) in terms of F n , G n functions. By substituting this result into (3.9) and using (3.13) one arrives at the final expression for the field v: Having the fields v and B one can construct the field u by Let us recall the definitions of the fields They satisfy the differential equations In the domain BDH the field v 0 (λ, μ) is given by It is seen that the second of Eqs (3.22) is fulfilled, since the functions F n , G n satisfy the same equation, see Eq. (LZRa 5). Let us recall that the function B(α, β) vanishes along the ( 3.25) which makes the computation of v H easier. We obtain By using (2.7), (2.8) one can rearrange the above formula to the form Note that the expression in the braces is equal to the Lamé field B in the domain ABDC, see Eq. (SL 2.4). Thus the formulae (3.3) and (3.27) determine the virtual work L. It is expressed in terms of r 1 , r 2 , θ 1 , θ 2 . The optimal volume is computed by performing the maximization operation in (1.3) over admissible values of r 1 , r 2 , θ 1 , θ 2 . More details connected with this topic are discussed in Section 7.  The parametric equations of the Hencky net in ABD are given by (2.1) in which while in the domain BDH the Hencky net is given by (2.1) with The Lamé coefficients A(α, β), B(α, β) in the domain ABD are and in the domain BDH: Moreover (cf. (2.7)-(2.8)) h = r 2 (F 0 (θ 1 , θ 1 + θ 2 ) − F 2 (θ 1 , θ 1 + θ 2 )) (4.5) and tan θ 2 = ξ/η. To find the fields u, v in the domains ABD, BDH it is sufficient to put r 1 = 0. In particular we can express the value of v at H The layout in the form as shown in Fig. 5 is determined by three numbers: h, x N , y N ; they determine the values of r 2 , θ 1 , θ 2 . The formula (4.5) is the transcendental equation to find θ 1 . We see that the substitution of r 2 , θ 1 , θ 2 into (4.6) determines the optimal volume by (1.3) uniquely, since maximization runs over one-parameter set.
The analysis here applies also to the peculiar case of θ 1 = 0; then h = r 2 and v H = −(1 + 2θ 2 )r 2 . The optimal layout has the shape of a sector of a circular domain.

Construction of the Hencky nets and virtual displacement fields. Case of the domain
Consider the case of ξ 2 + η 2 ≥ 1 corresponding to the point N lying outside the circle of the origin at R and of radius h. The optimal layout is composed of a circular fan RAB and the adjacent region BAN, see Fig. 6. The geometry of Hencky net and displacement field of the region BAN are derived in Lewiński et al. (1994b, see Fig. 1) and here we recall only the final formulae transformed to the new coordinate system.
The parameterization of the region BAN presented in Fig. 6 leads to the following relation between Cartesian and angular coordinates (cf. Eqs (LZRa 58), ): Note that the angles α and β in BAN are restricted by α ∈ [0, θ 1 ] and β ∈ [0, min(α, θ 2 )], and that parametric lines α are tangent to the line y = h. Moreover, the limiting angles θ 1 and θ 2 can be determined directly from (5.1) for a prescribed position of node N, i.e.
It should be noted here that the solution to the above transcendental system of two unknowns θ 1 and θ 2 is in general not unique, but the uniqueness holds for physically admissible angles θ 1 ∈ [0, π/2] and θ 2 ∈ [0, θ 1 ].
According to the applied parameterization, see Fig. 6, displacements of a point lying in the region BAN, referred to the local curvilinear coordinate system α, β are given by (cf. Eqs (LZRb 203-204)) Consequently, the displacements of the node N referred to the global Cartesian system x, y are expressed by and finally, the volume of the Michell structure of the layout presented in Fig. 6 is determined by the following formula Note that in the present case of ξ 2 + η 2 ≥ 1 the displacement field (5.4) of the region BAN satisfies the boundary condition at point B = H, i.e. w H x = u(α, β)| α=0, β=0 = 0. Thus contrary to the previous sections and especially to (3.1) we do not need to include the work Qw H x . Moreover, all formulae presented in this section cover the limiting case ξ 2 + η 2 = 1 giving θ 2 = 0, which means that the optimal structure reduces to a circular fan RAB.
Remark 1 The layouts of Figs. 4, 5 and 6 do not encompass the whole family of the exact solutions to Rozvany's problem discussed. The numerical experiments reveal that if point N lies close to the line RH, then point N becomes surrounded by the fibers and the layouts become different, more complicated. The subdomain of P of possible places of point N for which the layouts of Figs. 4-6 correspond to the exact solutions will be described in Section 7.

The optimal volume expressed by the stress fields
In Michell structures the mass is either distributed in the plane or along the lines or curves along which the reinforcing bars are placed. The tensor field N = (N ij ) represents the in-plane stress within the domains in which the mass is distributed in the plane, the eigenvalues N I , N II being called principal in-plane stresses. In the reinforcing bars lying along the curve where the mass is concentrated the stress resultant is represented by the axial force F . The limit stress is σ p . The condition of constant distribution of stress in all structural members links the volume of the structure with the stress field. The aim of the present section is just to express the volume of the optimal structure in terms of the stress fields. This formula for the volume will be compared with the kinematic formula (1.3). We shall show that for the optimal structures of the shapes outlined in Figs 4, 5 both the formulae produce the same results up to arbitrarily assumed accuracy. Thus the duality gap is zero, which confirms correctness of both the approaches: kinematic and static ones. Fig. 4 The volume of the optimal structure of the shape of Fig. 4 These fields represent the in-plane forces measured per a unit angle. They satisfy the equations

Optimal structures in the shape of the layouts of
The fibrous stress fields (6.2) and the axial force F in the reinforcing bar are linked along the edge by the known equilibrium equations of plane arches.
Remark 2 As in the theory of homogenization of elastic composites the theory of Michell structures introduces the stress fields at two different observation levels: at the macrolevel (or in the domain ) and at the microlevel (i.e. within the representative volume element, RVE). At the microlevel the stress state is strip-wise uniaxial and the stresses in the strips are equal either to −σ p or to σ p ; the units of this stress state is Pa=N/m 2 . On the other hand we deal with the non-uniformly distributed stress field N at the macrolevel (with units N/m), and, for the sake of convenience, we deal with the fibrous stress fields T 1 , T 2 (or Hemp's forces) of units N, defined by (6.2). Since the reinforcing bars inevitably occur to transmit point loads to point supports, the axial force F in these members appears (of unit N). The cross section of these bars is chosen such that the axial stress is equal to −σ p or to σ p , like in strips at the RVE level. Note that we do not face the problem of stress concentration: the norm of the stress state within the domain at the RVE level and in the reinforcing bars equals σ p . Thus the stress is not only bounded everywhere, but made uniformly distributed, like in the least compliant trusses of minimal weight.
The equilibrium problem of the structure of Fig. 4 is governed by the same equations as the Hemp's fields and the axial force F in the cantilever discussed in Sec. 6 of the paper by Graczykowski and Lewiński (2007). On the basis of the latter results we report below all the formulae for the stress fields in the structure of Fig. 4.

The interior of the fan RAB The radial forces read
The interior of the fan NAC The radial forces read The interior of the domain ABHDC Thus the formulae above are valid for two adjacent regions: ABDC and BDH. The line BD is a line of their continuous extension.
The edge bar HDCN The axial force in this bar is constant and equal to F = −Q.
The edge bar BH The force F in this bar is varying, since each line α approaches the line BH in a tangent manner giving a subsequent contribution to the value of F . The points along BH have the coordinates (λ, λ + θ 2 ), λ ∈ [0, θ 1 ]. The axial force along BH equals where T 1 (α, α + θ 2 ) |BH is given by (6.6) or The edge bar RB The force F in RB is constant and equals The edges RA and NA are not reinforced by bars. The point loads at R and N are transmitted by the axial forces in bars RB and NC as well as by the radial forces in the fans RAB and NAC. Let us write down the equilibrium equations of the nodes N and R. The equilibrium equations of the node N (along the directions x 0 , y 0 ) read where the resultants of the radial forces in the fan NAC are Upon performing integration in (6.11) one reduces (6.10) to the formulae Note that if both equations (6.12) are satisfied, then the condition (2.9) is fulfilled identically.
The equilibrium conditions of the node R, along the directions x and y, are − Q + S x 0 sin θ 2 − S y 0 cos θ 2 = 0 F RB − P + S x 0 cos θ 2 + S y 0 sin θ 2 = 0 (6.13) where are the resultants of the radial forces in the fans along the axes (x 0 , y 0 ). The integrals have been found, but the final formulae are fairly complicated and will not be reported. The above formulae determine the volumes of the subsequent segments of the structure of Fig. 4.

The volume of the interior of the domain ABDC
where T 1 , T 2 , A, B refer to the domain ABDC, see Eqs (6.6) and Eq. (SL 2.4).

The volume of the interior of the domain BDH
where T 1 , T 2 are given by (6.6) while A, B -by (2.3).

The volume of the bar NC
The volume of the bar CH The volume of the bar BH The volume of the bar RB The volume of the whole structure equals The volume V can be expressed in terms of P , r 1 , r 2 , θ 1 , θ 2 , but this formula is long and is not reported. Minimization of V over admissible parameters determines the volume of the optimal structure.

Optimal structures in the shape of the layous of Fig. 5
The optimal structure of Fig. 5 differs from the structure in Fig. 4: the fan NAC degenerates to a point A = N = C, but simultaneously a new stiffener arises: the bar RN (or RA) of finite cross section, transmitting the axial force denoted by F RN . Let us write down the equilibrium equations of the nodes: A = N = C and R. The equations of equilibrium of the forces acting at the node N = A = C, in directions x 0 , y 0 , read where U x 0 , U y 0 are resultants of the force T 2 along the line β in the curvilinear fan around the point A = N = C, see Fig. 5. These forces are equal to T x 0 , T y 0 given by (6.11), where T 2 (α) is given by (6.6) for β = 0. The second condition of (6.24) has the form of the second condition of (6.12) while the first of (6.24) assumes the form Let us write down the equilibrium conditions (along x and y) of the forces acting at the node R where S x 0 , S y 0 are given by (6.14). Thus we have all necessary formulae to determine the stress state within the structure of Fig. 5.
The volume of the bar RN equals The volume of the structure of Fig. 5 equals RN (6.28) and in the appropriate places one should substitute r 1 = 0. The equations (6.23), (6.28) refer to the case of ξ 2 + η 2 ≤ 1.

Discussion of the solutions generated by varying position of the lateral force P
The optimal layouts of Rozvany's problem strongly depend on the position N of the force P . The domain of possible points N has been denoted by P . The main division of this domain is its separation by the circle: (x N ) 2 + (y N ) 2 = h 2 or ξ 2 + η 2 = 1. The case of ext P given by ξ 2 + η 2 ≥ 1, ξ ≤ 1, η ≤ 1 has been fully solved, see Sec. 5. The case of int P defined by ξ 2 + η 2 ≤ 1 has turned out to be much more complicated. The numerical experiments reveal that the layouts of Figs. 4 and 5 describe the exact solutions correctly, if point N lies outside the shaded region * P of int P , see Fig. 7. The region int P \ * P determines the scope of applicability of formulae derived in Secs 2-4 and 6. This region can then be divided into two sub-regions 1 P , 0 P of non-vanishing (r 1 > 0) and vanishing (r 1 = 0) fan NAC. As explained before, the limiting case of ξ 2 + η 2 = 1 refers to the layouts composed of one circular fan RAB.
The line separating regions 0 P and 1 P is more difficult to obtain. Note that in the region 1 P the Eq. (6.12) hold, which means that angles θ 1 and θ 2 depend only on ξ . As before on the border line there holds θ 2 = arctan(ξ/η), thus for a fixed value of ξ ∈ [0, 1] one can calculate θ 1 and η from (6.12) giving the line η = η(ξ ). The angle θ 1 Fig. 8 Selected analytical and numerical solutions for: a ξ = 0.4, η = 0; b ξ = 0.6, η = 0.4; c ξ = 0.8, η = 0.6; d ξ = 0.9, η = 0.9 (here V a represents the value of the volume found analytically, while V n represents its numerical approximation) is redundant here but cannot be eliminated in an algebraic way because (6.12) are transcendental and implicit equations with regard to θ 1 . In other words, the procedure to obtain η = η(ξ ) can be performed only numerically.
For now the present authors cannot precisely determine the line separating regions * P and 1 P so in the Fig. 7 this line is marked 'approximately' by a dashed straight segment which, however, corresponds well with the line separating different layouts obtained numerically.
Let us recall that the data are: P , x N , y N , h, σ p . In case of P = int P \ * P the optimal solution is predicted to be of the form shown in Fig. 4. The geometry of the optimal structure can be found by solving the dual (maximization) problem (1.3) or by solving the primal (minimization) problem (6.1). It turns out that the former problem is easier. Thus the strategy is to solve the dual problem and then insert the optimal parameters into the stress-based expression (6.1), specified in the forms (6.23) or (6.28) to check whether the volumes found by the mentioned formulae coincide. Thus in the first step the optimal volume is computed by performing maximization over w in the expression (1.3) with L given by (3.3) and (3.27). The field w is parameterized by the geometric quantities r 1 , r 2 , θ 1 , θ 2 linked by the three equations: (2.5) and (2.9). Thus the virtual work L becomes a function (r 1 , r 2 , θ 1 , θ 2 ) and the problem (1.3) can be reduced to maximize 1 σ p (θ 1 , θ 2 , r 1 , r 2 ) subject to Eqs (2.5) and Eq (2.9) and r 1 , r 2 , θ 1 , θ 2 ≥ 0 (7.3) and then solved numerically with using for example Find-Maximum procedure available in Mathematica. It is worth noting that the solution to (7.3) is based only on the kinematically admissible displacement field. If we take into account the equilibrium conditions of the node N (or R) we can get an additional equation which eliminates the necessity of performing the maximization procedure. In the problem investigated in the present paper the design variables r 1 , r 2 , θ 1 , θ 2 can uniquely be determined from the system of the following equations: (2.5) 1,2 , (2.9) and (6.12) 1 .
The exact solutions described by the formulae derived in Secs 2-6 have been verified numerically using the adaptive ground structure method put forward in Gilbert and Tyas (2003), then gradually developed and improved (see Sokół 2011, 2014, and Sokół and Rozvany GIN 2015. The comparison of exact optimal volumes and layouts with the corresponding numerical results clearly indicates the correctness of the formulae derived. The series of optimal solutions parameterized by (ξ i , η j ), with ξ i = i/10, η j = j/10; i = 1, 2, . . . , 10, j = 0, 1, . . . , 10 were performed using the ground structures of 200 × 200 cells with almost 500 million potential non-overlapping bars. Exemplary solutions are shown in Fig. 8 and Table 1. It has been mentioned earlier that in the region 1 P the angles θ 1 and θ 2 depend only on ξ and this interesting property can easily be observed in the first three rows of Table 1, where the values of θ 1 and θ 2 are constant. For fixed ξ and increasing η the radius r 2 increases while r 1 decreases to zero, where the domain 0 P begins. For more detailed comparison of numerical and analytical results for various ξ and η the Here in the third column corresponds to the circle line separating the domains 0 P and ext P (see Fig. 7). Note also that for ext P the values of r 1 , r 2 and θ 1 , θ 2 are interchanged with respect to the notation in Fig. 6 to assure smooth changes of design parameters along the line Reader is referred to the Electronic Supplementary Material attached to this paper.
The numerical solutions corresponding to the shaded domain * P for selected ξ and η are presented in Fig. 9. They considerably differ from the layouts of Figs 4 and 5, since the point N is surrounded by a fibrous material domain.
Finally, the volumes of the optimal structures are checked by the stress-based formulae (6.23), (6.28),with the design parameters being maximizers of (7.3). The results thus computed compare favorably with those found by (7.3), which proves correctness of all resultsreported.
The virtual work equals L = 1.7021587177P h, which gives the optimal volume V = L/σ p according to the kinematic method. The optimal volume can alternatively be computed by using the formula (6.28) based on the static analysis. By using the formulae of Sec. 6.1, putting there r 1 = 0 and by (6.27) one computes all the components in the formula (6.28) contributing to the total volume of the structure. The volumes of the subsequent parts of the structure are: Summing up these values one finds V = 1.7021587177P h/σ p , exactly the same volume as that found by the kinematic formula.
Similar computations can be shown for all other optimal layouts reported in the present paper. This check is of vital importance taking into account that the formula (1.3) necessitates the construction of the virtual field w in each case in the whole feasible domain, while our constructions of the virtual field have been restricted to the domain, where the mass is present. Such extensions can be performed, but as they are non-unique they have only theoretical value.
Remark 3 Basing on the partial volumes presented above one can note that V CH = V RAB . This interesting property is not accidental and results from the relationship between integrals (6.15), (6.16) and (6.20). It can be proved that for a more general case of r 1 0 (i.e. for both domains Fig. 9 Selected numerical solutions beyond the applicability range of the theoretical layouts predicted in the present paper 1 P and 0 P ), the equality V CH = V RAB + V NAC holds. It means that the total volume of the infinitesimally thin fibrous bars in the fans RAB and NAC is equal to the volume of the rod CH of a finite cross section. The proof is omitted here.

Final remarks
The present paper shows that George Rozvany's problem can be viewed as a problem of finding the optimum structure to turn an applied load through 90 degrees onto the roller support. The paper shows that the layouts predicted in Figs 4, 5 are correct if the point load lies in appropriate distance to the line linking the supports. For closer placement of the point load, in the domain * P , the optimal layouts extend those predicted; the fibers around the load application point appear to strengthen the structure. The optimal layouts become non-convex.
The solutions concerning the domain * P are complicated due to the condition of the feasible domain lying to the right of the line of supports. Yet if one removes this condition and allows for placing the material left to this line, then immediately the optimal solutions become clear, easier for analytical description, see Fig. 10. The analytical formulae of Secs 2-4 suffice to construct the exact solutions; they will be the subject of a separate paper.
Topology optimization methods, and Michell's theory in particular, are helpful in designing the connection sections in thin-walled shells subjected to tangent point loads of large magnitudes, see the recent paper by Zhang et al. (2015). The aim is to transmit a given concentrated force towards the main part of the shell to make the interaction forces as uniform as possible. Let us note that the solution to the problem of Fig. 6 gives a hint of how to solve this optimal diffusion problem, since it delivers the construction of the analytical solution to the problem of equilibrating a given vertical force of magnitude P 0 by the system of n vertical forces of the vertical total resultant equal P 0 , acting at points along a horizontal line in a distance from the point where the force P 0 is applied, and forming a system of forces symmetric with respect to the line along which the force P 0 is acting. The number n can be arbitrarily large. The Michell truss being solution to this problem is constructed by superimposing the optimal solutions corresponding to the triplets of forces (P i , −2P i , P i ) acting at points N , R, N, where N is the image of point N by the reflection across the y 0 axis. In the problem considered, see Fig. 6, the points N and N lie on the x 0 axis. Consequently, the optimal structure carrying this i-th triplet of forces is constructed by the reflection of the Michell structure of Fig. 6 across the y 0 axis. Superposition of all the Michell trusses corresponding to the triplets of vertical forces (P i , −2P i , P i ) results in the optimal structure equilibrating the system of applied vertical forces of magnitude P i by one vertical force of magnitude P 0 . This construction comprises the case of the force systems in which one of the forces is applied at the middle point of the interval N N (or at B, cf Fig. 6); the only change in the layout is that one more straight bar appears linking the points B and R. The Michell solution described above has been confirmed numerically by the ground structure method, see Fig. 11 corresponding to 11 vertical forces of magnitudes: P /2, P, . . . , P, P/2 equilibrated by one vertical force of magnitude P 0 = 10P . To model a uniformly distributed load the magnitudes of the forces lying at extreme points are two times smaller than magnitudes of the other 9 forces. The optimal layout includes the Hencky nets strengthened by internal ribs linking point R with points of applications of the applied forces. Moreover, a horizontal non-prismatic bar appears to assure the equilibrium. Similar solutions for arbitrary number of vertical forces can be constructed thus The Michell-like solution of Fig. 11 neglects the interaction of the bars with the plate or the shell in which this structure is embedded. This interaction can be modelled by the methods used by Zhang et al. (2015) or those proposed by Zegard and Paulino (2013).
Let us note that the solution of Fig. 11 has much in common with the designs proposed in Zhang et al. (2015), assuring a uniform diffusion of concentrated forces applied along the generating lines of a thin cylindrical shell. The present authors do cherish the hope that the exact solutions to George Rozvany's problem will be also helpful for solving the optimal force diffusion problem in the case of the load acting parallel to the edges.