Design of optimum grillages using layout optimization

Grillages are often used to form bridge decks and other constructions. However, following a period of intensive research activity in the 1970s, comparatively little attention has been paid to optimizing the layout of grillages in recent years. In this contribution a new numerical procedure is proposed which takes advantage of the adaptive solution scheme previously developed for truss layout optimization problems, enabling very large scale problems to be solved. A key benefit of the proposed numerical procedure is that it is completely general, and can therefore be applied to problems with arbitrary loading and boundary conditions. Also, unlike some previously proposed procedures, the sizes of individual beams can readily be discerned. To demonstrate its efficacy the numerical procedure is applied to a range of grillage layout design problems, including load dependent problems which could not be solved using traditional methods. It is shown that important phenomena such as “beam-weaves” can be faithfully captured and new high-precision numerical benchmark solutions are provided.

viewed the plastic design of grillages in a continuous setting, considering a notional slab comprising an infinite number of fibre-like beams. An optimum fibrous slab can be considered as analogue to an in-plane Michell structure, which is well known in the structural optimization research community; further development of the theory of Michell structures (Michell 1904) has been described by workers such as Chan (1967), Hemp (1973) and Lewinski et al. (1994a, b). Also, a numerical means of identifying Michelltype structures using the "ground structure" approach was proposed by Dorn et al. (1964) and further developed by workers such as Gilbert and Tyas (2003), Sokol (2014), and Zegard and Paulino (2014).
Any structural design optimization problem can be posed in either equilibrium (primal) or kinematic (dual) form, where for a grillage the problem variables are usually moments and rotations in the equilibrium and kinematic forms respectively. However, in the paper by Rozvany (1972a) neither the equlibrium nor kinematic problem formulations are solved directly; instead a displacementbased, fully analytical method of finding the solution of the kinematic problem is proposed which essentially stems from the stress-strain optimality relation linking the solutions of the equilibrium and kinematic forms. The associated optimization problem is also confined to being applied to fully clamped slabs subject to an arbitrary, though always downward, loading. Significantly, for this class of problem the optimum layout is load-independent. This remarkable feature, combined with Rozvany's kinematic method, provided a means of obtaining universal exact optimum grillage layouts for problems involving downward loads for both single and multiple load cases. It should however be noted that this does not furnish the optimal distribution of beam widths. For this one obviously needs to know the magnitudes of the particular loads involved, and to use the governing equilibrium equation to determine the corresponding optimal bending moment field and thus the beam widths. However, the papers by Rozvany (1972a) and by Lowe and Melchers (1972) do not describe systematic means of recovering the beam width distribution.
In subsequent decades Rozvany's analytical kinematic method was applied to slabs with a range of other boundary conditions, including simply supported edges and combinations of free and simply supported edges, and also simply supported and clamped edges (Rozvany et al. 1973;Rozvany and Hill 1976;Prager and Rozvany 1977;Rozvany and Liebermann 1994). For each of the aforementioned cases the method proved capable only of solving problems involving exclusively downward loading, as it explicitly relies on the fact that the optimum layout is load-independent. The method was then implemented by Hill and Rozvany (1985) in a computer program which allowed automatic generation of analytical optimum layouts for arbitrary polygonal slabs with partially clamped and simply supported boundaries. The authors presented exact optimum layouts for an impressive range of complex polygonal domain shapes.
It should be noted that although Rozvany's method is capable of treating interior clamped supports, it cannot account for interior simple supports. This is because uplift may occur if such supports are present, which in turn means that there is no longer a universal kinematic solution, common for all types of downward load. This is also the case if mixed downward / upward loadings are present, or if point moment loadings are present. Less trivially, this also applies to slabs with partially clamped and free edges.
The load-sensitivity of many real-world problems encouraged researchers to seek general numerical methods. In the paper by Sigmund et al. (1993) the ground structure method was, apparently for the first time, used to obtain solutions to the grillage compliance minimization problem. (A "ground structure" comprises a network of structural members interconnecting nodes laid out on a grid from which the subset of members defining the optimum structure is sought, after Dorn et al. 1964.) By using the DCOC method in combination with linearly tapering beam finite elements the authors found a number of new grillage layouts, showing solutions for problems involving clamped and free edges; later the method was applied to problems involving mixed downward and upward loadings (Rozvany 1997). Low resolution ground structures were used, in part because of the available computing capabilities of the time. However because the method does not take advantage of modern adaptive solution schemes (e.g. Gilbert and Tyas 2003), the scale of problems that can be tackled even now appears to be limited. More recently Zhou (2009) proposed a method which involved recovering principal moment trajectories, but this is likely to be rather cumbersome in practice.
In summary, although the analytical approach initiated by Rozvany and co-workers had reasonably broad applicability, and allowed new insights to be drawn, it left a wide range of grillage optimization problems unsolvable, due to their inherent load sensitivity. Moreover, even for grillage problems which could be solved, the optimum beam width distribution was not identified. In the present paper the authors propose that a ground structure approach is adopted, and that techniques now well established in the field of truss layout optimization (e.g. see Gilbert and Tyas 2003;Sokol 2014) or limit analysis via discontinuity layout optimization (see Smith and Gilbert 2007;Gilbert et al. 2014) are applied. Here a plastic design formulation is used by posing two mutually dual linear programming forms: equilibrium and kinematic. The goal is to minimize the volume of material for specified applied loading. This leads to a simple linear formulation which can be used in conjunction with an adaptive solution scheme to solve problems involving ground structures consisting of many million beams, thus generating optimum layouts closely approximating the analytical solutions found by Rozvany et al., though with a far greater range of applicability.

Equilibrium formulation
Consider a ground structure (Dorn et al. 1964) consisting of a design domain discretized using n nodes and b beams, as shown on Fig. 1. For beam i, assume its cross sectional area varies linearly from a i1 to a i2 . The total volume V of the structure can be written as: [a 11 , a 12 , a 21 , a 22 , ..., a b1 , a b2 ] T are, respectively, vectors of beam lengths and areas. For each node, moment equilibrium needs to be enforced in the x and y directions and force equilibrium in the z direction. Denoting m i1 and m i2 as the moments at the Ground structure for a design domain, in this case a simple square domain discretized using n = 9 nodes and b = 36 interconnecting beams (including overlaps) two ends of beam i, the local equilibrium matrix can be expressed as: where θ i is the angle of beam i to the x + axis, and l i is its length. Also, assuming that the beam cross-sections are of uniform depth, let m + p and m − p denote the limiting moments per unit area. The yield condition of beam i can thus be written as: The grillage layout optimization problem can therefore be written as:

Adaptive solution scheme
When a fully-connected ground structure is used the number of beams b grows rapidly with the number of nodes n, limiting the size of problem that can be solved (since in this case b = n(n − 1)/2). This issue was addressed for truss layout optimization problems by Gilbert and Tyas (2003) who proposed an adaptive solution scheme, later further developed by Sokol (2014). This scheme employs an initial sparsely connected ground structure and uses constraints from the dual problem to check whether the solution could potentially be improved by adding additional members, as part of an iterative process. As problem formulation (4) is very similar to the truss formulation used by e.g. Gilbert and Tyas (2003), the same basic technique can be applied in this case, where the dual formulation of (4) involves maximizing virtual work: where, W is the virtual work and u collects the virtual rotations in the x and y directions and out-of-plane displacement in the z direction. Also the following constraint must be satisfied: which imposes limits on the maximum and minimum virtual rotation that can occur in each beam. Note that u is obtained automatically after solving (4), and (6) is only guaranteed to be satisfied in beams that are present in (4). This means that potential beams, not currently represented in the problem, may violate (6); in this case the beams most in violation should be added to problem (4) to prevent this violation in the next iteration. The process repeats until no violation is found in (6); for further details of the algorithm readers are referred to Gilbert and Tyas (2003) and Sokol (2014).

Commentary
The grillages considered herein are assumed to be rigidjointed, but with torsional resistance neglected. This assumption, also made by Sigmund et al. (1993), is justified by the fact that for most cross-sections used in practice the torsional resistance is low compared with the bending resistance. This is particularly true for open cross-sections, such as I-beams. In the latter case by varying only the flange width the linearity of the formulation is preserved. Now consider beam i such that m i1 · m i2 < 0, i.e. the bending moment changes sign across the length of the beam. Notice that restricting the cross sectional area function to vary linearly from a i1 = |m i1 | /m p to a i2 = |m i2 | /m p will not lead to an optimal beam being generated in this case, since each intermediate point along the length of the beam is overdesigned (e.g. consider the intermediate point where the bending moment vanishes). This can potentially be addressed in two ways: (i) via use of a nonlinear relation for the grillage volume (1) (see Bolbotowski 2018); (ii) ensuring that a large number of nodes and interconnecting beams are employed in the problem, such that any inaccuracy is small. A drawback with (i) is that it requires the use of computationally expensive non-linear optimizers and thus (ii) is adopted here. However, note that the single load case plastic design problem considered here is only equivalent to the corresponding compliance minimization problem when (i) is adopted; the same holds for the grillage-like continuum addressed by Rozvany et al.
The above argument suggests that, when m + p = m − p = m p , the volume V in (4a) approximates to the scaled (by m p ) integral of the absolute value of the bending moment diagram taken over the grillage. Thus when a high resolution ground structure is adopted the equilibrium form (4) can be viewed as a discrete version of the continuous problem addressed in Rozvany (1972a), Save and Prager (1985) and others.
In the field of truss optimization it is well-established that when a single load case is involved there must always exist a statically determinate optimum truss solution; see for example Achtziger (1997). The simplest proof of this statement revolves around existence of a basic solution to the underlying LP problem. As this also applies to the grillage optimization problem (4), it follows that there must always exist a statically determinate optimum grillage.
Finally, although thus far attention has focussed on single load case problems, it is well known that the plastic truss layout optimization can be extended to treat multiple load case problems (e.g. see Hemp 1973), and, though beyond the scope of the present contribution, it is worth pointing out that the grillage layout optimization formulation described herein can be similarly extended.

Numerical examples
The proposed numerical method was programmed independently in both Mathematica 11.1.0.0 and Matlab 2015a, respectively using the default Mathematica solver and Mosek 7 to obtain solutions to the LP problems involved. In all cases tried the results obtained from the two programs were identical for all quoted significant figures, though for the larger problems the Matlab / Mosek combination was favoured due to lower associated run times. All quoted CPU times are single core values obtained using a workstation equipped with Intel Xeon E5-2680v2 processors running 64-bit CENTOS Linux.
The efficacy of the method is demonstrated through application to a range of numerical example problems. For sake of simplicity beam moment capacities were in all cases taken to be equal for sagging and hogging, i.e. m + p = m − p = m p , and nodes were evenly distributed over each problem domain, with pressure loads approximated using point loads applied at these nodes. In this case the magnitude of each point load was calculated by taking into account the area of the surrounding domain, taking the load applied at an intermediate node along an edge to be half that applied at an interior node, and the load applied at a corner node as one quarter. For example, considering the domain shown in Fig. 1, and assuming a uniform pressure load of total magnitude pL 2 is applied, the loads applied to nodes A, B and I would be pL 2 /16, pL 2 /8 and pL 2 /4 respectively. Note that because of the presence of supports the grillage designed would in this case only need to carry a load of pL 2 /4, leading to an underestimate in the volume of material required. To address this load discretization error, and also the nodal discretization error that limits the range of layouts that can be identified, and hence tends to overestimate the required volume of material, most problems described were solved using a sequence of increasing nodal divisions, enabling approximations of the exact values to be obtained via extrapolation (see Appendix A for details); these latter approximations are quoted in the main text, whilst tabulated results are presented in Tables 1-3 of Appendix A. However, in the interests of visual clarity the graphical results presented correspond to problems with a moderate number of nodal divisions. Beams are drawn in blue and red to indicate sagging and hogging respectively in the graphical solutions, with line widths proportional to beam cross-sectional areas. In the interests of visual clarity, beams with very small cross-sectional areas have been filtered out. Symbols used in the paper are illustrated in Fig. 2. The symbols used by e.g. Rozvany (1972a) for region type are used herein to describe the analytical optimum layouts; see for example Fig. 3a. Specifically, a design domain can be divided into regions where each region is labelled to indicate the optimum directions of beams of possibly non-zero cross section, whether sagging ("+" symbol) or hogging ("−" symbol) moments are involved. The circle symbols denote so called "indeterminate regions", where the optimum beam direction is arbitrary.
Due to the presence of indeterminate regions it was found that the numerical layouts obtained often became complex in form. This is because the interior point method used to solve the underlying LP problem (4) will normally identify a solution that combines all possible designs in  Fig. 4b. To address this, the length vector l can be modified by adding a constant joint cost / length (j = ±10 −6 unit length), i.e.l i = l i + j . Joint costs were first introduced by Parkes (1978) as a simple means of rationalizing optimum trusses. Here a very small joint cost is used to ensure the numerical layout is pushed towards a basic LP solution to increase visual clarity. Furthermore, numerical tests showed that the clearest visual results could be obtained when j is taken as a small positive value for hogging beams and a small negative value for sagging beams. Notwithstanding this, all optimum volumes presented herein were computed without employing a joint cost.

Benchmark examples
A range of example problems are presented, starting with problems for which closed-form analytical solutions are available. It should however be noted that although countless analytical optimum layouts were presented in the papers by Rozvany et al., optimum volume values were rarely quoted. This is due to the fact that loads were generally not specified, since the layouts derived had universal applicability for arbitrary (downward) load. However, since the optimum displacement field can be recovered from an analytical layout, the optimum volume V can be computed from the virtual work done by given load W , although the process can be laborious and hence analytical volumes will only be provided for selected example problems.

Square domain with simple supports
The first example considered herein involves a square design domain with simple supports, as shown on Fig.  3a. This problem is one of the oldest and simplest to derive analytically, e.g. see Morley (1966). The solution is optimum for arbitrary downward load; one square R ++ region is present along with four triangular R +− regions. For comparison an optimum layout for a uniform pressure load was generated via the new layout optimization method, see Fig. 3b. Since the R ++ region is indeterminate the numerical solution presented is in fact one of an infinite number of possibilities, where here the pressure load is transferred to two beams of significantly larger cross section. The four triangular regions appear to be R + type rather than R +− as proposed analytically since only sagging beams are present, with the orthogonal hogging beams vanishing. This apparent discrepancy, along with other subtle issues associated with numerical layout optimization of grillages, will be considered in the next section.

Square domain with clamped supports at corners
The second example involves a square design domain with clamped supports at the corners, as shown in Fig. 4. This serves to illustrate the effect of the joint cost used to rationalize the solution. The analytical layout shown in Fig.  4a is proposed based on the approach described by Rozvany (1972b) for arbitrary downward load. Aside from four R +−

Fig. 3
Square domain with simple supports: a optimum layout derived analytically by e.g. Morley (1966) Square domain with clamped supports at corners: a analytical optimum layout according to Rozvany (1972b); new result obtained by numerical layout optimization for uniform pressure load; b without joint cost regions the optimum grillage comprises five indeterminate regions, four R −− regions and a single R ++ region. For comparison an optimum layout for a uniform pressure load was generated by layout optimization, initially without using a joint cost, as shown on Fig. 4b. The numerical representation of the R +− regions coincides perfectly with the analytical design, whereas the indeterminate regions involve numerous overlapping beams in different orientations, thus rendering the numerical solution of little practical value. However, by re-running the problem with joint costs the solution is greatly simplified, as shown on Fig. 4c. The R ++ region is now transformed to a regular grid and the R −− regions to cantilever fans radiating out from each of the four point supports; similar fans will be observed in the vicinities of clamped point supports (or concave corners of supported boundaries) in subsequent examples. However, it is evident that the introduction of a joint cost has appeared to transform the R +− regions into R + regions, as occurred in the previous example. This is because the optimum beam width distribution is not necessarily unique for a given applied load. Thus in the example shown in Fig.  4 the particular representation of the indeterminate square R ++ region influences whether or not hogging beams are present in the adjoining R +− region. However, the lack of a unique beam width distribution can also be demonstrated in simpler problems, without R ++ or R −− regions. For example, consider the case of two opposing cantilever beams of equal length subjected to a shared load at their tips; also Fig. 12 serves as a further example. Although a statically determinate optimum grillage is guaranteed to exist, the structure shown in Fig. 4b is clearly not statically determinate, due to the non-uniqueness of the solution. The use of a joint cost does not necessarily fully remedy this, e.g. see Fig. 4c.

Square domain with external clamped and interior point supports
The next example involves a square design domain with clamped external supports and four interior clamped point supports. In the paper by Rozvany et al. (1973) the problem was deemed load-independent and consequently an analytical layout was given for all downward loads; see Fig. 5a where sagging or hogging indeterminate regions are depicted in solid ink. Figure 5b shows the new numerical solution obtained, assuming a uniform pressure load is  Table 1 of Appendix A, with the number of adaptive member adding iterations required to obtain a solution for a given nodal discretization shown together with associated CPU times.

Square domain with four column supports
The next example involves a square design domain with free external boundaries and four supporting columns in the interior, represented by square clamped supports. Assuming arbitrary downward load, Fig. 6a shows the optimum layout derived analytically by Rozvany (1972a) for this particular problem. An analytical solution u of the kinematic form can be uniquely derived based solely on the layout from Fig. 6a which is independent of the load, provided this is always downwards. For example, if a uniform pressure load of magnitude p is applied then, based on duality arguments, the exact volume V exact of the optimum grillage can be computed from the virtual work done by the load as follows: where denotes the design domain. Both the function u and the integral are computed in Appendix B.
The numerical solution for this case is shown in Fig. 6b; the optimum volume V num = 13.33pb 4 /m p is derived from the values tabulated in Table 2 of Appendix A. The close correlation between the numerical and analytical solutions is clearly evident, both in terms of computed volume and grillage layout. Note that for this example only two adaptive member adding iterations are required to obtain a solution (see Table 2 in Appendix A) because the initial ground structure already contains most critical members; this also leads to relatively low associated CPU times.

Domains with free and simply supported edges
Problems involving domains with both free and simply supported edges were investigated by Rozvany and Liebermann (1994). The associated optimization problems were challenging mathematically, with the formulas describing the optimum directions of the beams given in implicit integral form which had to be solved numerically.
Here a right-angled isosceles triangle domain with a simply supported base edge and a simple point support in the right-angle corner is considered, as shown in Fig. 7. The optimum grillage layout for this problem found by Rozvany and Liebermann (1994) for arbitrary downward loading is given in Fig. 7a. Note that the layout is not trivial as the beams do not radiate from the simple point support. A numerical solution was obtained for the uniform downward pressure load case and is presented in Fig. 7b. The close resemblance between the analytical and numerical solutions is clear.

Square domain problem
The next example involves a square domain comprising two clamped and two free edges, as shown in Fig. 8a. This problem was previously considered by Rozvany (1972b), though an exact analytical solution was not derived (even then the class of optimum grillage problems for clamped / free boundary conditions was recognized as being difficult). The same problem was revisited by Sigmund et al. (1993), who presented numerical solutions obtained using a ground structure-based approach combined with FEM. Despite the insights these generated a general analytical method for grillages with clamped and free edges has still not been (a) (b) Fig. 7 Triangular domain with free and simply supported edges: a optimum layout derived analytically by Rozvany and Liebermann (1994); b new result obtained by numerical layout optimization for a uniform pressure load, V num = 0.09505ph 4 /m p found, highlighting a clear gap in the grillage optimization theory developed by Rozvany et al. In Fig. 8 a range of problems are solved, for cases involving point and pressure loads.

Domain with hole problem
The next example involves a domain with a hole and free and clamped edges, as shown in Fig. 9a. Solutions were obtained for three different loading scenarios, involving either point or pressure loads. The optimum layouts for the problems involving point loads, presented in Fig. 9b,c, are perhaps of particular interest since they clearly indicate how the load finds its way through an optimum grillage around the hole back to the supports.

Uplift effect
One of limitations of the computer software tool produced by Hill and Rozvany (1985) was that it could not model internal simple supports because of the potential for uplift, rendering the optimum grillage layout dependent on the position of the load(s) involved. An example is shown in Fig. 10. In this case the internal support divides the design domain into two parts: one subjected to pressure load p 1 , and the other p 2 . When p 1 and p 2 are both applied, the optimum grillage layout is shown in Fig. 10b. If only one of them is applied, different results are obtained, as shown in Fig. 10c and d, indicating the load dependant nature of the problem. Uplift effects are present in the problems shown in Fig. 10b and d, where in (b) the load p 1 effectively cancels out some of the bending effects caused by p 2 , leading to a lower volume than in (d).

Partially downward and partially upward load
Applying mixed downward and upward loads yields yet another class of problem for which load-independent optimal layouts cannot be found. Analytical results for a modest range of such problems were published in a short paper by Rozvany (1997); however a general method of treating such problems analytically has not yet been developed. The example presented in Fig. 11 gives a good insight into the nature of such problems; here two point loads are to be transferred to four simple point supports. Figure 11a shows the solution for all-downward load; in this case two separate simply supported beams of total volume V dd = 0.25P L 2 /m p prove to be optimal. However the optimum layout shown in Fig. 11b for the case when one of the loads is upward clearly involves interaction between the two forces, thus considerably reducing the optimum volume, to V du = 0.1875P L 2 /m p .
(e) (f) Fig. 8 Square domain with two free and two clamped edges: a problem definition (domain has dimensions L×L, with point loads applied on the diagonal); b P 1 = P , P 2 = 0, V 1 = 0.1408P L 2 /m p (which can be shown to coincide with the exact solution); c P 1 = 0, P 2 = P , V 2 = 0.4093P L 2 /m p ; d P 1 = P , P 2 = P , V 1+2 = 0.5315P L 2 /m p ; e uniform pressure load p = P /L 2 , V = 0.07067P L 2 /m p ; f "beam-weave" phenomenon

Point moment load
Formulation (4) permits point moments to be applied directly, thus yielding another class of load-dependent problem. An illustrative example involving a rectangular domain with a point moment load remote from a support is shown in Fig. 12; domains of constant width and varying height are considered. The key observation is that, despite filling the entire height of the domain with an optimum layout, the optimum volume remains constant, at V = M L/m p . This indicates the indeterminacy of the layout in each case (since e.g. design (c) is also a viable solution to problems (a) and (b)).

Non-optimal design of beams with end moments of different signs
As mentioned in section 2.3 the solutions obtained via the proposed numerical method will overestimate the true solution in cases where the bending moment function changes in sign along the length of one or more beams. However, none of the optimum layouts presented in section 3 contained any beams where this was the case. Numerical experiments involving other problems showed that when such beams were present, use of a higher nodal refinement remedied this. This is to be expected, since two shorter beams can always be chosen to meet at the point of contraflexure in a long beam, at least approximately.

Load dependent layouts in grillage optimization
The class of grillage optimization problems solved fully and analytically by Rozvany et al. share an essential property: independence of the layout from load, providing the latter is always applied in a downward direction. This means that there is a displacement vector u that solves dual problem (5) for every downward load, or alternatively, that there exists a displacement vector u that maximizes the out-of-plane displacement of every point (node) simultaneously. In contrast the optimal layouts for the problems considered in Sections 3.2 to 3.5 are load dependent, and there are currently no analytical methods that can be applied to such problems. In this context it is worth revisiting the triangular domain problem initially investigated in Section 3.1.5. According to Rozvany and Liebermann (1994) the analytical layout shown in Fig. 7a should be universal for all downward loads. However, this can be checked by using the numerical method developed herein to explore a range of different loading scenarios. Thus consider for example the case of  Fig. 13a. Comparing Fig. 13a with Fig. 7a it is evident that the beam directions differ, suggesting that this problem is not load independent after all. To verify this finding the theory of grillage-like slabs can be invoked; see e.g. Rozvany (1972a). With the given point load P duality theorems can be used to show that the analytically derived layout of Fig. 7a yields a lower bound volume V ≈ 0.18P h 2 /m p . Conversely the numerical solution of Fig. 13a is associated with a oneline bending moment field M that furnishes an upper bound volume V = 0.25P h 2 /m p . In order to prove that the exact volume V exact = V , and the associated exact moment field M exact = M, it is sufficient to guess a displacement function u such that curvature constraints are met and the optimality relation between M and u holds, i.e.
-the principal curvatures κ I , κ II produced by u satisfy the point-wise inequalities: −1/m p ≤ κ I , κ II ≤ +1/m p and -the left free edge is one of the principal trajectories of curvature field κ and the principal curvature κ I is equal to +1/m p along this edge respectively. Naturally the function u must also satisfy the support conditions. It can easily be verified that the function in question can be given by an extremely simple closedform expression (where x, y are Cartesian coordinates, as indicated in Fig. 13a): which implies M exact = M, u exact = u and V exact = V = 0.25P h 2 /m p , and further that the numerical solution given in Fig. 13a is in fact the exact solution for the grillage-like slab problem with a single point load. The u field resembles a slab being twisted around the y axis, as shown in Fig.  13b. In fact, the same field u is also found when a point moment is applied at the point support; Fig. 13c shows the corresponding layout. (It now becomes clear that a function of the form of (8) also furnishes a solution to the problem described in Section 3.5.) It is also of interest to now consider the case of two point loads applied symmetrically midway along each the free edges; this yields a volume V = 0.354P h 2 /m p and the layout shown in Fig. 13d. Here the optimum layout appears to be inscribed within the analytical layout proposed by Rozvany and Liebermann (1994); see Fig. 7a. Note that the optimum volume is considerably smaller than double the volume of the one-beam solution. However, if the magnitudes of the applied loads are changed, a non-symmetrical numerical layout is obtained; see Fig. 13e. Here the orientation of the sagging beams forming the fans noticeably diverge from the analytical layout shown on Fig. 7a.
These numerical experiments, together with the analytically proposed function u ultimately show the loaddependence of the triangular domain problem initially investigated in Section 3.1.5. The load-dependence appears to be due to the same uplift effect that occurs in the problems considered in Section 3.3, where in that case the axis of uplift was an internal line of simple support. In the triangular domain problem every straight line passing through the point support and the interior of the domain is a potential uplift axis. From this argument it can be concluded that the presence of a simple point support, either placed on the boundary or in the interior of the design domain, is likely to lead to load-dependence in the grillage optimization problem. Note that the triangular domain problem was the only example given in Rozvany and Liebermann (1994) that considered a simple point support; the authors' focus was originally a class of problems with free and simply supported edges only, so all other optimum layouts derived therein can be assumed to be truly load-independent and hence correct.
Finally, suppose that the triangular domain of this problem is transformed into a trapezium to allow the point support to be replaced with a very short simply supported edge. The solution when a single point load applied midway along the free edge is shown in Fig. 13f. It is evident that the new optimum layout now appears to be in agreement with the analytical layout derived by Rozvany and Liebermann (1994). This suggests that the anomaly in this case stemmed from an assumption that an infinitely short line of simple support could be taken to be equivalent to a point support. However, the former prevents rotation about the y axis, and allows a reaction moment about the same axis to be generated. This appears to be crucial in order for the optimum grillage to comprise beams which coincide with the analytical solution shown in Fig. 7a.

Beam-weave phenomenon and including torsion
In the solutions shown in e.g. Fig. 8e,f a thin region of orthogonally intersecting sagging and hogging beams occurs along each free edge. This phenomenon has previously been identified in the optimum grillage layouts found analytically for problems involving mixed free and simply supported edges (Rozvany and Liebermann 1994); the result from Fig. 7 with R +− -type free edges serves as an example. This "beam-weave", as it was called therein, is particularly difficult to approximate using the ground structure approach since, theoretically, it is supposed to be infinitely thin. Similarly, a beam-weave turns out to be an optimum means of transferring load along free edges; e.g., see Fig. 9. Here the role of the beam-weave becomes more apparent; essentially it attempts to mimic a single beam capable of transferring torsion. Consequently, one can observe that limiting the height of the domain in the problem shown in Fig. 12 would provide an infinitely thin, beamweave-like design which is essentially equivalent to a single member in pure torsion.
The above suggests that it may be worthwhile to include torsion in the problem formulation after all, since the beamweave regions degrade the quality of the numerical layouts. However, the underlying problem formulation then becomes nonlinear, and means of obtaining a suitable linearized approximation of the problem will be the subject of future research.

Conclusions
A new numerical layout optimization method capable of identifying the minimum volume and associated optimal layout of a grillage has been proposed. Beam members which are tapered along their lengths between nodes have been employed to maintain the linear character of the problem. This means that highly efficient linear programming algorithms can be used to obtain solutions, with the adaptive "member adding" technique previously applied to truss layout problems enabling solution of large-scale problems, containing large numbers of nodes and interconnecting members. A key feature of the new method is its generality; it can be applied to problems involving arbitrary domain geometries and loading and support configurations, and can faithfully capture important phenomena such as "beam-weaves", which provide resistance to torsion when individual beams have negligible torsional resistance.
When applied to problems for which exact analytical solutions exist it has been found that close approximations of these solutions can be found. However, analytical methods developed to date by workers such as Rozvany et al. can only be applied to problems for which the optimal layout is independent of loading. Thus the proposed method has also been applied to a range of load dependent problems, for which analytical solutions are currently not available.
Interestingly the new method revealed that one problem in the literature which had been thought to be load independent (providing the load was always applied in a downward direction), is in fact load dependent, rendering the proposed analytical solution less generally applicable than previously thought. Acknowledgements The first author would like to thank the National Science Centre (Poland) for financing the Research Grant no 2015/19/N/ST8/00474 entitled: Topology optimization of thin elastic shells -a method synthesizing shape and free material design. The second and third authors acknowledge the funding provided by the UK Engineering and Physical Sciences Research Council, under grant no. EP/N023471/1.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: Computing extrapolated volumes
As described in Darwich et al. (2010), numerical solutions obtained from numerical layout optimization runs appear to follow a relation of the form: where V n is the numerically computed volume for n equally spaced nodal divisions, V ∞ is the volume when n → ∞, and k and α are constants. Using (9), a weighted least-square approach can be used to find the best-fit values for V ∞ , k and α, with the weighting coefficient taken as n. Numerical solutions are given in Tables 1-3.    By duality arguments one can compute the volume of the optimum grillage (grillage-like slab) from the following equality where u is an out of plane displacement function that solves the dual displacement form and p is a load function. As implied in Rozvany (1972a) the four column slab problem considered in section 3.1.4 enjoys a solution u that is independent of the load p provided the latter is always downwards. The analytical layout given in this paper and presented in Fig. 6a furnishes principal trajectories of curvature field κ associated with u and of principal curvatures being equal to ±1/m p . To facilitate comparison of the volume V exact with numerical results, the load p is assumed to be uniformly distributed. Now, by making use of information on the curvature function, the solution u can be recovered region by region; in addition the function u must be continuous and continuously differentiable in each point of the domain. As the layout is bisymmetrical, 1/4 of the domain is considered; for region partition and coordinate systems see Fig. 14. For sake of simplicity m p is taken as unity.