New analytical benchmarks for topology optimization and their implications. Part I: bi-symmetric trusses with two point loads between supports

Extending some results by Sokół and Lewiński (Struct Multidisc Optim 42:835–853, 2010), the optimal topology of bi-symmetric trusses with two symmetrically placed pointloads is determined analytically, and verified by highly accurate numerical calculations. It is rather remarkable that Michell’s best-known classical solution had to wait over hundred years for its generalization from one to two point loads. Some implications of these solutions, including properties of so-called O-regions, are also discussed.

connecting two supports of given location (Fig. 1a). The feasible domains will be chosen as rectangles of various heights h; the limiting case of h being infinite. As can be seen from Fig. 1b, h 1 denotes the vertical distance of the highest point of the optimal truss from the horizontal axis of symmetry, if the structural domain is the full plane (or if the optimal height of the truss is smaller than the height of the feasible domain (h 1 < h in Fig. 1b)). Only one quarter of the structural domain, together with symmetry conditions, is shown in Fig. 1b. By recently introduced symmetry principles (Rozvany 2011), the optimal topologies to be discussed are equally valid for two pin supports, or for a pin and a roller support.
The present paper differs from those by Lewiński (2010, 2011a), because in the latter the feasible domain was a half-plane (above or below the line of supports) and there was only a vertical axis of symmetry. In the present paper (i) the feasible domain is a rectangle or the full plane, and we have two axes of symmetry, (ii) the adjoint strain field is also given for the memberless central region of some of the optimal layouts.
Various cases of truss topology problems with two symmetrical loads between two supports (pin or roller) and for different domains (full or half-plane) were recently reviewed briefly by Sokół and Lewiński (2011b).
The line notation in this text is as follows. Thick and thin continuous lines denote concentrated and 'distributed' truss members, broken lines domain boundaries, dotted lines region boundaries and dash-dot lines axes of symmetry or skew-symmetry.
For trusses with a size-and sign-independent stress constraint and a single load condition, the specific cost function is where A is the member cross-sectional area, F is the member force, k is a constant, and σ p is the constant permissible stress. The known optimality conditions (Michell 1904;Prager and Rozvany 1977) for the above problem arē ε = ksgnF (for F = 0) , |ε| ≤ k (for F = 0) , (2) whereε provides the 'adjoint strain field', which must be kinematically admissible (satisfying kinematic support and continuity conditions). The optimal structural volume can also be calculated from the 'dual formula': where the vector P denotes the external forces and u the adjoint displacements at these forces. The optimal adjoint strain fields for plane Michell trusses may consist of the following types of regions (e.g. Prager and Rozvany 1977; similar regions were defined and used for grillages, see (Rozvany 1972)): T-region with a tensile and a compression member at right angles:ε 1 = −ε 2 = k, S-region with members having forces of the same sign in any direction:ε 1 =ε 2 , |ε i | = k (i = 1, 2), R-regions with only one member at any point: |ε 1 | = k, |ε 2 | ≤ k, or O-region with no members: |ε 1 | ≤ k, |ε 2 | ≤ k, where subscripts 1 and 2 indicate principal strains. The concept of O-regions was introduced already in the 1970s (e.g. Rozvany and Hill 1976). A special case of O-regions, in which both principal strains are zero (rigid region withε 1 = ε 2 = 0), will be termed a Z-region. Another special case of O-regions, so-called V-regions (with |ε 1 | ≤ k,ε 2 = 0 ) is introduced in this paper (see Section 2.1).
Referring to Fig. 1, the following types of optimal topologies exist: Exact analytical solutions with numerical confirmation will be presented for Topologies 1 and 2, but only high resolution numerical solutions for Topology 3, for which the analytical solution is not yet available.

Details of optimal topologies
2.1 Topology 2/3 (limiting case between topologies 2 and 3) This topology with h = d/ √ 2 is discussed first, because for this case we have a complete adjoint strain fields for all analytical solutions, some of which are shown in Fig. 2. The optimal solution for Topology 2/3 depends on the relative magnitude of (L − d) and h.
If (L − d) = 0, then we have the solution in Fig. 2a, which is the classical solution by Michell (1904). The strain fields for the quarter domain consist of two T-regions: one has constant principal directions, the other one consists of a circular fan. This solution is in fact valid even when h > h 1 , because the T-regions shown in Fig. 2a can be extended even to the entire plane, satisfying the optimality condition of Michell in (2).
The layout in Fig. 2a will become part of all other topologies in Fig. 2b-d (shown partially on the left hand side of these diagrams). The right hand side of these layouts consists of concentrated horizontal bars along the domain boundaries. It is important to note that these optimal topologies are only valid if the domain has a limited height, i.e. h = d/ √ 2. This is because in Fig. 2b, for example, kinematic continuity would not be satisfied if we moved upwards the present domain boundary (the strains in the two T-regions would cause an overlapping of the horizontal displacements). However, since the adjoint strain field in O-regions may be non-unique (see Section 4.2), there should be some other adjoint displacement field for the present problem of In Fig. 2b and d the central part of the adjoint strain field consists of T-regions and Z-regions, see the definitions above. In Fig. 2c, however, there is also an R-region, with a horizontal principal strain ofε 1 = −k. The value ofε 2 can be calculated from the condition that the strains must be zero in the direction of the region boundary with the Z-region. Then elementary calculations giveε 2 = k tan 2 α. This impliesε 2 = 0 andε 2 = k, respectively, for Fig. 2b (with α = 0) and Fig. 2d (with α = π/4).
The solution in Fig. 2b was actually presented by the second author previously (Rozvany 2011).
Extending the region pattern in Fig. 2b to d, we can also construct optimal strain fields for any value of L − d > 2h.  Fig. 3 Optimal adjoint strain fields for Topology 2/3 with 0 < L−d < h For 0 < L − d < h the optimal topology is discussed in greater detail subsequently. For this case the optimal adjoint strain field is shown in Fig. 3a. Each quarter of the central, memberless O-region consists of a T-region and a V-region (see Section 1). The state of adjoint strains in these two regions, respectively, is represented by the Mohr-circles in Fig. 3b and c. It can be seen from Fig. 3b that the strain along the boundary of the T-and V-regions is Moreover, one can infer from Fig. 3c that for the V-region, we have It can be seen that for α = 0 and α = 45 • (5) gives the correctε 1 values for the limiting cases in Fig. 2a and b. The above adjoint strain field satisfies all optimality criteria of Michell (1904) (see also (2) in this paper). With this addition, we have complete analytical solution for all aspect ratios of Type 2/3 topologies, and also for the O-region of all Type 3 topologies.

Topology 1
The topology of the optimal truss for h = ∞ (or h ≥ h 1 , see Fig. 1b) can be constructed as follows. We can infer from numerical results that the optimal layout is that of a 'Michell cantilever' 1 and a horizontal bar. The geometry of this structure is fully determined by three angles: θ 1 , θ 2 and γ 2 , shown in Fig. 4. To make the paper self-consistent, it is useful to outline the so-called Lommel functions U n (·, ·) that play the fundamental role in deriving coordinates and the displacement field in the Hencky net of mutually orthogonal members in region ABDC (see Fig. 4). Following the notation by Lewiński et al. (1994) we will use in the present paper two functions: G n (·, ·) and F n (·, ·), defined as and where I n (·) denotes the modified Bessel function of the first kind. Many important properties of functions G n and F n can be found in the paper mentioned before (Lewinski et al. 1994). Now, the coordinates of the point D, connecting the horizontal bar with the Michell continuum, can be written as: Equations (8) are obtained by rearranging (2.9-11) of the paper by Sokół and Lewiński (2010). In the latter paper the displacement field of the domain RBDCNR (divided by appropriate sub-regions) was derived in detail. In particular, the adjoint displacements of points N and D (normalized with k = 1) are defined by: and where The (8)-(11) are general and valid for any angles: θ 1 , θ 2 and γ 2 . They can be simplified considerably for a specific problem. For the problem under consideration, one can show that the angles θ 1 , θ 2 and γ 2 must satisfy γ 2 = π/4 and θ 1 = θ 2 + γ 2 .
The first relation under (12) comes directly from the zero value of the horizontal displacement at point N: w N x = 0 and (9) 1 . The second constraint can be derived from the fact that the upper chord BD should connect smoothly to horizontal bar DS at the point D, see (5.1) by Sokół and Lewiński (2010) and the detailed explanation given there. Now, by (12) the coordinates of point D can be expressed as while the displacements of the point D can be written as where the auxiliary function is introduced to shorten the notation. Note that both the horizontal and vertical displacements of the point N defined in (9) are equal to zero for γ 2 = π/4. The reason for this strange result comes from the fact that the displacements (9) and (10) are derived by means of a rigid rotation around point R. The magnitude of this rotation is not known in advance but can be determined from the boundary conditions. This procedure was described in detail by Sokół and Lewiński (2011a, see (23) and (24) there). The displacements of points N and D adjusted in this way (without normalization) can be written as and where ψ denotes the angle of the rigid rotation. It can be determined from the zero value of the horizontal displacement of point S, lying at the vertical axis of symmetry (see Fig. 4). We know that for fully stressed truss the virtual elongation of the bar DS is equal to its length (L − x D ), thus From (18) and (17) 1 one can easily deduce that Now, the adjusted (properly rotated) displacement w N y can be written as where x D , y D and w D x are given in (13) and (14). The forces for the layout in Fig. 4 are given by: Q y = P/2, F x = Q y d/y D and F y = 0. The horizontal reaction Q x at point N can be obtained from Chan's formula (see (2.47) in the paper by Sokół and Lewiński (2010)), while the reactions at point R can easily be obtained from the equilibrium equations. These reaction forces, however, do not generate any virtual work in (3), because the corresponding virtual displacements are zero (points R and N in Fig. 4).
The volume of one quarter of the full structure can be calculated by adding the volume of the Michell's continuum and the volume of the horizontal bar DS (see Sokół and Lewiński 2010). However, it is more elegant and convenient to use the dual method (see (3)) The volume in (21) is a function of one variable V = V (θ 2 ) and its minimum can be determined from the necessary condition V (θ 2 ) = 0. This leads to the transcendental equation where ξ = d/L and function q(θ 2 ) is defined by Equation (22) uniquely defines the optimal angle θ 2 because the function q(θ 2 ) is monotonic (it is a decreasing function for θ 2 ≥ 0 which starts from q (0) = √ 2 and then asymptotically approaches 0 for θ 2 → ∞). The lower limit of ξ for which the solution of (22) is valid for the rectangular domain shown in Fig. 1a is equal to ξ = 0.182027. This corresponds to θ 2 = π/4 which means that the upper chord RBDS starts vertically from the support. For lower values of ξ the external fans extend beyond vertical lines drawn above the supports and the solution is formally infeasible. However, if we allow the horizontal expansion of our rectangle domain outside the supports, we can obtain the feasible solution for ξ < 0.182027. In this case the lowest limit of ξ for which the solution of (22) makes physical sense is equal to ξ = 0.0477491. This corresponds to θ 2 = π/2 and θ 1 = 3π/4. For lower values of ξ the internal circular fan goes outside the symmetry line connecting two supports, and that is obviously infeasible. Thus we can conclude that for ξ < 0.0477491 the optimal solution is not known.
Examples of optimal layouts of Topology type 1 are shown in Fig. 5. For increasing ξ the angle θ 2 decreases and at the same time the height and the volume of the structure increase. Obviously, the final solution for ξ = 1 is identical with the well known Michell's solution (1904) (with θ 2 = 0, θ 1 = π/4, and 4V = P L/σ p (2 + π)). All solutions have been confirmed by numerical calculations. The exact layouts are given on the left side and their numerical equivalents on the right.
The exact and numerically calculated volumes are compared in Section 3.

Topology 2
In this case y D = h and then by (13) 2 we have Obviously this case is simpler than the previous one (compare (22), (23) and (24)). For the limiting case h = d/ √ 2 one can easily obtain θ 2 = 0 and then the total volume of the optimal truss becomes The above formula gives the correct volume for special cases. For ξ = 1 we obtain 4V = P L/σ p (2 + π), which is the solution by Michell (1904). For ξ = √ 2/ 1 + √ 2 we obtain 4V = P L/σ p (4 + π) 2 − √ 2 , which can also be readily derived from the second author's solution (Rozvany 2011, here Fig. 2b).

Topology 3
The exact solutions for Topology 3 are not known at present, but the optimal layouts can be predicted on the basis of numerical solutions (e.g. those in Fig. 7). Here the value ξ = d/L = 0.5 was assigned to the length ratio. The layout in Fig. 7a was calculated for 353553L. It is the limiting case between topology types 2 and 3, hence this layout is similar to the one presented in Fig. 6e and f (see also Fig. 2 in Section 2.1). By progressively decreasing the height h of the structural domain for the numerical solutions, we find that additional T-regions develop in the optimal topology, similarly to a long cantilever (Lewiński et al. 1994), but following a different geometry (see Fig. 8).
In Fig. 8 the upper horizontal bar is connected to the chord of a circular fan having an angular width of 45 • and a curved three-sided region, which is similar to that derived by Chan (for details see the paper by Lewiński et al. (1994)). The upper border of this region can be a straight line if the height is sufficiently small (Fig. 8b) or it starts with a straight section and smoothly passes into a curved section (Fig. 8a). It is to be noted that contrary to Topologies type 1 and 2 the external (upper) chord has not a constant cross section in the straight segment. The three-sided domain above the circular fan is connected with a region with straight members in one direction. This is the reason why the rest of the regions with Hencky nets are different from those derived for the long cantilever problem by Lewiński et al. (1994). Nevertheless, there are also some similarities in forming subsequent new regions if the permissible height is decreased. On the basis of many additional numerical tests we have established that switching from the first sub-type (Fig. 8a) to the second one (Fig. 8b) occurs at a height value in the range of h ∈ (0.24d, 0.25d).
The adjoint strain field in empty regions of Topology 3 can be filled the same way as described in Section 2.1 (Fig. 2).

Confirmation of the analytical results by numerical calculations
All exact solutions presented in this paper have been confirmed by numerical calculations. They were performed using an improved version of the program developed by Sokół (2011a). This newer version of the program makes use of the 'adaptive ground structure' approach, similar to that proposed by Gilbert and Tyas (2003), however, it applies a different strategy of adding-removing of the active bars from a huge number of potential members of the full ground structure. In the present version the program is capable of solving large-scale problems with the number of potential members exceeding one billion (10 9 ). This program was announced in 2011 (see Sokół 2011b) and will be presented in detail in a separate paper in a near future. It should be noted that this program uses a linear programming formulation, and therefore the solution obtained guaranteed to be the absolute minimum volume for the given discretization of nodes. As was shown in Figs. 5 and 6, the numerical layouts are very close to the analytical ones. A comparison of the volumes obtained by numerical and analytical calculations is given in Table 1. Even in the worst cases ( Fig. 5a and b) the numerically optimized volume differs from the exact one by no more than 0.1%.

R, S, O and Z regions in optimal Michell topologies
Various types of optimal regions were defined-in the context of grillages-already in the early 1970s (e.g. Rozvany 1972), and the present notation (including O-regions) was used a few years later (e.g. Rozvany and Hill 1976, p. 208).
In a highly creative recent research paper, Melchers (2005) used the notation N and B for our O-and Z-regions. He pointed out that for Michell trusses relatively few known solutions contain regions other than T-type. This is possibly so, because much of the theory of T-regions is based on the mathematically similar theory of slip lines in plasticity, which was developed earlier. However, for Michell structures with line supports, there are in general R-regions in the solution, and T-regions are rather an exception (see Rozvany et al. 1997, e.g. Figs. 5 and 6). In the same paper Z-regions (which are special cases of O-regions) were used in several solutions.
Since, in the definition of O-regions (see above), we have non-strict inequalities, and the essential point is that we have no members in these regions, the central memberless parts of the topologies in Fig. 2b to d can be regarded as O-regions, although they contain T-, Z-and R-regions.  The optimal O-regions containing T-, Z-and R-regions are not necessarily unique for a given optimal Michell truss topology, as was demonstrated earlier (Rozvany et al. 1997).
For the simple problem in Fig. 9, the domain has a restricted height of a and two hinge supports, and is loaded by a horizontal couple. The fairly obvious optimal topology consists of two parallel bars. A simple optimal adjoint strain field is shown in Fig. 9, but an alternative strain field consists of cycloids which was given by Hemp (1973, p. 95). Hemp's corresponding displacement field is which gives u = v = 0 at the two hinges, and the displacements u = ±k L at the forces, which implies the correct structural volume by (3).
The above example demonstrates that the adjoint strain field may be non-unique in a memberless part (i.e. O-region) of an optimal Michell topology.
An even much simpler example of non-uniqueness of the adjoint strain field is shown in Fig. 10. In this problem we have two vertical line supports at a distance of 3L from each other, and a horizontal point load at a distance L from the right hand support. The optimal layout obviously consists of a single horizontal bar between the load and the nearer support. The principal adjoint strain in the vertical direction is everywhere zero:ε y = 0. In Fig. 10a we have an R-region on the right hand side and an O-region on the left hand side, where the inequality in (2) admits a horizontal strain of −k/2. An alternative, discontinuous adjoint strain field is given in Fig. 10b, in which we have a Z-region and an R-region on the left, both admissible by (2).
In Figs. 2, 9 and 10 the optimal truss layout is unique, only the adjoint strain field is non-unique. The above nonuniqueness is, therefore, not to be confused with the finding (Rozvany 2011) that a Michell truss problem may have either one or an infinite number of optimal solutions (of the same weight).
In a recent paper (Rozvany and Sokół 2012, Fig. 6) another example of an O-region is given, which covers a half plane and has principal adjoint strains that have a constant absolute value which is greater than zero and smaller than k. Both Fig. 10a, and this latter example disprove the notion that in O-regions the adjoint principal strains are either zero or k.

Relaxation of optimality criteria
The authors believe that, if there exists at least one feasible solution for a Michell problem, then there exists also an optimal solution satisfying Michell's optimality criteria in (2) exactly. Nevertheless temporary relaxation of optimality criteria (as suggested by Melchers 2005), followed by an optimization within the set of solutions obtained by the relaxed conditions, can be very useful in finding optimal topologies.
An important example for this was Morley's (1966) contention that no optimal solution satisfying known optimality criteria exists for a square clamped domain in flexure (either reinforced plate or grillage). It was the greatest breakthrough in the optimal grillage theory when Melchers (Lowe and Melchers 1972) found the exact solution of the above problem (by temporarily relaxing and then enforcing the relevant optimality criteria). As indicated above, we do have the optimal adjoint strain fields for Type 2/3 and Type 3 topologies for the memberless central parts of these problems, but it will be a highly challenging task to find these for Type 1 and Type 2 topologies, because of their partially curved boundaries.

Concluding remarks
In this paper we have presented exact analytical solutions for a new class of Michell truss problems: two symmetrically placed point loads in between supports, the feasible domain being symmetric to the line of supports. The exact solutions have been confirmed by very close numerical results. A rather remarkable historical aspect of this development is the fact that it has taken over hundred years to f ind the exact (analytical) extension of Michell's bestknown classical (1904) solution to two point loads. Some important general implications of these solutions have also been pointed out.