Optimal design of a long and slender compressive strut

This article deals with the optimal design of long and slender compressive struts. The main objective is to minimize the mass of the struts under certain non-failure constraints and thus find the optimal material. We show that the main failure mode of the struts is Euler buckling. The results clearly show that the struts should be constructed from unidirectional carbon fiber composites. A Monte-Carlo model for random microstructure homogenization of unidirectional composites is developed. We finish by performing a numerical computation of the effective properties of the chosen carbon fiber/epoxy composite using COMSOL MULTIPHYSICS software.


INTRODUCTION
Recent developments in material technology have introduced new computational challenges, such as the determination of thermo-elastic properties in a composite material or the calculation of flow through a porous medium.Several of these problems have in common that the governing equations involve rapidly oscillating functions due to the heterogeneity of the microstructure.These rapid oscillations render a direct numerical treatment impossible.Instead, one is forced to do some type of averaging or asymptotic analysis, which is the starting point for the concept of homogenization.
In this article, we start by analyzing a simple optimization problem, namely to minimize the mass of a long and slender compressive strut.We base the minimization on four basic load situations and determine material indices for each case.It quickly turns out that the best choice of material will be a unidirectional carbon fiber composite.However, to check the validity of the investigated load cases, we also find the dominating failure mode in section three.In section four, we give a brief introduction to the theory of homogenization and show that the effective material depends on the solution of a certain cell problem.Finally, in section five we develop a method to approximate a random microstructure by solving an ensemble of pseudorandom cell problems.We finish by giving a numerical example where we find the effective elastic tensor for a specific carbon fiber composite.

MATERIAL INDICES
Consider a long and slender compressive strut of length l and constant cross section A. The strut is pin-jointed in both ends, see Figure 1 below.We analyze the strut using the method of material indices described in [1].
We want to find the optimal material of the strut in the sense of minimum weight, under the four constraints: 3. The strut will not buckle.4. The strut will not fail due to impact damage.
In short, we have the following design requirements:

Objective:
Minimize the mass of the strut.
Will not buckle, F < F crit .Toughness against impact, σ < σ crit .Assume that the strut is made of a material with density ρ.Then we find the cross section A by (1) The shape of the cross section will greatly affect the behavior in for instance buckling due to different second moments of area.We therefore introduce a parameter ϕ that we call shape factor in order to determine the effectiveness of the chosen shape.For a given cross-section area A, take ϕ to be the ratio between the second moment of area I of the chosen shape to that of a solid circular section, I 0 , with the same cross section A. Now (2) where A 0 = A is the cross section of a solid circular section of radius r.Hence (3)

MATERIAL INDEX FOR FAILURE STRENGTH
If F is the applied force, the compressive stress σ in the strut will be (4) The constraint σ < σ f thus yields By rearranging, we find that (6) This means that in order to find the minimal mass of the strut under the failure constraint, we need to minimize ρ /σ f .We should thus maximize the material index (7) Taking logarithms on both sides, we see that we should choose the materials lying above and to the left of the line (8) in a ρ − σ f diagram using the software Cambridge Engineering Selector.See Figure 2 below.

MATERIAL INDEX FOR STIFFNESS
Hooke's law in one dimension states that (9) where E is Young's modulus for the material and the corresponding strain.Hence we get the elastic compression (11) The constraint δ < δ max thus yields (12) which implies that (13) This means that in order to find the minimal mass of the strut under the stiffness constraint, we need to maximize the material index (14) Materials with high values of M 2 are plotted in a ρ -E diagram using Cambridge Engineering Selector in Figure 3 below.

MATERIAL INDEX FOR EULER BUCKLING
The strut will buckle if the applied compressive load F exceeds the critical buckling load F crit given by (15) where the constant k depends on the end supports.For the case of pin-jointed ends, k is 1.
Using the shape factor ϕ to express I, the constraint F < F crit thus implies (16) Hence we get that (17) By rearranging and taking the square root, we find that (18) This means that in order to find the minimal mass of the strut under the Euler buckling constraint, we need to maximize the material index (19) Materials with high values of M 3 are plotted in a ρ -E diagram using Cambridge Engineering Selector in Figure 4 below.Note that the slope of the selection line is different compared to Figure 3, due to the changed exponent on E in M 3 .
Of course we can further minimize the mass by maximizing the shape factor ϕ.However, this will also affect the failure modes of the strut.This is analysis is carried out in the section Failure modes below.

MATERIAL INDEX FOR TOUGHNESS AGAINST IMPACT
We assume that the surface of the strut can be modeled as a thick plate where all energy from impact will be absorbed by contact indentation.This is in general not true since energy also will be absorbed in shearing, bending and membrane effects.However, the general situation is way too complicated to deal with here so we need to simplify the problem.On the other hand, we are not interested in a quantitative estimate of the absorbed energy, but rather the qualitative behavior in order to find the correct material index for this physical situation.We assume that the impacting object is spherical with radius R, mass m and velocity v.This gives an impact energy (equal to the kinetic energy) of (20) We also assume that all impact energy is absorbed by the indentation.A reasonable assumption for the contact stress is given by the Hertzian law (see e.g [2], [3], [4] or [5]) which states that (21) where F is the impact force, b the radius of the indented section, y the relative motion of the center of the ball compared to the indented surface and E the corrected modulus (22) where indices 1 and 2 corresponds to the ball and plate, respectively.We thus get that (23) and ( 24) . Optimal design of a long and slender compressive strut The Hertzian theory then predicts a maximal (tensile) contact stress in the plate of (25) when y = y max .Now the absorbed kinetic energy must be equal to the work done by indenting the plate, that is (26) Inverting this relation gives us the maximal indentation y max as (27) Hence the maximal stress becomes (28) The Griffith criterion for crack growth states that a crack of size a will grow if the stress exceeds σ crit given by (29) where Y is a geometric parameter close to one and K IC the fracture toughness of the material, see [6].To avoid fracture, we thus need to fulfill the inequality σ crit < σ max which means that (30) that is, (31) This means that in order to find the material that will absorb maximal impact energy under the no crack growth constraint, we need to maximize the material index π π

OPTIMAL MATERIAL CHOICE FROM MATERIAL INDICES
The graphs of the three material indices M 1 , M 2 , and M 3 suggest that the strut should be constructed from (unidirectional) carbon fiber composites.Possible other choices include aramid fiber composites, certain woods and certain extra durable steels.Moreover, these graphs also show that certain ceramics have high values of the material indices, but we exclude them due to their brittle nature as well as their somewhat complicated manufacturing process.We also exclude diamond for obvious reasons.
Material index M 4 suggests that carbon fiber composites with fiber direction 90°are especially good to absorb impact energy, due to their ability to capture cracks growing transversally to the fiber direction.Therefore, in order to get a good composite that will withstand fracture, we need to laminate the composite with plies having fiber directions other than 0°.We will see in next chapter that material index M 3 is of special importance, since the dominating failure mode will be Euler buckling.

FAILURE MODES OF A TUBULAR COLUMN
For simplicity, we consider a thin walled tubular column with circular cross section loaded in compression by a force of F. Suppose that the column is made from a material having stiffness E, Poisson's number ν and yield stress σ f .Moreover, it has length l, outer radius R o and inner radius R i .We also assume that both ends are pin-jointed.Let t be the wall thickness of the tube and r its mean radius, i.e.We then have that (34) which gives us the shape factor for a thin walled tube (t<<r) as (35) We consider four different failure mechanisms: 1. General (Euler) buckling, where the critical buckling load is Using the shape factor ϕ, we get that Hence the failure stress for mechanism 1 is where α is a parameter in the range 0.4-0.6,see [7].see [8].

Boundary 3-4:
(47) We write the first three boundaries in logarithmic form The last three boundaries only depend on ϕ, but not on the load factor F/(σ f l 2 ).These three boundaries coincide when the shape factor ϕ equals the value ϕ 34 , that is, which occurs for a failure strain σ f /E of about (50 when α is taken to be 0.5.Since this value is much higher than the failure strain of all interesting materials, we conclude that we always will have the same order of the boundaries 2-3, 2-4 and 3-4, independent of the failure strains we will consider.
The dominating failure mechanism will be the one that occurs for the lowest stress σ.Thus we make the following comparisons: 1. Comparison of σ 1 and σ 2 .
(56) which implies that failure mechanism 3 dominates mechanism 4 left of the 3-4 line and that mechanism 4 dominates mechanism 3 right of the 3-4 line.
We now assume that the column is made by a carbon fiber composite with a typical failure strain of ε f = σ f /E = 1.4%.We construct a diagram with logarithmic scales where the four failure mechanisms are shown.Let the y-axis be the load factor F/(σ f l 2 ) and the x-axis be the shape factor ϕ. First we note that the failure strain ε f = 0.014 yields that (57) which means that ϕ 34 < ϕ 23 < ϕ 24 .Considering all details above, we can construct the map shown in Figure 6 above.For a typical long and slender strut with a load factor of about 10 −6 and shape factor in the range 5-20, we see that the dominating failure mechanism will be general Euler buckling.

Remark:
The actual profile of the strut is of course more always complicated than a circular tube.However, we can do a similar analysis using the corresponding shape factor ϕ for the direction in which the strut is mostly prone to buckle (that is, in the direction in which the second area moment of inertia is smallest).Then the failure of the strut will be described by a similar diagram as above and the general conclusion will be the same.

HOMOGENIZATION OF HETEROGENOUS MEDIA
We have come to the conclusion that we need to construct the strut from a carbon fiber composite.However, composites introduce a heavy computational problem when analyzed by typical finite element software due to their inhomogeneous microstructure.More specific, the elastic constants for fibers and matrix will oscillate fast in space, making the problem virtually impossible to solve numerically.Instead one looks for a so called effective material that behaves like the composite from a macroscopic point of view.The mathematical problem of finding the effective material is called homogenization.
Let us consider a linear elastic body which occupies a region Ω in R 3 and introduce a Cartesian coordinate system x = (x i ).Moreover, we introduce σ = (σ ij ), f = (f i ), t = (t i ), u = (u i ) and n = (n i ) as the stress tensor, the internal force field, the surface force field, the displacement field and the outer normal to the boundary ∂Ω of Ω, respectively.The governing equilibrium equation is then where Γ 1 ∪Γ 2 =∂Ω and Γ 1 ∩Γ 2 =∅.We define the strain tensor e=(e ij ) by (59) Hooke's law states that the stress σ is related to the strain e by the relation σ ij = C ijkl e kl , where C ijkl is the fourth order elastic tensor (or stiffness tensor) of the material that occupies Ω.It is well known that the stiffness tensor has the symmetries C ijkl = C jikl = C ijlk = C klij , effectively reducing the number of independent constants in the stiffness tensor.Further symmetries in the stress and strain tensors, σ ij = σ ij and e ij = e ji , respectively, then imply that we can represent Hooke's law in the (symmetric) matrix form where the engineering strains γ ij are defined as (61) The equilibrium equation thus becomes (62) Let us now assume that the body consists of two or more different linear elastic materials which are periodically distributed in the sense that we can define unit cells Y which are periodically repeated, see Figure 7 below.
For such a material, the stiffness tensor C ijkl =C ijkl (x) will oscillate periodically.To model the heterogeneous material, we therefore introduce a local variable y=x/ε and assume that (63) is Y-periodic.By Y-periodicity we mean that C ijkl (y 1 ) = C ijkl (y 2 ), whenever y 1 and y 2 occupy the same positions in their corresponding cells.This means that ε is a parameter for varying the fineness of the microstructure, see Figure 8.
We also assume that the functions C ε ijkl satisfy the following coercivity and continuity conditions (64) for every symmetric real-valued tensors ξ ij where 0 < α ≤ β < ∞.Physically this means that the strain energy is positive and bounded.
We now study the following class of problems, one for each choice of ε, The main idea in the homogenization theory is to approximate the solutions u ε of our model problem by means of a function u which solves the problem corresponding to a homogeneous material, (66) where C * ijkl is a constant tensor and g i is defined below.Mathematically speaking, the solutions u ε converge weakly to the solution u as ε → ∞.The homogenized tensor may be interpreted as the physical parameters of a homogeneous material whose overall response is "close'' to that of the heterogeneous material, when the size of the cell tends to zero.We say that C * ijkl is the effective stiffness tensor.Using some rather deep results from functional analysis and the theory of partial differential equations, it is possible to show that the effective tensor C * ijkl can be given explicitly by the formula where U mn =U mn (y) are the weak solutions of the Y-cell problem (68) Examples of Y-cells for a periodic material.
Figure 8 The microstructure becomes finer as the parameter ε becomes smaller (left to right).
See for instance the books by [9], [10], [11], [12], [13], [14] or [15] for proofs of this result.From the proof of the theorem, we also get the result that (69) Let us now consider this cell problem.We have that which can be rewritten as where V i mn =δ im y n which, written out in full, says that (72) etc.This gives the strains (73) Hence we have that (74) etc., where Thus if we put W mn =U mn +V mn , our cell problem becomes 1 0 0 0 0 0 0 0 0 0 0 1 ( , , , , , ) , e 22  0 1 0 0 0 0 y y with the boundary conditions that points on opposite faces (with normals parallel to the y m -direction) of the cell are coupled to each other such that they move identically, except in the y m -direction where they will move equally up to the difference ∆y m (the length of the cell in the y m -direction).We now solve the cell problem for all W mn using the corresponding boundary conditions.We then find the effective stiffness coefficients as (77 where σ mn =(σ ij ) mn are the corresponding stress tensors from the cell problem for W mn .In other words, we get the effective stiffness coefficients C * ijmn as the average stresses σ ij from the solutions of the cell problems with the boundary conditions described above.See for instance the articles [16], [17], [18] and [19], for more applications of the homogenization method to composite materials.

NUMERICAL COMPUTATION OF EFFECTIVE STIFFNESS
We showed earlier that we need to solve a family of cell problems in order to find the effective properties.We will here give a detailed description of a model to generate a composite with pseudorandom microstructure and then give a numerical example on homogenization of it.

BACKGROUND
The cell problem is defined over a periodic cell with the assumption that the underlying microstructure is periodic.However, the choice of microstructure will greatly influence the numerical values we get.Should we for instance use rectangular or hexagonal stacking of the fibers?
The rectangular cell (see Figure 9) has the advantage that it is easy to implement and the hexagonal cell (see Figure 10) has the advantage that it will give a transversely isotropic material.However, most people will probably agree that the proper choice of unit cell is neither of them -fibers are usually randomly distributed in a composite.
So how can we model such a random structure with our homogenization procedure?A popular method is to use some kind of Monte-Carlo algorithm to generate a "large enough'' random microstructure and impose periodic boundary conditions, see e.g.[20], [21], [22] and [23].See Figure 11 for a pseudorandom unit cell with 1000 fibers.
What does then large enough mean?One would argue that the more fibers, the better agreement with reality.On the other hand, the accuracy of the solution declines as the number of fibers increase.A compromise of these both extremes can be to generate a large family of random structures, each with a reasonable number of fibers, and take the average of all generated coefficients, see e.g.[24], [25], [26] or [27].
We chose to follow this latter principle for our computations.We used a Metropolis type algorithm (see [28]), starting from a periodic rectangular array of 64 fibers.Each fiber is then given a small tentative displacement in a random direction.The move is accepted or rejected whether or not the move will cause the fiber to overlap a neighboring fiber.One iteration step consists of trying to move each fiber once.The length of the displacement is chosen such that the ratio of acceptance is 30-50%.The procedure is repeated several times until an Each cell problem was then solved using COMSOL MULTIPHYSICS and the corresponding effective properties for each realization were computed, see Figure 13 below for an example with 9 fibers.Finally, the mean values were taken for each coefficient.The choice of 64 fibers in each cell is motivated by the results in the articles [25] and [27], where it is shown that the standard deviation (as expected in a Monte Carlo simulation, see e.g.[29]), is reciprocally proportional to the square root of the number of fibers.However, in this case when we deal with only a few (2-3) significant figures, the mean values for the effective properties are stable already at one fourth of this number of fibers.So there is no need to lose accuracy by choosing too many fibers in each cell.

A NUMERICAL EXAMPLE
Let us consider a typical carbon fiber composite consisting of transversely isotropic carbon fibers with a volume fraction of Vf=60% having the mechanical properties E L =230 GPa, E T =40 GPa, v TT =0.20, v LT =0.256, G TT =16.7 GPa, G LT =24 GPa, and an isotropic polymer matrix having the properties E=2.9 GPa, v=0.30.
It is a well known fact (see e.g.[30], [31], [32] or [33]) that orthotropic materials described by the nine constants E 1 , E 2 , E 3 , v 12 , v 13 , v 23 , G 12 , G 13 , G 23 have the stiffness tensor (in matrix form)      (79) Furthermore, isotropic materials have infinitely many planes of isotropy, yielding Thus we easily find the corresponding stiffness matrices for the matrix and the fibers.Solution of the 1000 cell problems in COMSOL MULTIPHYSICS then yielded an average effective stiffness tensor C=C ijkl with the elements Moreover, since (83) we see that the resulting homogenized composite (not surprisingly) is transversely isotropic, having the x 2 x 3 -plane as the plane of isotropy and x 1 as the longitudinal direction L. Inverting the stiffness matrix gives us the corresponding compliance matrix, from which we finally extract the effective engineering constants

Figure 2
Figure 2 Materials with high index M 1 lie above the indicated line.

Figure 3
Figure 3 Materials with high index M 2 lie above the indicated line.

Figure 4
Figure 4 Materials with high index M 3 lie above the indicated line.
Int. Jnl. of Multiphysics Volume 3 • Number 3 • 2009 Materials with high values of M 4 are plotted in an E − K IC diagram using Cambridge Engineering Selector in Figure 5 below.

Figure 5
Figure 5 Materials with high index M 4 lie above the indicated line.

3 .
General yield (crushing) of the column occurs at the failure stress (40) 4. Static axial plastic buckling occurs at the failure stress (41) which implies that failure mechanism 1 dominates mechanism 3 below the 1-3 line and that mechanism 3 dominates mechanism 1 above the 1-3 line.3. Comparison of σ 1 and σ 4 .
which implies that failure mechanism 1 dominates mechanism 4 below the 1-4 line and that mechanism 4 dominates mechanism 1 above the 1-4 line.4. Comparison of σ 2 and σ 3 .
which implies that failure mechanism 3 dominates mechanism 2 left of the 2-3 line and that mechanism 2 dominates mechanism 3 right of the 2-3 line. 5. Comparison of σ 2 and σ 4 .(55)which implies that failure mechanism 4 dominates mechanism 2 left of the 2-4 line and that mechanism 2 dominates mechanism 4 right of the 2-4 line.

4 F f l 2 σFigure 6 A
Figure 6 A map of dominating failure mechanisms.
obtained, see Figure12.We used 2000 iterations for each of totally 1000 different realizations of pseudorandom cells.

Figure 9
Figure 9 Stress distribution (in x 1 -direction) in a rectangular unit cell.

Figure 10
Figure 10 Stress distribution (in x 1 -direction) in a hexagonal unit cell.

Figure 11 A
Figure 11 A pseudorandom unit cell with 1000 fibers.

Figure 12
Figure 12 One realization of a pseudorandom unit cell with 64 fibers.
isotropic materials are special cases of orthotropic materials having one plane of isotropy (here the x 2 x 3 -plane) where                                                                                                                                           .

Figure 13
Figure 13 Stress distribution in the x 1 -direction for a pseudorandom cell with 9 fibers. C