Optimal Design of Fibre Reinforced Membrane Structures

A design problem of finding an optimally stiff membrane structure by selecting one-dimensional fiber reinforcements is formulated and solved. The membrane model is derived in a novel manner from a particular three-dimensional linear elastic orthotropic model by appropriate assumptions. The design problem is given in the form of two minimization statements, reminiscent of a Nash game. After finite element discretization, the separate treatment of each of the two minimization statements follows from classical results and methods of structural optimization: the stiffest orientation of reinforcing fibers coincides with principal stresses and the separate selection of density of fibers is a convex problem that can be solved by optimality criteria iterations. Numerical solutions are shown for two particular configurations. The first for a statically determined structure and the second for a statically undetermined one. The latter shows related but non-unique solutions.


Introduction
A finite element membrane shell model was recently derived by Hansbo and Larson [8] using tangential differential calculus, meaning that the problem is set in a Cartesian three dimensional space as opposed to a parametric plane, thereby generalizing the classical flat facet element shell model to higher order elements. The present study further extends this membrane model by allowing for nonisotropic materials. In particular, one-dimensional fibers are added to a base material, modeling, e.g., the reinforcements seen in modern racing boat sails. The plane stress property, as well as the membrane property of complete out-of-plane shear flexibility, is shown to be exact consequences of certain material parameter selections for a three-dimensional transversely isotropic base material. This together with a displacement assumption reduces the three dimensional model to the surface model. Based on this finite element model we formulate a design problem where we seek to find the best fiber reinforcements of the membrane, meaning that we find the stiffest structure by both rotation and sizing of the fibers. The formulation is reminiscent of a Nash game [1,10], consisting of two minimization statements. However, the two players of the game have the same objective, i.e., stiffness, as opposed to standard game theory. Nevertheless, since the two minimization statements relates to rotation and sizing of the fibers, respectively, such a formulation ties directly to the sequential iterative treatment suggested for similar problems previously [3]. The optimal rotation is found by identifying the material as a so-called low shear material, implying that the optimal orthotropic principal directions coincides with the principal stress directions [13,14], while the optimal thickness distribution is found by a classical optimality criteria iteration formula.

The model
We consider a material that is a mixture of a transversely isotropic linear elastic base material and n reinforcing fibre materials. The transversely isotropic material has material constants that satisfy the plane stress assumption as well as the membrane behaviour of having complete flexibility when sheared perpendicularly to the membrane surface.

Geometry
The geometry of the membrane is defined by an orientable smooth surface Σ with outward normal n. For any point x ∈ R 3 we denote the signed distance function relative to Σ by ζ(x). The membrane with thickness t then occupies Note that ∇ζ(x) = n for x ∈ Σ. For a sufficiently small t, the orthogonal projection point p(x) ∈ Σ of x ∈ Ω t is unique and given by Moreover, for x ∈ Ω t , the linear projection operator of vectors onto the tangent plane of Σ at p(x) is where I is the identity tensor and ⊗ denotes exterior product. In the sequel we will also need the projection operator onto the one-dimensional subspace spanned by n, i.e., Note that P Σ N Σ = N Σ P Σ = 0. The directions of the reinforcing fibers are given by vector fields s i , i = 1, . . . , n, such that s i · n = 0. Projections onto these directions are then defined by

The material
The base material is transversely isotropic with respect to an axis defined by n. Such a material can be described by an elasticity tensor expressed in terms of five material constants δ 1 , δ 1 , δ 1 , γ and µ according to [11,12], so that the fourth order tensor of elastic moduli of the base material E base can be written Here dyadic products of second order tensors are defined by their action on a third second order tensor, i.e., where a double dot indicates inner product of second order tensors. The reinforcing fibers have elasticity tensors of the form where α i are Young type elasticity coefficients. The constitutive law of the membrane material is now taken as being composed of a constrained mixture of base material and reinforcing material. The amount of each material is defined by fractions t b and t i , i = 1, . . . , n, of the membrane thickness t such that the total constitutive tensor is given as and the linear constitutive law is where σ and ε are the stress and strain tensors, respectively.

Membrane stress assumptions
We define a membrane material by the requirements that it is always in a state of plane stress and no shear stress perpendicular to the membrane surface exists, i.e., The zero bending stiffness behaviour of membranes will be a result of a kinematic assumption introduce subsequently. Inserting (3) into (4) gives Thus, we conclude that the constitutive constant γ needs to be zero and that the strain perpendicular to the membrane is controlled by the in-plane strain as Moreover, the in-plane stress can be calculated from (3) as follows: and when using (7) we get where The elasticity coefficient µ equals the in-plane shear modulus, while δ is a plane stress Lamé coefficient. The two elasticity moduli δ and µ can be expressed in terms of in-plane Young and Poisson moduli E and ν as The volumetric specific strain energy can, due to (4) be written as Inserting (8) we get where the membrane elasticity tensor is defined by

Potential energy
The strain is derived as usual as the symmetrized gradient of the displacement vector u: Therefore, we can regard the volume specific strain energy as a function of the displacement field, i.e., W s = W s (u). We now introduce the basic kinematic assumption that all material points in Ω t that lie along a normal to the surface Σ have the same displacement vector, i.e., This kinematics imply that bending of the membrane is essentially eliminated and no bending stiffness, despite the finite thickness, is present.
The total strain energy, which is the volume integral of W s can then be written: where dΣ ζ is an area element for a surface parallel to Σ at the distance ζ, which reads where dΣ is the area element of Σ, and H and K are the mean curvature and Gaussian curvature, respectively. For a membrane that is thin compared to its curvature we can use the approximation The total potential energy is now taken as where the force F is a member of the dual space of displacement fields on Σ and ·, · Σ is a duality paring.
We define the membrane forces (per unit length) as Stationarity of the potential energy gives the following principle of virtual work: for all kinematically admissible fields v. Such fields will generally be restricted in the tangential direction on a subset of ∂Σ. We will assume that loading on the membrane can be written as where f is a force per area over Σ, and p is a force per unit length over the part S of ∂Σ where the displacement is not prescribed. Using now Lemma 2.1 of Gurtin and Murdoch [6], i.e., an integral theorem for surfaces, we obtain the equilibrium equations where div Σ is the surface divergence, and ν is a unit vector of ∂Σ, tangential to Σ. Since M ν will also be a vector tangent to Σ we conclude that p can have no component perpendicular to the surface.

Design problem
From now on we will consider the special case of an orthotropic material consisting of two orthogonal families of fibers, consisting of the same material, i.e., α 1 = α 2 = α. From now on we use the notation s = s 1 and s ⊥ = s 2 .
The orientation of the fibers in the tangent plane of the membrane, i.e., s and s ⊥ , can be defined by an angle θ belonging to This angle will be a design variable in the optimal design problem. Other such design variables are t 1 and t 2 , i.e., the fiber contents in the two orthogonal directions. The field t = (t 1 , t 2 ) belongs to the set where t α and t α are non-negative upper and lower bounds and V is a limit for the total amount of material that can be used for the fibers.
The potential energy is seen as a function where V is the set of kinematically admissible displacements. Minimizing Π with respect to the first argument gives the equilibrium displacement as a function of the design variables, i.e., u = u(t, θ). As a measure of stiffness we use the so called compliance Our design goal is to find a design that minimizes the compliance. We choose to split this into two parts as follows: find t * ∈ T and θ * ∈ Θ such that This is reminiscent of a Nash equilibrium problem [1,10], but the two players have the same objective in contrast to natural Nash games where objectives are opposing, which in the present case would happen if minimization is changed for maximization in one of the two lines of (P). The second sub-problem of (P), i.e., finding an optimal orientation for an orthotropic material, has been extensively discussed by Pedersen [13,14] and Hammer [7], see also Bendsoe and Sigmund [3] for further discussion and related references. It turns out that the problem can be solved locally, i.e., the orientation of the material is determined by the local stress state only, and in particular the orientation of principal stresses and strains. Due to the plane stress assumption there are only two possibly non-zero principal components of the stress tensor σ, denoted σ I and σ II , such that |σ I | ≥ |σ II |. The corresponding principal directions (eigenvectors) are tangent to the membrane plane. Obviously, these facts also holds for the principal components of M , i.e., M I and M II , such that |M I | ≥ |M II |. For a so-called low shear orthotropic material, the solution θ * of the second sub-problem of (P) represents an orientation where the orthotropic principal directions coincide with the principal stress or membrane force directions, which are also the principal strain directions. Moreover, the orthotropic principal direction having that highest stiffness should be in the direction corresponding to σ I and M I . In the Appendix we show that the particular orthotropic material defined above, having by two families of fibers in orthogonal directions s and s ⊥ , is indeed a low shear material and, therefore, the optimal directions of s and s ⊥ are in the directions of principal stress. Moreover, if t 1 > t 2 then s is in the direction of σ I .
The first sub-problem of (P) is a classical stiffness optimization problem, albeit having two design fields, one for each fiber orientation. This is a convex problem and can be solved by satisfying the optimality conditions. The surface elasticity tensor S memb = t E memb is regarded as a function of the design, i.e., S memb = S memb (t, θ). The optimality conditions of the first sub-problem of (P) become [3,4]: where Λ, λ + α and λ − α are Lagrangian multipliers, t ∈ T and u = u(t, θ) is the displacement solution, i.e., the minimum field with respect to v of Π(v, t, θ).

Discretization and algorithm
For the numerical treatment of (P) we need to introduce a discrete approximation. The discretization of the state problem, i.e., the problem of finding the minimum displacement u ∈ V of the potential energy Π for a given design θ ∈ Θ and t ∈ T , follows Hansbo and Larson [8]. This implies introducing a triangulation of Σ resulting in a discrete surface, with corresponding discrete normal vector field and projections. The displacement field is approximated using the same triangulation but is possibly of different order.
In addition to the approximation of the state problem we also need to approximate the design fields t ∈ T and θ ∈ Θ. This is achieved by using point values: these are denoted t i = (t 1i , t 2i ) and θ i for point i. In particularly, we use superconvergence points of the finite elements [2]. Such a discretization means that (12) and (14) are imposed at these evaluation points and the integral in (13) is replaced by a sum. Let be the left hand side of (12) evaluated at point i and for a displacement field u k . Also, let B k αi = (Λ k ) −1 A k αi where Λ k is a currant iterate of the Lagrangian multiplier Λ. For a given displacement iterate u k and rotation θ k the following fixed point iteration formula is suggested by the optimality conditions (12) through (14): where t αi and t αi are point values of the upper and lower bounds and 0 < η ≤ 1 is a damping coefficient.
The following algorithmic steps, the convergence of which gives satisfaction of a discrete version of the optimality conditions of (P), is now suggested: 1. For a given design θ k and t k , solve the state problem, i.e., find the minimum displacement field of Π(v, t k , θ k ) so as to obtain the currant displacement iterate u k . 2. Obtain new fiber thickness distributions by the optimality criteria formula (15) where -Λ k is determined such that A local iteration is needed for this. 3. For each integration point, calculate principal stresses (and/or principal membrane forces). Take s to correspond to the main material direction, i.e., to t 1i , such that t 1i ≥ t 2i , and chose θ k+1 so that this s aligns with the main principal stress direction. 4. Let k = k + 1 and return to the first step.
Steps 1 and 2 can be iterated several times before continuing with calculation of fiber directions in Step 3. In fact, in the examples the fixed point iteration (15), for newly calculated displacement u k , is repeated until convergence before continuing with the fiber directions in Step 3. Note that step 3 assumes distinct principal stresses. Numerically coalescence of such stresses occur with close to zero probability but may show up as nonconvergence issues. For statically determined structures, i.e., when M is uniquely determined by (10) and (11), this may be of particular concern. For such cases that have distinct principal stresses, step 3 above needs to be performed only ones since these principal stress are independent of t. Such problems essentially becomes convex since the first part of (P) is a convex problem. The first problem of the Section 6 is statically determinate but has not everywhere distinct principal stresses.

Oblate spheroid
An oblate spheroid, where Σ is defined by was solved by different finite elements and triangulations in Hansbo and Larson [8]. Here we treat the same geometry but use an internal pressure p as loading. We seek for optimal fiber distribution as described in previous sections. The data are E = 1, ν = 0.3, t b = 0.005, p = 10, V = 0.01, t 1 = t 2 = 0.004, t 1 = t 2 = 0 and α = 1. The initial fiber thickness is uniform and chosen so as to satisfy the volume constraint as an equality. We use 3072 bilinear 4-node fully integrated isoparametric elements, implying one superconvergent point per element and, thus, three design variables per element. Symmetry is utilized and only half of the spheroid is modeled. The problem converged in 36 optimality criteria updates and 7 updates of the fiber orientations. As convergence criteria an objective value change below 0.001 % and a change of θ such that cos θ > 0.999 are used. Note that the problem is statically determinate, but at the poles of the spheroid symmetry implies that the principal stresses coincide for an exact solution. This is the reason for the need of several updates of fiber orientations before convergence, despite the problem being statically determinate.
What concerns the general features of the solution one finds, on examination of Figure 1, that close to the equator both fiber families are present, with a compressive stress in the latitudinal direction. As we move towards the poles only the longitudinal fiber family is present, while at the very poles the principal stresses coincide and the direction of fibers becomes indeterminate.

Membrane strip
A rectangular membrane of shape 1 × 0.5 is fixed along one of its short sides and loaded by a force q per unit length on a part of length 0.1 of the other short side, as shown in Figure 2. The date are E = 1, ν = 0, t b = 0.005, q = 0.001, V = 0.01, t 1 = t 2 = 0.008, t 1 = t 2 = 0 and α = 2. As in the previous example, the initial fiber thickness is uniform and chosen so as to satisfy the volume constraint as an equality.
The left hand solution of Figure 2 is found using initial fiber directions defined by the rectangle sides. The right had solution, on the other hand, uses initial directions defined by principal stress directions found in an initial calculation where fibers are excluded. The left hand problem converged, using the same tolerances as in the previous problem, in 28 optimality criteria iterations and 12 updates of the fiber directions. The right hand problem converged in 15 optimality criteria iterations and 6 updates of the fiber directions. The slightly difference between the two solutions is likely the result of a possible non-uniqueness of the solution of problem (P). How-ever, the objective function values for the two cases are essentially the same. Fig. 2 Optimal fiber distribution for a rectangular membrane using two different initial fiber directions.

Conclusions
The classical facet approach to membrane shells was recently extended to curved elements by Hansbo and Larson [8]. Here we make a further extension by showing how orthotropic material, of fiber type, can be treated in a similar way, partly inspired by exact plate theory of Nardinocchi and Podio-Guidugli [12]. Based on this orthotropic membrane shell theory we formulate a stiffness design problem, where we seek an optimal structure by both rotation and sizing of reinforcing fibers. The two design variables -representing rotation and sizing -naturally splits the formulation into two minimum statements, reminiscent of a Nash game, which suggests a sequential numerical treatment, previously used for similar problems [3]. This type of formulation also makes clear the distinct character of statically determined problems, which occur for large classes of membrane shells [5]. For such problems, the material independent stress state implies that the two minimization statements of (P) become decoupled, and since the sizing of fibers is a convex problem, the full problem (P) essentially inherits this property.
The approach presented in this paper has several intriguing extensions, that would be important for applications such as the design of racing boat sails. Inclusion of pre-stress and wrinkling states related to negative stresses are examples of this. Extension to large deformations, based on the model of Hansbo et al. [9], should also be of clear interest.