Form-finding of shells containing both tension and compression using the Airy stress function

Abstract

Pure-compression shells have been the central topic in the form-finding of shells. This study investigated tension-compression mixed-type shells by utilizing a NURBS-based isogeometric form-finding approach that analyses Airy stress functions to expand the possible plan geometry. A complete set of smooth version graphic statics tools is provided to support the analyses. The method is validated using examples with known solutions, and a further example demonstrates the possible forms of shells that the proposed method permits. Additionally, a guideline to configure a proper set of boundary conditions is presented through the lens of asymptotic lines of the stress functions.

Keywords

Force diagram Airy stress function graphic statics form-finding membrane shell tension-compression mixed type shell asymptotic lines wave equation

Introduction

Shell structures, such as concrete, masonry, and metal shells, steel or timber gridshells, cable nets, and fabric structures, are elegant and light weight (Figure 1). Their dominating load-carrying action is membrane action, which is a combination of tensile, compressive, and shear stresses or forces acting in a plane tangential to the surface of the structure. Shell structures containing only tensile or only compressive stresses cannot have unsupported boundaries with overhang but their boundary curves must draw curves shrinking inward (Figure 2).

Figure 1.

Compression shells and lightweight tensile structures: (a) highway service area Deitingen South by Heinz Isler©, 2009, (b) British Museum Great Court Roof by fosters and partners in collaboration with Chris Williams, Andrew Dunn©, 2005, (c) Multihall in Mannheim by Frei Otto, Immanuel Giel©, 2006, and (d) olympic stadium in Munich by Frei Otto, Meister Eiskalt©, 2014.

Figure 2.

Plans (left) and bird’s-eye view (right) of pure compression gridshells that have open edges on the boundary.

To allow for overhangs, a mix of tension and compression stresses is needed (Figure 3). Whenever there are compressive stresses, some bending stiffness is required to avoid buckling. Bending stiffness can also be added to allow for the shape of the shell to deviate from the pure membrane action shape (Figure 4).

Figure 3.

Félix Candela’s hyperbolic-paraboloidal reinforced concrete shell for the Los Manantiales restaurant at Xochimilco.

Figure 4.

Sydney Opera house by Jorn Utzon (left) and JKF Terminal by Eero Saarinen (right).

During the design of shell structures, membrane action is commonly secured using structural form-finding whereby a geometry is determined in such a way that no bending action is needed for the transfer of the dominating load. Most numerical methods for form-finding of shells simulate a physical model that might involve hanging chains or fabric that will be then inverted and form a compressive structure. The models are usually discretized using line elements connected to mutual nodes where external loading also is applied, and the geometry is iteratively updated based on the previously computed stress state. As a result, such methods are restricted to tensile only or compressive only shell structures.

In this study, we revisit, clarify, and develop the numerical form-finding method presented by Miki et al.¹ and the smooth version Graphic Statics tools discussed by Miki et al.² The method takes as input a pre-determined membrane stress state specified using the scalar-valued Airy³ stress function $ϕ$ defined on a plane. With $ϕ$ given, the shape of the shell $z$ is solved such that the vertical equilibrium equation

\frac{\partial^{2} ϕ}{\partial x^{2}} \frac{\partial^{2} z}{\partial y^{2}} - 2 \frac{\partial^{2} ϕ}{\partial x \partial y} \frac{\partial^{2} z}{\partial x \partial y} + \frac{\partial^{2} ϕ}{\partial y^{2}} \frac{\partial^{2} z}{\partial x^{2}} = - ρ,

(1)

is satisfied^4,5 given an appropriate set of boundary conditions. $ρ$ is the upward vertical loading per unit plan area, which may be constant or may depend on the shape. Equation (1) is equivalent to equation (f) on page 462 in Theory of Plates and Shells,⁶ whose authors attribute the first use of the Airy stress function in this manner to Pucher.⁷ In Chapter 10.3 and the subsequent chapters of Theory and Practice of Membrane Shells, Csonka⁵ discusses stress functions for shells on polygonal plans in pure compression and star-shaped plans with a mix of compression and tension. The stress function is reintroduced in the Appendix D.

Some shell structures consist of several smooth surfaces joined via sharp kinks where concentrated forces arise, for example, several of Félix Candela’s concrete shells (Figure 3) and Jorn Utzon’s Sydney Opera house (Figure 4, left). Such force concentrations can be represented in the Airy stress function by narrow regions with a high degree of curvature. If these narrow regions are taken to the limit, “folds” form on the stress function.⁸ Then, the piecewise smooth stress function can be seen as a hybrid between a smooth stress function and a stress polyhedral.⁹

In graphic statics for bar frameworks such as planar trusses, stress polyhedrons can be used to establish a discrete force diagram,^10–14 which has a reciprocal relation to a discrete form diagram. Graphic statics originated in the 18th and 19th century^15–20 and have recently regained attention^21–30 with generalization into higher dimensions^31–38 and applications for the form-finding of shells.^14,39–44 Attempts have been made to establish a similar relation between a piecewise smooth stress function and the shape of the shell,¹ which require the computation of the Christoffel symbols of the second kind. To avoid the often tedious computation of the Christoffel symbol, a reciprocal relation between the piecewise smooth continuous stress function, a pre-computed continuous force diagram, and continuous form diagram (i.e. plan geometry) has been established.²

If $ϕ$ has a positive Gaussian curvature, it represents a state of pure compression or pure tension, whereas if the Gaussian curvature is negative, it represents a state of compression and tension. Thus, the method, in principle, allows for the form-finding of shells with any membrane stress state. The stress functions considered in this paper all have negative Gaussian curvatures, resulting in mixed tension-compression shell structures. For such cases, equation (1) becomes a hyperbolic partial differential equation,⁴ which is often challenging to solve analytically.

In line with the methodology of Miki et al.,¹ equation (1) is solved numerically using isogeometric analysis with NURBS surfaces as finite elements.^45–48 In isogeometric analysis, the computational model and the geometry model are the same, eliminating the need to discretize the smooth shell surface into flat panels or a network of straight-line elements. Without the need to mesh the surface, designers can concentrate on the form in early design stages, leaving discretization of the surface into structural elements such as panels and bars for later stages.

In this study, we only consider isogeometric models consisting of multiple NURBS surfaces, with their control points on the boundary edges shared at their intersections. This strategy makes the computational method simpler because there is no need to consider additional boundary terms to glue surfaces. More general isogeometric models may contain B-rep representations, where each surface is defined by a combination of base surface—usually a NURBS surface—and arbitrary boundary curves and holes. In this case, extra “weak” boundary terms need to be incorporated to stitch surfaces together. Nitsche’s method is known as one of the standard methods to do this.^49–52 However, we leave it to the future work and only discuss simple cases in which each surface is rectangular in its UV-space and adjacent surfaces simply share control points.

Contributions

In this study, we breakdown the original method presented by Miki et al.¹ to make it more accessible and reproducible. While the earlier paper¹ restricted its attention to shells of pure compression, we concentrate on shells containing both compression and tension stresses. The contributions in this paper are as follows:

Provides a complete set of smooth version graphic statics tools through stress function, form diagram, and force diagram (Section 2).

Revisits the numerical form-finding method discussed by Miki et al.¹ (Section 2). Minor revisions include

(a) The stress function and the shell can have different parametric representations.

(b) Linear springs are provided to constrain an arbitrary point to a desired height.

Known solutions of tension-compression mixed type shells are used to validate the method (Section 3).

A guideline to choose proper boundary conditions for tension-compression mixed type problem is addressed (Section 4) through the lens of asymptotic lines.

Introduces an Airy stress function (Section 5) that can take an arbitrary plan geometry. An example problem is provided to demonstrate the possible forms enabled by the proposed stress functions.

Theoretical background and numerical method

In this section, starting from the traditional discrete graphic statics tools, their smooth versions are derived. Then, a NURBS-based isogeometric form-finding method is introduced. The method solves a linear system of equations when the surface area is evaluated on the projected plane. Its nonlinear version that accounts for the surface area accurately is also presented.

Discrete force diagrams and Airy stress polyhedron

Figure 5 shows an example of a force and form diagram pair. If each internal node in the form diagram is balanced by axial forces acting along the incoming lines, a reciprocal force diagram consisting of closed polygons can be constructed. If the form diagram is not in equilibrium, one or several polygons in the force diagram can not be closed. In the force diagram, the lengths of the edges of the polygons represent the magnitude of the corresponding force and are rotated by 90 compared to the corresponding line in the form diagram.⁵³ The reciprocal relation is such that a point, a polygon, and a line in a form diagram maps to a polygon, a point, and a line in the force diagram, respectively.

Figure 5.

An example of form and force diagrams: (a) form diagram, (b) force diagram, and (c) overlay of (a) and (b).

The existence of a force diagram is equivalent to the existence of a planar-faced polyhedron (i.e. stress polyhedron⁹ or Airy stress polyhedron) whose planar projection is the form diagram.^10–14 If the plane of a face in such a polyhedron is expressed as

a x + b y + c = ϕ,

(2)

its corresponding point in the force diagram is located at $(a, b)$ . At the same time, the “normalized” normal vector of the same face is given by

N = (a, b, - 1) .

(3)

Thus, a force diagram can be obtained by flattening all normal vectors of the polyhedron and connecting the points if their corresponding faces are adjacent (Figure 6).

Figure 6.

Stress polyhedron and face normal vectors: (a) form diagram, (b) force diagram, (c) stress polyhedron and normal vectors, and (d) normal vectors.

Continuous force diagram

For a stress function $ϕ = ϕ (x, y)$ , the normalized normal vector is given by

N = (\frac{\partial ϕ}{\partial x}, \frac{\partial ϕ}{\partial y}, - 1),

(4)

and the second Piola-Kirchhoff (PK) stress tensor⁵⁴ for the projected geometry of the shell can be computed as

S = [\begin{matrix} S_{x x} & S_{x y} \\ S_{y x} & S_{y y} \end{matrix}] = [\begin{matrix} \frac{\partial^{2} ϕ}{\partial y^{2}} & - \frac{\partial^{2} ϕ}{\partial y \partial x} \\ - \frac{\partial^{2} ϕ}{\partial x \partial y} & \frac{\partial^{2} ϕ}{\partial x^{2}} \end{matrix}] .

(5)

When partial differentiation is performed to a function multiple times, the order of the differentiation is insignificant. This often-overlooked condition, that is

\frac{\partial^{2}}{\partial y \partial x} = \frac{\partial^{2}}{\partial x \partial y} .

(6)

allows any choice of $ϕ$ automatically satisfies the horizontal equilibrium equations,

\begin{array}{l} \frac{\partial S_{x x}}{\partial x} + \frac{\partial S_{x y}}{\partial y} = \frac{\partial}{\partial x} \frac{\partial^{2} ϕ}{\partial y^{2}} - \frac{\partial}{\partial y} \frac{\partial^{2} ϕ}{\partial x \partial y} & = 0 \\ \frac{\partial S_{x y}}{\partial x} + \frac{\partial S_{y y}}{\partial y} = \frac{\partial}{\partial y} \frac{\partial^{2} ϕ}{\partial x^{2}} - \frac{\partial}{\partial x} \frac{\partial^{2} ϕ}{\partial x \partial y} & = 0 \end{array}} .

(7)

Similarly, we choose any pair of functions $X$ and $Y$ that satisfy

\frac{\partial Y}{\partial x} = \frac{\partial X}{\partial y} \Leftrightarrow \frac{\partial Y}{\partial x} - \frac{\partial X}{\partial y} = 0 .

(8)

For such a pair, the stress tensor is given by

S = [\begin{matrix} S_{x x} & S_{x y} \\ S_{y x} & S_{y y} \end{matrix}] = [\begin{matrix} \frac{\partial Y}{\partial y} & \frac{- \partial X}{\partial y} \\ - \frac{\partial Y}{\partial x} & \frac{\partial X}{\partial x} \end{matrix}] .

(9)

Thus, from equations (5) and (9),

X = \frac{\partial ϕ}{\partial x}, Y = \frac{\partial ϕ}{\partial y},

(10)

which are the components of the normal vector in equation (4). Hence, in analogy with the relationship between the normal vector of a stress polyhedron and the discrete force diagram, $(X, Y)$ is the continuous force diagram corresponding to a continuous stress function.²

Note that the existence of a stress function $ϕ$ is not necessary to construct a force diagram $(X, Y)$ . In fact, it is sufficient to require the conditions given in equation (8). This condition is known as “curl-free” and is essentially an integrability condition of gradient vectors. It is also explained in a structured manner in exterior calculus.⁵⁵

Curl-free means that the deformation gradient

F = [\begin{matrix} F_{x x} & F_{x y} \\ F_{y x} & F_{y y} \end{matrix}] = [\begin{matrix} \frac{\partial X}{\partial x} & \frac{\partial X}{\partial y} \\ \frac{\partial Y}{\partial x} & \frac{\partial Y}{\partial y} \end{matrix}],

(11)

is symmetric and consequently has no components of rotation. The force diagram can be considered as a deformed shape of the form diagram, and the local deformation at each point can be described by a deformation gradient. The deformation gradient maps a unit circle in the form diagram (Figure 7(a)) to an ellipsoid in the force diagram (Figure 7(b)). Without any rotational deformations, there are two orthogonal directions where the local deformation has stretches only (the red lines in the figure) and these coincide with the principal directions of the stress tensor; but the first and the second directions are swapped (Figure 7(b)).

Figure 7.

Rotation-free local deformation: (a) form diagram and (b) force diagram.

Introduction of curvilinear coordinate system

So far, an orthonormal coordinate system $(x, y)$ has been used. However, for shells with arbitrary boundary shape, a curvilinear coordinate system $(u, v)$ is often preferred. The mapping between the two systems can be described well by the terms of differential geometry.⁵⁶ Books such as Ciarlet.^57–59 are a good reference because they discuss differential geometry from the viewpoint of the theory of elasticity.

Let $(u, v) \to (x, y)$ and $(u, v) \to (X, Y)$ be two symmetrical maps. Then, the maps $(u, v) \to (x, y)$ and $(u, v) \to (X, Y)$ define a pair of parametric surfaces representing the form and force diagrams that can be expressed as

x (u, v) = (x (u, v), y (u, v)),

(12)

respectively

X (u, v) = (X (u, v), Y (u, v)) .

(13)

Because the inverse of $(u, v) \to (x, y)$ exists (or we assume singularity/overlapping does not occur in our surfaces), $(x, y) \to (X, Y)$ exits. If it is given by equation (10), it can be computed by

\begin{array}{l} X & = \frac{\partial ϕ}{\partial x} = \frac{\partial ϕ}{\partial u} \frac{\partial u}{\partial x} + \frac{\partial ϕ}{\partial v} \frac{\partial v}{\partial x}, \\ Y & = \frac{\partial ϕ}{\partial y} = \frac{\partial ϕ}{\partial u} \frac{\partial u}{\partial y} + \frac{\partial ϕ}{\partial v} \frac{\partial v}{\partial y}, \end{array}

(14)

where

[\begin{matrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \end{matrix}] = {[\begin{matrix} \frac{\partial x}{\partial u} & \frac{\partial x}{\partial v} \\ \frac{\partial y}{\partial u} & \frac{\partial y}{\partial v} \end{matrix}]}^{- 1} .

(15)

The deformation tensor $F$ on a $(u, v)$ system is known⁶⁰ as

F = G_{u} \otimes g^{u} + G_{v} \otimes g^{v},

(16)

where $g^{u}$ and $g^{v}$ are the dual base vectors of $x (u, v)$ and $G_{u}$ and $G_{v}$ are the base vectors of $X (u, v)$ . For details on the computation of the base vectors and the extraction of the components $F_{x x}$ , $F_{x y} = F_{y x}$ , and $F_{y y}$ , see the Appendix A.

The deformation gradient tensor $F$ given by equation (16) is equivalent to the one given by equation (11), and both are symmetric. Therefore,

G_{u} \cdot g_{v} = G_{v} \cdot g_{u},

(17)

which is an equivalent condition to the curl-free condition. A special case of the curl-free condition is when $G_{u} \cdot g_{v} = G_{v} \cdot g_{u} = 0$ (Figure 8), which is an interpretation of the perpendicular condition of the discrete force diagrams to the continuous ones.

Figure 8.

Perpendicular base vectors: (a) form diagram and (b) force diagram.

In general, a computation of the second derivatives of the stress function on a curvilinear coordinate system requires the computation of the Christoffel symbols of the second kind, $Γ_{i j}^{k}$ , describing the curvature of the isocurves of a parametric surface (see e.g. Green and Zerna,⁴ Eisenhart,⁵⁶ and Ciarlet^57–59). However, through the above formulation, the first differentiation is computed on $(u, v)$ , yielding a vector on the $(u, v)$ system, which is transformed to the $(x, y)$ system using equation (14). Then, a differentiation of $(X, Y)$ with respect to $(u, v)$ is computed using equation (16), but $(X, Y)$ is already a vector in the $(x, y)$ system. Therefore, there is no need to compute the Christoffel symbols. Moreover, the precomputed force diagram, $(X, Y)$ , assists the evaluation of the crease curves in a piecewise smooth stress function (Section 2.7). This is where the proposed method takes advantage of the original method by Miki et al.,¹ substantially simplifying the computations.

Linear method

The horizontal equilibrium of the shell is ensured by equation (7). To complete the equilibrium, the shape of the shell $z$ must be determined so that the vertical equilibrium given by equation (1) is fulfilled. Using the principle of virtual work, a finite number of simultaneous equations can be obtained as an approximation of the vertical equilibrium equation. This approach is equivalent to the Galerkin method,⁶¹ and the details are provided in the Appendix B.

Assume the vertical loading $ρ$ is constant, so that the solution $z$ of equation (1) becomes a linear problem. Later on, in Section 2.5, $ρ$ will be taken as the self-weight of the shell. Then, $ρ$ is no longer constant but depends on the actual area of the shell, turning equation (1) into a non-linear problem.

We introduce a curvilinear coordinate system $(ξ, η)$ on the shell that is different from that of the force and form diagrams $(u, v)$ . Let $h_{ξ}$ and $h_{η}$ be the base vectors and $h_{i j}$ be the first fundamental form of the $(ξ, η)$ system, and we denote the reference states of the properties of the shell that are not affected by differential operators by a bar ⁻. Then, from the principle of virtual work, we get:

δ w = \sum_{s u r f a c e s} \int (\frac{1}{2} {\bar{S}}^{i j} δ h_{i j} - \bar{ρ} δ z) d \bar{a} = 0,

(18)

where

d \bar{a} = \sqrt{\det {\bar{h}}_{i j}} d ξ d η,

(19)

is a small area element, ${\bar{S}}^{i j}$ the components of the second PK stress tensor re-evaluated on $(ξ, η)$ , and $\bar{ρ}$ a constant loading per unit area, all measured in the “plan” of the form diagram. The Einstein summation convention applies, and the indices $i$ and $j$ are either 1 or 2, where 1 corresponds to the $ξ$ -direction and 2 corresponds to the $η$ -direction. Details on the extraction of the stress components on the $(ξ, η)$ system are given in the Appendix A.

If the shell is represented by multiple finite elements or parametric surfaces (e.g. triangles, quadrilaterals, or NURBS-surfaces) defined by n independent control points, all functions defined on the shell, such as $h_{i j}$ and $z$ , become functions of the control point coordinates. By packing the coordinates of these control points in a 3n-dimensional column vector $r$ and replacing the variation of functions $δ f$ by $\nabla f δ r$ , where $δ r$ is a 3n-dimensional arbitrary column vector and the gradient operator is

\nabla = [\begin{matrix} \frac{\partial}{\partial r_{1}}, \dots, \frac{\partial}{\partial r_{3 n}} \end{matrix}],

(20)

equation (18) becomes

\sum_{s u r f a c e s} \int (\frac{1}{2} {\bar{S}}^{i j} \nabla h_{i j} - \bar{ρ} \nabla z) d \bar{a} = 0 .

(21)

Since the horizontal equilibrium is automatically ensured, solving equation (21) in terms of $z$ is enough to obtain a complete equilibrium. In case equation (21) is solved in terms of $x$ and $y$ , the plan of the solution should match the form diagram.

Note that equation (21) is linear in $r$ . By distributing integration points and performing a numerical integration (see Section 2.8), a system of linear equations is obtained. With a sufficient number of fixed points given, this system of equations can be solved easily by computing an inverse matrix. Although presented quite differently, Pauletti and Pimenta⁶² have already pointed out that having a constant second PK stress tensor results in a linear system of equations.

Nonlinear iterative method

If $\bar{ρ}$ is the self-weight of the shell, equation (18) only gives an approximate solution. The exact solution is given by

δ w = \sum_{s u r f a c e s} \int (\frac{1}{2} {\bar{S}}^{i j} δ h_{i j} d \bar{a} - \bar{ρ} δ z d a) = 0,

(22)

where

d a = \sqrt{\det h_{i j}} d ξ d η

(23)

is the actual element area measured on the shell. Since $d a$ depends on $h_{i j}$ which is a property of the sought solution $z$ , equation (22) is a non-linear problem requiring an iterative solution process.

In each iterative step, the residual forces of the system are computed by

R = \sum_{s u r f a c e s} \int r^{T} (\frac{1}{2} {\bar{S}}^{i j} \nabla^{2} h_{i j} d \bar{a} - \bar{ρ} \nabla z d a),

(24)

where $\nabla^{2} = \nabla^{T} \nabla$ and $\nabla^{2} h_{i j}$ is a constant matrix. Thus, the matrix $S^{i j} \nabla^{2} h_{i j}$ does not change during the iterations.

Then, with $h^{i j}$ , the inverse of $h_{i j}$ , and using the relation

d \sqrt{\det h_{i j}} = \frac{1}{2} h^{i j} d h_{i j} \sqrt{\det h_{i j}},

(25)

where $d$ is a generic differential operator that can be either a $δ, d, o r \nabla$ operator, the system of equations to solve in each step becomes

\begin{array}{l} \frac{1}{2} [\sum_{s u r f a c e s} \int ({\bar{S}}^{i j} \nabla^{2} h_{i j} d \bar{a} - \bar{ρ} h^{i j} \nabla z^{T} \nabla h_{i j} d a)] Δ r \\ = - R . \end{array}

(26)

Note that in the first step, solving equation (26) is the same as solving equation (21). The convergence of the solution is often improved when the second term in the parentheses of equation (26) is omitted, that is, it is equivalent to repeating the linear method but with the updated loading. Another good strategy is to scale down $δ r$ by $0.5$ , which contributes to the stability of the method.

Point and area supports

A control point can be easily fixed to a specific height by eliminating the control point parameters completely from the system of equations. However, if NURBS surfaces are used as finite elements, the control points do not pass through the surface except at the corner points. Thus, to constrain an arbitrary point of such a surface to a specific height, a linear spring can be used.

The elastic energy of such a spring is

\frac{1}{2} k {(z (r) - \bar{z})}^{2},

(27)

where $\bar{z}$ is the prescribed height of the point, $z (r)$ is a linear function that converts the control points to the height of the point, and $k$ is a penalty factor.

When such springs are incorporated into the system, the equilibrium of the entire system is

\sum_{s p r i n g s} k (r^{T} \nabla z^{T} \nabla z - 2 \bar{z} \nabla z) + m e m b r a n e t e r m = 0,

(28)

where the spring term is linear in $r$ and the membrane term is discussed in Section 2.4 or Section 2.5 depending on the loading condition.

Similarly, area supports defined by closed curves can be imposed on the shell by the use of a set of linear springs distributed along the curve. This approach is used later on in Section 5 when discussing some examples.

Piecewise smooth stress functions

Whenever the shell is made up of several smooth surfaces joined by sharp kinks along their intersection curves, the stress function must be a piecewise smooth function. For such shells, the equilibrium is sustained by concentrated forces acting along the kinks on the shell. The axial force along a kink in the shell surface can be calculated by measuring the “jump” of the normal vectors on either side of the corresponding kink in the stress function. The same entity can be found by measuring the “jump” in the force diagram. Denoting the “jump” by $\bar{L}$ , the second PK stress tensor along a kink in the surface¹ becomes

{\bar{n}}^{11} = ({\bar{h}}_{11})^{- 1} \bar{L},

(29)

and the principle of virtual work adding up to equation (18) or equation (22) is

δ w = \sum_{c u r v e s} (\int \frac{1}{2} {\bar{n}}^{11} δ h_{11} d \bar{ℓ}),

(30)

where

d \bar{ℓ} = \sqrt{\det {\bar{h}}_{11}}

(31)

is a length element on the intersection curve as seen in the plane of the form diagram.

Numerical integration

In this work, a standard Gauss quadrature was used to numerically compute the integrals over the surfaces and curves.⁴⁵ Because the stress components are calculated from a stress function that is not a polynomial in general, a sufficient number of integration points is hard to estimate. Hence, the number of integration points was controlled by an integration multiplier, denoted by $N I$ .

We denote the number of control points in the $u$ and $v$ directions by $N U$ and $N V$ , respectively, and the degree of the NURBS surface by $D i m$ . Then, a NURBS surface is typically divided into $(N U - D i m) \times (N V - D i m)$ smaller patches, and each of them is dependent on $(D i m + 1) \times (D i m + 1)$ control points.⁶³

In this study, $(D i m \times D i m \times {(N I)}^{2})$ Gauss integration points were distributed on each patch (i.e. $n = D i m \times N I$ in the commonly used Gaussian quadrature rule tabulation is used). NURBS surfaces with $D i m = 3$ and $N I = 1 o r 2$ work for most cases, resulting in $n = 3 o r 6$ for the Gauss quadrature rule. In the worst case $n = 4$ is enough.

Convergence study

A variational operator $δ$ has infinite degrees of freedom (DoF) by definition. The Galerkin method restricts the variational operator to a finite number of DoF equal to the DoF of the finite element model. If a solution exists in the original PDE, by increasing the DoF in the system, the numerical solution should converge to the exact solution.

In this study, the convergence towards the exact solution is evaluated by studying the total energy of the system given by

E = \sum_{s u r f a c e s} \frac{1}{2} \int {\bar{S}}^{i j} h_{i j} d \bar{a} + \sum_{c u r v e s} \frac{1}{2} \int {\bar{n}}^{11} h_{11} d \bar{ℓ} .

(32)

The elastic energy of the springs is omitted because they will have high energy if used as hard constraints.

It is known that there are incompatible boundary conditions (BCs) that make hyperbolic PDEs ill-posed (i.e. solutions do not exist). In such cases, numerical solutions hardly ever converge. Since a mix of tension and compression stresses results in a hyperbolic PDE, attention must be paid to such situations since the proposed numerical method returns a solution regardless of the problem is well-posed or not. A guideline to choose proper BCs is provided in Section 4.

Benchmark problems

In this section, we present two simple benchmark problems. We only consider Dirichlet type BCs and BCs involving higher order derivatives will not be discussed. With more complex problems, by taking the advantage of the fact that the NURBS surfaces have higher order non-zero derivatives, methods that require computations of higher order derivatives, such as the least squares method, may be considered. However, we leave it to future work.

Hypar

Consider a simple stress function

ϕ = - \frac{1}{2} x^{2} + \frac{1}{2} y^{2},

(33)

that is defined on a square that spans $(- 1, - 1)$ , $(1, - 1)$ , $(1, 1)$ , and $(- 1, 1)$ . In Figure 9, the purple surface corresponds to $ϕ$ , and the gray flat planes are the stress-free region outside the structure apart from the reaction forces seen as the kinks at the top.

Figure 9.

Stress function (purple) and corresponding shell (green) of the hypar problem.

The force diagram of $ϕ$ is given by $(X, Y) = (- x, y)$ , which is curl-free. Figure 10 depicts the force diagram and the “jumps” of the normal vectors between the smooth pink surface of $ϕ$ and the surrounding flat gray planes. By adding edge beams, the shell can be supported by only two points, and the magnitudes of the “jumps” represent the axial forces in the beams.

Figure 10.

Force diagram and the “jumps” of the hypar problem.

Substituting equation (33) into equation (1) gives

- \frac{\partial^{2} z}{\partial y^{2}} + \frac{\partial^{2} z}{\partial x^{2}} = - ρ,

(34)

for which there exists a known solution,

z = \frac{ρ}{4} (x^{2} - y^{2}) .

(35)

The solution satisfies the boundary equilibrium between the edge beams and the shell.

Note that equation (35) is only a particular solution. With more complex BCs, the solution is a sum of the particular solution and a general solution (a solution of equation (34) where the right-hand side is set to zero), which satisfies the BCs. While a general solution always exists in an elliptical problem, in a hyperbolic problem, there exist incompatible BCs, where general solutions do not exist. In general, a hyperbolic PDE only accepts compatible BCs, otherwise the problem becomes ill-posed. This issue is further addressed in Section 4.

This problem can be used as a useful benchmark problem, and the proposed method was tested. The green surface in Figure 9 represents the solution obtained using the proposed method with $ρ = 4.0$ as load. The height of the obtained green shell is $1 / 4$ of $ρ$ , which matches equation (35).

The blue surface in Figure 11 is a plot of the left-hand side of equation (1) multiplied by $1 / ρ$ obtained numerically with a NURBS surface with $32 \times 32$ control points, $D i m = 3$ , and an integration point multiplier $N I = 1$ . It clearly shows a drop of $- 1$ from the orange surface representing zero loading. This indicates that the method provides the sought solution.

Figure 11.

Left-hand side of vertical equilibrium obtained using $32 \times 32$ control points (blue surface) and zero-loading reference (orange surface) for the hypar problem.

Figure 12 shows a convergence analysis of this example. As shown in Figure 12, the total energy remains flat even if the number of control points was increased. This is because the hypar can be accurately represented with even a NURBS surface of $2 \times 2$ control points and with $D i m = 3$ .

Figure 12.

Convergence analysis of the hypar problem.

Bowl

Consider a stress function $ϕ$ and a shell $z$ that both have a radial symmetry defined in a polar coordinate system. Then, the vertical equilibrium becomes

\frac{1}{r} (\frac{\partial ϕ}{\partial r} \frac{\partial^{2} z}{\partial r^{2}} + \frac{\partial^{2} ϕ}{\partial r^{2}} \frac{\partial z}{\partial r}) = - ρ,

(36)

where $r \geq C \geq 0$ is the parameter in the radial direction in the polar coordinate system.

By assuming the same shape for the shell and stress function, so that $ϕ (r) = A z (r)$ where $A$ is a constant, the solution to equation (36) is given by

z (r) = \frac{1}{2 \sqrt{2 A}} (r \sqrt{B^{2} - r^{2}} + \arcsin (r)),

(37)

with $r \leq B \leq 1$ . Thus, at the circumference where $r = B$ , there are concentrated stresses parallel to the circumference direction. If vertical support is provided at $r = C$ , the edge at $r = B$ can “float” without any support and the circumferential concentrated stress is tensile.

Let $A = 1 / 2$ , $B = 1$ , $C = 0.2$ , and $ρ = 1$ . Then, the total height of the shell is given from the analytical solution as $z (B) - z (C) = z (1.0) - z (0.2) = 0.5867$ , and the cross section of the shell is plotted in Figure 13. Solving $z$ numerically with the proposed method gives the shape shown in Figure 14 whose profile and total height, $0.58$ , match the analytical solution.

Figure 13.

Cross section of an analytical solution of the bowl problem.

Figure 14.

Elevation of the numerical solution of the bowl problem.

The blue surface in Figure 15 is a plot of the left-hand side of equation (36) multiplied by $1 / ρ$ obtained numerically with a NURBS surface with $32 \times 32$ control points, $D i m = 3$ , and integration point multiplier $N I = 1$ . Similarly as in the hypar example, the blue surface sits at a level $- 1$ below the red surface representing zero loading, again indicating that the method provides the sought solution.

Figure 15.

Left-hand side of vertical equilibrium obtained using $120 \times 60$ control points (blue surface) and zero-loading reference (orange surface) for the bowl problem.

A convergence analysis was also conducted. Figure 16 shows plots of the energy as a function of the number of control points, proving the solution converges to a single solution.

Figure 16.

The convergence of the total energy of the bowl problem.

Unlike pure-compression or pure-tension zones, there is a restriction on the BCs in tension-compression mixed zones. In this example, because we assumed a radial symmetry when deriving the stress function, fixing the rim of the oculus along a flat circle should work, that is, it is compatible to the problem. This compatible boundary condition may not be the unique one and there may be more, but, in general, an arbitrarily determined boundary condition is incompatible, making the problem ill-posed.

As such, in general, a hyperbolic partial differential equation is well-posed only when a proper set of compatible boundary conditions is given. This makes solving the hyperbolic problems substantially different from solving elliptic problems, that is, pure-compression or pure-tension ones. The next section further elaborates the topic about compatible and incompatible BCs of hyperbolic PDEs.

Asymptotic lines and boundary conditions

Within areas of pure compression or tension, support BCs cause singularities that dissipate rapidly and only affect the solution locally, allowing free placement of the supports in such areas. However, in areas with a mix of compression and tension, BCs give rise to singularities that are transferred through characteristic lines with no dissipation, often making the hyperbolic PDE ill-posed.^5,64 Such incompatible BCs should be avoided.

The four green surfaces in Figure 17(a) are the numerical solutions obtained by imposing two compatible sets of BCs and two incompatible ones on the hypar example introduced in Section 3.1. As in Figures 11 and 15, the blue surfaces of Figure 17(a) are supposed to drop by 1.0 below the reference surface if the PDE is well-posed. However, the incompatible BCs cause disturbances seen as bumpy ridges that run across the solution, and the blue surfaces do not drop by 1.0 but form spikes.

Figure 17.

Compatible and incompatible BCs: BCs (red), numerical solution (green), and left-hand-side of the equilibrium multiplied by $1 / ρ$ (blue), which should drop by 1.0 below the references surface (gray) for proper solutions: (a) hypar example of Section 3.1, (b) bowl example of Section 3.2, and (c) Félix Candela example briefly mentioned in Appendix D.

A second-order hyperbolic PDE is often called a wave equation, and waves propagate along so-called characteristic lines. Thus, a perturbance transfers from one end of a characteristic line to the other end with no dissipation, and if no compatible BC exists at the other end, the problem becomes ill-posed (i.e. no solution exists). Even if solutions do not exist, the proposed method returns a numerical solution with no errors. Hence, extra attention must be paid to this issue.

In the problems discussed in this paper, the characteristic lines are the projections of the asymptotic lines of the stress function to the ground plane (see Chapter 4 in Csonka,⁵ Chapters 1 and 2 in Sanchez-Palencia et al.,⁶⁵ and Appendix C). An asymptotic line is a line on a surface whose normal curvature is zero, and it only exists in areas where the Gaussian curvature is negative. Typically, the asymptotic lines are a group of diagonal lines that run across the negative Gaussian curvature area, intersecting the principle curvature lines at roughly 45 degrees.

The general solution of the hypar problem is $z (x, y) = F (x + y) + G (x - y)$ with $F$ and $G$ arbitrary functions that are determined such that boundary conditions are satisfied, and the asymptotic lines are such that they result in characteristic lines given by $x - y = c o n s t .$ and $x + y = c o n s t .$ , which run parallel with the edges of the shell. As seen in Figure 17(a), the BCs that contradict the general solution are incompatible or, in other words, BCs that do not have compatible BCs at the end of the characteristic lines are incompatible. Figure 17(b) and (c) show more examples of compatible and incompatible BCs imposed on the bowl example of Section 3.2 respectively on a problem reproducing Félix Candela’s Xotimilco restaurant (see Appendix D). The BCs for these problems were guessed based on the knowledge of radial symmetry. However, in general, identifying compatible BCs for hyperbolic PDEs is a challenge.

To avoid this challenge, “safe” areas with positive Gaussian curvature may be added to the stress function in which the supports are placed. This strategy is illustrated in the examples in Figure 18 where the purple surfaces are the used Airy stress functions with their asymptotic lines drawn in white. Though not shown, the “blue surfaces” dropped properly for all cases. The safe-area-strategy is used in the examples discussed in Section 5.

Figure 18.

Moving supports away from the asymptotic line zones is recommended to avoid compatibility issue of BCs. (Green): shells. (Purple): stress functions.

(a) Hypar example of Section 3.1.

(b) Bowl example of Section 3.2.

More complex example

A stress function that can take a general plan geometry

Given a closed flat guide curve defined by $(c_{x} (u), c_{y} (u))$ , a simple pair of stress function $ϕ$ and form diagram $x$ may be obtained as

\begin{array}{l} ϕ (u, v) = Q (v) \\ x (u, v) = (v c_{x} (u), v c_{y} (u)) \end{array}},

(38)

where $v$ is bounded by $v_{\min} > 0$ and $v_{\max}$ , and $Q (v)$ defines the generatrix profile curve of the stress function. Then, from equation (10), the force diagram follows as

X = Q^{'} (v) \frac{1}{c_{x} c_{y^{'}} - c_{y} c_{x^{'}}} (c_{y^{'}}, - c_{x^{'}}),

(39)

where a prime is used to denote differentiation.

If $Q^{'} (v) = 0$ at an edge of the shell, the stress function tangentially touches the horizontal plane intersecting the edge of the stress function. Such points on the form diagram $x$ are mapped to the center of the force diagram $X$ . Therefore, the distances between any proximity points of the edge as measured in the force diagram are always zero. Hence, there is no horizontal thrust acting on the edge. In contrast, if $Q^{'} (v) \neq 0$ at an edge, the stress function forms a kink with the horizontal plane intersecting the edge.

As pointed out before, for example, when discussing the kinks between the hypar stress function and the surrounding planes in Figure 10, a kink represents concentrated stresses acting along the edge, which are balanced by stresses acting in the perpendicular direction. The magnitude of the axial force can be obtained by measuring the “jump” from the corresponding points of the edge in the force diagram to its center, which represents the axial force of the reinforcing beam placed along the perimeter of the edge.

The proposed stress function, form diagram, and force diagram of equations (38) and (39) are valid not only for smooth guide curves but also for piecewise smooth guide curves. In those cases, the stress function also becomes a piecewise smooth function, with creases along lines intersecting the discontinuities in the guide curve. As discussed in Section 2.7, such stress functions result in “jumps” in the force diagram representing concentrated axial forces running along the corresponding creases in the shell, and their magnitudes are given by equation (29).

Figure 19 shows a form diagram generated using a piecewise smooth guide curve with $v_{\min} = 0.2$ and $v_{\max} = 1.0$ . Vertical supports are provided along three short stretches of the exterior perimeter of the form diagram. Along these stretches, the form diagram is modified with small semi-circles that ensure a local area of pure compression close to the support avoiding problems of incompatible BCs.

Let the stress function $ϕ = Q (v)$ be such that

Q^{'} (v) = (v - 1) + 12 λ {(v - 1)}^{2},

(40)

where $λ$ is a parameter embedded to the equation to control shape of the guide curve (generatrix). Then, a stress function that comprises three smooth areas and three crease curves is obtained, as shown on the left in Figure 20, resulting in a force diagram shown to the right. As illustrated in Figure 20, if one travels from point a to b to c on the stress function, corresponding travel from A to B to C takes place on the force diagram where point B exists on both sides of the “jump,” representing the crease running through point b on the stress function.

Figure 19.

Dimensions of the used form diagram (unit = [m]).

Figure 20.

A piecewise-smooth stress function (left) and the corresponding force diagram (right).

Figure 21 shows two numerical form-finding results obtained using the discussed functions $ϕ$ , $x$ , and $X$ , one using the linear method (Section 2.4) and the other the nonlinear method (Section 2.5). Each leaf is represented by 2 NUBRS surfaces with $32 \times 16$ control points, $D i m = 3$ and $N I = 2$ , and $λ$ of equation (40) was set to 0.

Figure 21.

Solutions using the linear method (left) and the nonlinear method (right), both with $ρ = - 0.05$ .

The resulting shells are rather similar, with the nonlinear solution a little higher than the linear solution. The difference is due to the accurate account of the surface area in the nonlinear method. The shells were obtained with loading $ρ = - 0.05$ , and the nonlinear method converged, as proved by the history of the residual norm shown in Figure 22. However, when the load was increased to $ρ = - 0.15$ , the nonlinear method did not converge. As such, the nonlinear method does not converge when the loading is large. However, fortunately, in the cases we tested, although it depends on the problems, the height of the shell at the maximum loading $ρ$ that gave a convergence gives forms that are usable as architectural spaces.

Figure 22.

History of the residual norm of the vertical equilibrium equation solved with nonlinear method and $ρ = - 0.05$ .

Verification

A convergence study was conducted by increasing the number of control points from $20 \times 10$ to an extreme of $400 \times 200$ . $D i m = 3$ NURBS surfaces were used as finite elements, and the integration point multiplier was set to $N I = 2$ . As shown in Figure 23, the solutions were very stable even if the number of control points was increased. Figure 24 is a plot of the energy versus the number of control points, and the solutions converged.

Figure 23.

From left to right: solutions obtained with $40 \times 20$ , $160 \times 80$ , and $360 \times 180$ control points.

Figure 24.

Convergence study, (blue) total energy, and (orange) increments of the energy from the last step .

So far, in the form-finding process, only an ideal membrane stress state is considered, and thus, structural properties such as thickness, stiffness, and deformation are not studied.

Using the shell obtained from the nonlinear method, a series of deformation analyses was performed by varying the thickness of the shell from a realistic one to an extremely small one. The three crease curves where modeled as solid beams with cross-section width × height = 0.2 m × 0.3 m, and standard steel was used for both the shell and beams. The analysis was performed using NURBS-based shell and beam elements provided by Kiwi3D,⁶⁶ which consider both membrane and bending stiffness. The NURBS surfaces and curves obtained in the form-finding stage was used as finite elements without meshing.

Table 1 shows the maximum deflection with varied thicknesses. Even for a small thicknesses, the deformations are kept rather low. This indicates that the shell is working primarily in membrane action, just as intended.

Table 1.

Maximum deflection with varied thickness.

Thickness $t$ (mm)	Max deflection $δ$ (mm)	Deflection/cantilever depth $δ / L$ (%)
0.1	76	0.54
0.5	80	0.57
1	82	0.59
5	83	0.60
10	81	0.58
50	62	0.44

Percentage based on cantilever depth L = 13,860 mm.

Another way to check if the shell works in pure membrane action is to study the ratio between the membrane energy and the bending energy of the shell. As can be seen in Figure 25, the membrane energy becomes more dominant as the thickness of the shell is decreased.

Figure 25.

Ratio between membrane energy and bending energy for thicknesses 10, 5, 1, and 0.1 mm. Blue = 100% membrane energy and red = 100% bending energy.

Figure 26 displays the left hand side of the vertical equilibrium multiplied with $1 / ρ$ obtained for one third of the shell with $240 \times 120$ control points, $D i m = 3$ and $N I = 2$ . As expected, the blue surface drops by $1.0$ , except at two areas which correspond to the locations of the supports.

Figure 26.

A blue surface obtained with $240 \times 120$ control points, $D i m = 3$ , and $N I = 2$ .

From these observations, it can be concluded that the obtained solution is an approximation of the original equilibrium equation and that the shell works by membrane action.

Remarkably, few asymptotic lines of the stress function (Figure 27, left) emerge as waves running on the form-finding result (Figure 27, right). Compared to the bumpy ridges in the numerical solutions obtained with incompatible BCs discussed in section 4, these waves look clean, and there are no spiky noises in the blue surface except the areas near the supports. Thus, these waves can be distinguished from the numerical errors caused by the incompatible BCs; rather, it can be considered that a solution of a wave equation is obtained.

Figure 27.

An overlay of asymptotic lines (left) to the obtained solution (right).

Implementation scenario

Figure 28 shows an architectural rendering that illustrates an implementation scenario of the example.

Figure 28.

An architectural implementation scenario of the example discussed in this section.

Conclusion

The numerical form-finding of membrane action shells containing a mix of tensile and compressive stresses involves solving a hyperbolic second-order partial differential equation. In this study, an existing NURBS-based isogeometric approach originally designed for the form-finding of pure compression or pure tension membrane shells was developed further. It is demonstrated that the developed method can be used to solve tension-compression mixed type hyperbolic equilibrium equations with a satisfactory accuracy.

The method takes as input a form diagram and a stress function from which a force diagram is computed, and their reciprocal relations are discussed.

In areas where the stress function describes a state of both tension and compression, it has negative Gaussian curvature and thus asymptotic lines, and it was recommended to put supports away from such areas. Otherwise, the boundary conditions must be compatible with the equilibrium equation to make the problem well-posed.

The proposed method was tested using simple problems whose analytical solutions are known, and the numerical and analytical results are shown to have a good agreement. A solution to a more complex problem was also demonstrated and verified, and the formulation of the problem can easily be modified to suit a wide range of architectural applications.

Footnotes

Appendix

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by Skidmore, Owings, & Merrill from 2015 to 2020.

ORCID iDs

Masaaki Miki

Alexander Sehlström

References

Miki

Igarashi

Block

. Parametric self-supporting surfaces via direct computation of Airy stress functions. ACM Trans Graph 2015; 34(4): 1–12.

Miki

Mitchell

Baker

. Continuum form and force diagrams in NURBS-based shell form-finding. In: Lázaro

Bletzinger

Oñate

(eds) Form and force 2019. Madrid, Spain: International Association for Shell and Spatial Structures (IASS), 2019, pp.1–8.

Airy

GB.

IV . On the strains in the interior of beams. Philos Trans R Soc Lond 1863; 153: 49–79.

Green

Zerna

. Theoretical elasticity. 2nd ed. Oxford, England: Oxford University Press, 1968.

Csonka

. Theory and practice of membrane shells. Duesseldorf, Germany: VDI-Verlag., 1987.

Timoshenko

Woinowsky-Krieger

. Theory of plates and shells. New York, NY: McGraw-hill, 1959.

Pucher

. Die berechnung der dehnungsspannungen von rotationsschalen mit hilfe von spannungsfunktionen. Int Vereinig Brücken-u Hochbau Abh 1938; 5: 275.

Williams

McRobie

. Graphic statics using discontinuous airy stress functions. Int J Space Struct 2016; 31(2–4): 121–134.

Klein

Wieghardt

. Über spannungsflächen und reziproke diagramme mit besonderer berücksichtiging der maxwellschen arbeiten. Arch Math Phys III/VIII 1904; 8(3): 95–119.

10.

Crapo

Whiteley

. Statics of frameworks and motions of panel structures: a projective geometric introduction. Struct Topol 1982; 6: 43–82.

11.

Fraternali

Angelillo

Fortunato

. A lumped stress method for plane elastic problems and the discrete-continuum approximation. Int J Solids Struct 2002; 39(25): 6211–6240.

12.

Fraternali

. A thrust network approach to the equilibrium problem of unreinforced masonry vaults via polyhedral stress functions. Mech Res Commun 2010; 37(2): 198–204.

13.

Fraternali

Carpentieri

. On the correspondence between 2D force networks and polyhedral stress functions. Int J Space Struct 2014; 29(3): 145–159.

14.

Vouga

Höbinger

Wallner

, et al. Design of self-supporting surfaces. ACM Trans Graph 2012; 31(4): 1–11.

15.

Varignon

. Nouvelle Mecanique Ou Statique: Dont Le Projet Fut Donné En M. DC. LXXXVII. Vol. 2. Paris, France: Jombert, 1725.

16.

Maxwell

. Xlv. on reciprocal figures and diagrams of forces. Lond Edinb Dublin Philos Mag J Sci 1864; 27(182): 250–261.

17.

Culmann

. Die graphische statik. Vol. 1. Zürich, Switzerland: Meyer & Zeller, 1875.

18.

Rankine

. Xvii. principle of the equilibrium of polyhedral frames. Lond Edinb Dublin Philos Mag J Sci 1864; 27(180): 92.

19.

Cremona

. Le figure reciproche nella grafica statica. Tipografia Giuseppe Bernardoni, Milano 1872.

20.

Cremona

. Le figure reciproche nella statica grafica. Milan, Italy: Hoepli, 1879.

21.

Micheletti

. On generalized reciprocal diagrams for self-stressed frameworks. Int J Space Struct 2008; 23(3): 153–166.

22.

Zanni

Pennock

. A unified graphical approach to the static analysis of axially loaded structures. Mech Mach Theory 2009; 44(12): 2187–2203.

23.

Van Mele

Block

. Algebraic graph statics. Comput Aided Des 2014; 53: 104–116.

24.

Lee

Mueller

Fivet

. Automatic generation of diverse equilibrium structures through shape grammars and graphic statics. Int J Space Struct 2016; 31(2–4): 147–164.

25.

Lee

Van Mele

Block

. Disjointed force polyhedra. Comput Aided Des 2018; 99: 11–28.

26.

Mazurek

Beghini

Carrion

, et al. Minimum weight layouts of spanning structures obtained using graphic statics. Int J Space Struct 2016; 31(2–4): 112–120.

27.

Mitchell

Baker

McRobie

, et al. Mechanisms and states of self-stress of planar trusses using graphic statics, part I: the fundamental theorem of linear algebra and the airy stress function. Int J Space Struct 2016; 31(2-4): 85–101.

28.

Alic

Åkesson

. Bi-directional algebraic graphic statics. Comput Aided Des 2017; 93: 26–37.

29.

Baranyai

. On the admissible dualities in Maxwell’s reciprocal construction. Int J Solids Struct 2019; 171: 10–16.

30.

McRobie

. The geometry of structural equilibrium. Royal Soc Open Sci 2017; 4(3): 160759.

31.

Akbarzadeh

Van Mele

Block

. On the equilibrium of funicular polyhedral frames and convex polyhedral force diagrams. Comput Aided Des 2015; 63: 118–128.

32.

Fivet

. Projective transformations of structural equilibrium. Int J Space Struct 2016; 31(2–4): 135–146.

33.

Ole Ohlbrock

Schwartz

. Combinatorial equilibrium modeling. Int J Space Struct 2016; 31(2–4): 177–189.

34.

Ohlbrock

D’Acunto

. A computer-aided approach to equilibrium design based on graphic statics and combinatorial variations. Comput Aided Des 2020; 121: 102802.

35.

McRobie

. Maxwell and Rankine reciprocal diagrams via Minkowski sums for two-dimensional and three-dimensional trusses under load. Int J Space Struct 2016; 31(2–4): 203–216.

36.

Konstantatou

D’Acunto

McRobie

. Polarities in structural analysis and design: n-dimensional graphic statics and structural transformations. Int J Solids Struct 2018; 152: 272–293.

37.

Bolhassani

Akbarzadeh

Mahnia

, et al. On structural behavior of a funicular concrete polyhedral frame designed by 3D graphic statics. Structures 2018; 14: 56–68.

38.

D’acunto

Jasienski

Ohlbrock

, et al. Vector-based 3D graphic statics: a framework for the design of spatial structures based on the relation between form and forces. Int J Solids Struct 2019; 167: 58–70.

39.

Block

Ochsendorf

. Thrust network analysis: a new methodology for three-dimensional equilibrium. J Int Assoc Shell Spatial Struct 2007; 48(3): 167–173.

40.

Liu

Pan

Snyder

, et al. Computing self-supporting surfaces by regular triangulation. ACM Trans Graph 2013; 32(4): 1–10.

41.

Rippmann

Lachauer

Block

. Interactive vault design. Int J Space Struct 2012; 27(4): 219–230.

42.

Rippmann

Block

. Funicular funnel shells. In Rethinking prototype. Proceedings of the design modelling symposium, Berlin, pp.75–89.

43.

De Goes

Alliez

Owhadi

, et al. On the equilibrium of simplicial masonry structures. ACM Trans Graph 2013; 32(4): 1–10.

44.

Panozzo

Block

Sorkine-Hornung

. Designing unreinforced masonry models. ACM Trans Graph 2013; 32(4): 1–12.

45.

Hughes

Cottrell

Bazilevs

. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005; 194(39–41): 4135–4195.

46.

Kiendl

Bletzinger

Linhard

, et al. Isogeometric shell analysis with Kirchhoff–Love elements. Comput Methods Appl Mech Eng 2009; 198(49–52): 3902–3914.

47.

Breitenberger

Apostolatos

Philipp

, et al. Analysis in computer aided design: nonlinear isogeometric B-Rep analysis of shell structures. Comput Methods Appl Mech Eng 2015; 284: 401–457.

48.

Bauer

Wüchner

Bletzinger

. Innovative CAD-integrated isogeometric simulation of sliding edge cables in lightweight structures. J Int Assoc Shell Spatial Struct 2018; 59(4): 251–258.

49.

Nguyen

Kerfriden

Brino

, et al. Nitsche’s method for two and three dimensional nurbs patch coupling. Comput Mech 2014; 53(6): 1163–1182.

50.

Guo

Ruess

Schillinger

. A parameter-free variational coupling approach for trimmed isogeometric thin shells. Comput Mech 2017; 59(4): 693–715.

51.

Apostolatos

Schmidt

Wüchner

, et al. A nitsche-type formulation and comparison of the most common domain decomposition methods in isogeometric analysis. Int J Num Methods Eng 2014; 97(7): 473–504.

52.

Nitsche

. Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind. Abh Math Semin Univ Hambg 1971; 36: 9–15.

53.

Allen

Zalewski

. Form and forces: designing efficient, expressive structures. Hoboken, NJ: John Wiley & Sons, 2009.

54.

Bonet

Wood

. Nonlinear continuum mechanics for finite element analysis. Cambridge: Cambridge University Press, 1997.

55.

Cartan

. Differential forms. Chelmsford, MA: Courier Corporation, 2006.

56.

Eisenhart

. Riemannian geometry. Princeton, NJ: Princeton University Press, Vol. 51, 1997.

57.

Ciarlet

. Mathematical elasticity: volume I: three-dimensional elasticity. Amsterdam, The Netherlands: North-Holland, 1988.

58.

Ciarlet

. Mathematical elasticity: volume II: theory of plates. Amsterdam, The Netherlands: Elsevier, 1997.

59.

Ciarlet

. An introduction to differential geometry with applications to elasticity. J Elast 2005; 78(1–3): 1–215.

60.

Bletzinger

Wüchner

Daoud

, et al. Computational methods for form finding and optimization of shells and membranes. Comput Methods Appl Mech Eng 2005; 194(30–33): 3438–3452.

61.

Connor

Brebbia

. Finite element techniques for fluid flow. Oxford, England: Newnes, 2013.

62.

Pauletti

Pimenta

. The natural force density method for the shape finding of taut structures. Comput Methods Appl Mech Eng 2008; 197(49–50): 4419–4428.

63.

Rogers

. An introduction to NURBS: with historical perspective. Amsterdam, The Netherlands: Elsevier, 2000.

64.

Colton

. Partial differential equations: an introduction. Chelmsford, MA: Courier Corporation, 2004.

65.

Sanchez-Palencia

Millet

Béchet

. Singular problems in shell theory: computing and asymptotics. Berlin, Germany: Springer Science & Business Media, Vol 54, 2010.

66.

Baue

Längst

Magna

, et al. Exploring software approaches for the design and simulation of bending active systems. In: Proceedings of the international association for shell and spatial structures (IASS) symposium 2018: creativity in structural design (eds Mueller

Adriaenssens

), Boston, USA, 16–20 July, pp.1–8. Madrid, Spain: International Association for Shell and Spatial Structures.