Buckling of Planar Micro-Structured Beams

Abstract

In this paper, a Timoshenko beam model is formulated for buckling analysis of periodic micro-structured beams, uniformly compressed. These are planar grid beams, whose micro-structure consists of a square lattice of equal fibers, modeled as Timoshenko micro-beams. The equivalent beam model is derived in the framework of a direct one-dimensional approach and its constitutive law, including the effect of prestress of the longitudinal fibers, is deduced through a homogenization approach. Accordingly, micro–macro constitutive relations are obtained through an energy equivalence between a cell of the periodic model and a segment of the equivalent beam. The model also accounts for warping of the micro-structure, via the introduction of elastic and geometric corrective factors of the constitutive coefficients. A survey of the buckling behavior of sample grid beams is presented to validate the effectiveness and limits of the equivalent model. To this purpose, results supplied by the exact analyses of the equivalent beam are compared with those given by finite element models of bi-dimensional frames.

1. Introduction

Research on periodic micro-structured systems and architected materials have recently received great attention in several fields of applied sciences [1,2], also due to significant advances in 3D printing techniques [3]. Indeed, micro-structured bodies reveal some exotic characteristic which far exceed those of classical bulk materials, compared to weight, like, e.g., higher mechanical performances, resilience, impact resistance and, in general, energy absorption capability [4,5,6,7,8]. The lattice of fibers, i.e., the micro-structure, with its mechanical properties and topology, confers to the resulting (macro) structure such peculiar characteristics.

However, periodic lattices are often characterized by high slenderness and can be particularly sensitive to local and global buckling phenomena. For example, in [9,10,11,12,13,14,15], buckling- and plasticity-driven failure mechanisms are investigated in honeycombs, under different loading conditions. In [16], buckling and post-buckling of elastic square honeycombs, under in-plane biaxial compression, is analyzed using a two-scale approach; a simple formula for the evaluation of the very-long-wave buckling stress is also provided. In [17], buckling of a square lattice is studied also to address the possibility of tuning the buckling shape by varying the cell geometry. In this framework, a very promising aspect in dealing with buckling of cellular structures, is the idea to transform the usually (historically) undesirable effects of buckling phenomenon from a negative into a positive, for their potential contribution to ‘smart’ applications in several fields [18] as, e.g., in biology, mechanics, physics, engineering, aerospace and so on.

One of the most powerful tools in modeling micro-structured lattices as equivalent continua is homogenization [19,20,21,22,23,24,25,26,27,28,29,30,31] and several methods, see, e.g., [20,26], have been developed. In this field, generalized (continuous) beam models have been proposed in [32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54] to analyze statics, dynamics, buckling and aeroelasticity of different micro-structured bodies, as lattices, and civil structures as tower buildings and multi-storey frames.

In a set of the mentioned papers, [43,44,45,46,47,49,50,51,52,53], buildings and planar frames have been modeled, through a homogenization approach, as continuous beams, namely as: shear-shear-torsional beams, by neglecting macro-bending, or Timoshenko beams. In this framework, the floors and the columns of the periodic frames have been identified with the cross-sections and longitudinal fibers of the continuous beam, respectively. The reader is referred to [53], where a review of these papers has been provided together with a systematization of the homogenization procedure. In [43,44,45,46,47,49,50], the floors have been assumed to behave as rigid bodies. More details regarding the shear-shear-torsional beam model: in [43,44], this model has been proposed to analyze the aeroelastic behavior of tower buildings; in [47], the main assumptions and the limits of applicability have been discussed; in [50], the beam model, under small displacement assumption, has been proved to be effective in analyzing linear free and forced dynamics of multi-storey buildings. On the other end, with reference to the Timoshenko equivalent beam model: in [45], such a model has been proposed to properly account for macro-bending of multi-storey buildings in linear statics and dynamics; in [46], the model, embedded in a 3D space, is improved to account for buckling analysis of uniformly compressed multi-storey buildings; in [49], an extension of [46] is proposed, where the same equivalent beam model is shown to be effective in describing the flexural-torsional buckling behavior of not-uniformly compressed tower-buildings, also accounting for the effect of an elastic soil. Moreover, in [51,52,53], it is shown that it is possible to account for the flexibility of the cross-section, in a direct 1D (unwarpable) beam model, as, indeed, the shear and the Timoshenko ones. To this end: (i) corrective factors have been analytically derived in [51], which correct the constitutive coefficients of a planar equivalent model, thus accounting, on (energy) average, for the neglected warping; (ii) an energy-based numerical algorithm is presented in [52], where the constitutive properties of the planar equivalent beam model have been determined through a finite-element (FE) analysis of a single cell.

In a companion paper [54] the Timoshenko beam model so far developed in [51,52] is extended to the 3D space in order to describe statics and dynamics of spatial grid beams, namely grid cylinders with cubic micro-structure (see, e.g., [5]), consisting of a periodic cubic pattern of orthogonal micro-beams.

In this paper, the buckling behavior of uniformly compressed planar grid beams is analyzed through an equivalent Timoshenko beam model (coarse model), embedded in 2D space and derived through a homogenization procedure. The formulation follows the lines of [46,49,53] and accounts for the effect of prestress of the longitudinal fibers, both in the equilibrium and in the constitutive law. In addition and differently from these papers, the proposed equivalent Timoshenko beam model takes into account the warping of the micro-structure, via the introduction of elastic and geometric corrective factors of the constitutive coefficients, by extending what done in [51]. Buckling analyses of sample grid beams are presented aiming at the validation of the equivalent model against benchmark finite-element (FE) solutions of bi-dimensional frames (fine models).

The paper is organized as follows: In Section 2, the equations of the Timoshenko beam model are recalled, while in Section 3, the identification procedure for the constitutive constants is detailed. Section 4 shows the results of numerical analyses. Finally, in Section 5, some conclusions and perspectives of the present work are drawn.

2. Model

A planar grid beam is considered, consisting of a planar frame, made of two order of micro-beams (namely the fibers), which are arranged in a periodic square pattern of side h, as displayed in Figure 1a. The families of fibers are: (i) that of the longitudinal x-fibers, which are m in y-direction; (ii) that of the transverse y-fibers, which are $n+1$ in x-direction. The nodal points of the y-fibers constitute the cross-section ${\mathcal{D}}_i$. Throughout the paper it is assumed that the x-axis is a symmetry axis for the grid beam and the origin of the coordinate system is placed at the left end A (see Figure 1a). All the fibers are modeled as Timoshenko micro-beams, having axial, shear and bending stiffnesses $EA$, $G{A}^{*}$ and $EI$, respectively. The x-fibers of the grid beam are uniformly compressed by axial forces, whose resultant magnitude is P.

The coarse model of the grid beam is (heuristically) assumed to be the Timoshenko beam model, embedded in a 2D space (Figure 1b), whose constitutive law needs to be suitably defined, as described in what follows. Accordingly to this choice, the internally unconstrained model is referred to the material abscissa $x\in \left[0,\ell \right]$, belonging to the interval $\left(A,B\right)$, of length $\ell :=nh$, and the following strain–displacement relationships hold [55,56]:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \epsilon & =u^{\prime },\hfill \\ \hfill \gamma & =v^{\prime }-\theta ,\hfill \\ \hfill \kappa & ={\theta }^{\prime },\hfill \end{array}\end{array}$$

(1)

where a prime denotes space-differentiation and, moreover, $u\left(x\right)$, $v\left(x\right)$ and $\theta \left(x\right)$ are the longitudinal and transverse displacements, and the rotations of the cross-section; these displacements are linked by Equation (1) to the elongation $\epsilon \left(x\right)$, the shear strain $\gamma \left(x\right)$ and the flexural curvature $\kappa \left(x\right)$. Different restraints at the boundary $H=A,B$, namely at $x=0$ and $x=\ell$, are considered, for which the relevant boundary conditions read as follows:

• clamp (C) at H: $u_H=v_H={\theta }_H=0$;
• hinge (H) at H: $u_H=v_H=0$;
• horizontal roller (R) at H: $v_H=0$;
• horizontal slider (S) at H: $v_H={\theta }_H=0$.

To analyze the buckling problem, equilibrium of the Timoshenko beam is enforced in the adjacent (slightly varied) configuration. This is done via the virtual work principle (VWP), which, as usual, also accounts for the virtual work spent by the prestress $\mathring{N}=-P$ on the second-order part of the unit-extension ${\epsilon }^{\left(2\right)}$. Once the exact unit-extension $\epsilon$ is expanded in series, it reads [55,57,58,59,60,61]:

$$\begin{array}{cc}\hfill \epsilon & :={\epsilon }^{\left(1\right)}+{\epsilon }^{\left(2\right)}+\dots =u^{\prime }-\frac{1}{2}{\theta }^2+v^{\prime }\theta +\dots .\hfill \end{array}$$

(2)

In the absence of incremental loads, the VWP reads:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \underset{0}{\overset{\ell }{\int }}\left(N\delta u^{\prime }+T\delta \left(v^{\prime }-\theta \right)+M\delta {\theta }^{\prime }\right)dx-\underset{0}{\overset{\ell }{\int }}P\delta \left(v^{\prime }\theta -\frac{1}{2}{\theta }^2\right)dx=0,\forall \left(\delta u,\delta v,\delta \theta \right),\end{array}\end{array}$$

(3)

which supplies the following equilibrium equations:

$$\begin{array}{cc}& N^{\prime }=0,\hfill \\ \hfill & T^{\prime }-P{\theta }^{\prime }=0,\hfill \\ \hfill & M^{\prime }+T+P\left(v^{\prime }-\theta \right)=0,\hfill \end{array}$$

(4)

$N\left(x\right)$, $T\left(x\right)$ and $M\left(x\right)$ being the (incremental) axial and shear internal forces, and bending moment, respectively. The mechanical boundary conditions, considered in what follows, are determined from the principle (3); they are:

• free (F) at H: $N_H=T_H-P{\theta }_H=M_H=0$;
• hinge (H) at H: $M_H=0$;
• horizontal roller (R) at H: $N_H=M_H=0$;
• horizontal slider (S) at H: $N_H=0$.

A linear hyper-elastic law is assumed, linking incremental stresses to strains, in which the elastic constants depend on the prestress, i.e., the elastic operator is $\mathbf{C}=\mathbf{C}\left(\mathring{N}\right)$. Accordingly to what was done in [49], a linear function is chosen, namely (It is important to remark that the linear constitute law (5) can be derived by the Green law, once the existence of a suitable strain energy function, namely $\varphi =\varphi \left(\epsilon ,\gamma ,\kappa ,\mathring{N}\right)$, is postulated):

$$\sigma =\mathbf{C}\left(\mathring{N}\right)\epsilon =\left({\mathbf{C}}^0-P{\mathbf{C}}^{*}\right)\epsilon ,$$

(5)

where $\sigma :={\left(N,T,M\right)}^T$, $\epsilon :={\left(\epsilon ,\gamma ,\kappa \right)}^T$ and:

$$\mathbf{C}:={\mathbf{C}}^0-P{\mathbf{C}}^{*}=\left[\begin{array}{ccc}{c}_{11}& 0& 0\\ & c_{22}& 0\\ \mathrm{Sym}& & c_{33}\end{array}\right]-P\left[\begin{array}{ccc}0& 0& 0\\ & c_{22}^{*}& c_{23}^{*}\\ \mathrm{Sym}& & c_{33}^{*}\end{array}\right].$$

(6)

It is worth noticing that the particular structure of ${\mathbf{C}}^0$ and ${\mathbf{C}}^{*}$ is due to the regularity of the lattice (see, also [49,53]). Moreover, the dependency of the elastic matrix on the prestress entails the existence of two critical values of P, which render it singular, namely $det\left({\mathbf{C}}^0-P{\mathbf{C}}^{*}\right)=0$, the smallest of which is the cell critical load $P_{cell}$. When a cell is loaded at $P=P_{cell}$, a local buckling mode occurs; however, when the cells are assembled, as in the grid beam, they can elastically interact, this entailing that bifurcations is expected to manifests itself at a critical load, $P_c<P_{cell}$.

Another important aspect relies on how the geometric effects are taken into account in the equivalent model. First, as it has been previously done to enforce equilibrium, the adjacent configuration was considered, in order to describe global geometrical effects. Second, the constitutive matrix $\mathbf{C}$, governing Equation (5), also accounts for the prestress to describe local geometrical effects, inside the cell [53].

A displacement formulation of the elastic problem is, finally, obtained by re-writing the equilibrium (4) in terms of displacements, and making use of Equations (1) and (5). The elastic problem reads:

$$\begin{array}{cc}& {\mathbf{K}}_2{\mathbf{v}}^{″}+{\mathbf{K}}_1{\mathbf{v}}^{\prime }+{\mathbf{K}}_0\mathbf{v}=\mathbf{0},\hfill \\ \hfill & {\mathbf{A}}_1{\mathbf{v}}_A^{\prime }+{\mathbf{A}}_0{\mathbf{v}}_A=\mathbf{0},\hfill \\ \hfill & {\mathbf{B}}_1{\mathbf{v}}_B^{\prime }+{\mathbf{B}}_0{\mathbf{v}}_B=\mathbf{0},\hfill \end{array}$$

(7)

where $\mathbf{v}={\left(v,\theta \right)}^T$ and $u=0\forall x$ is accounted for. Moreover, the following definitions are introduced:

$$\begin{array}{cc}\hfill {\mathbf{K}}_2:=& \left[\begin{array}{cc}{c}_{22}-P{c}_{22}^{*}& -P{c}_{23}^{*}\\ -P{c}_{23}^{*}& c_{33}-P{c}_{33}^{*}\end{array}\right],\hfill \\ \hfill {\mathbf{K}}_1:=& \left[\begin{array}{cc}0& -c_{22}+P\left(c_{22}^{*}-1\right)\\ c_{22}-P\left(c_{22}^{*}-1\right)& 0\end{array}\right],\hfill \\ \hfill {\mathbf{K}}_0:=& \left[\begin{array}{cc}0& 0\\ 0& -c_{22}+P\left(c_{22}^{*}-1\right)\end{array}\right],\hfill \end{array}$$

(8)

where ${\mathbf{K}}_i\left(i=0,1,2\right)$ are symmetric or skew-symmetric stiffness matrices. Equation (7) defines a boundary value problem, whose field equation is (7)-a, which is sided by geometric and mechanical boundary conditions at A and B, (7)-b,c. Moreover, the following restraints at the ends are considered, for which the correspondent definitions of boundary operators in Equations (7)-b,c read:

• hinge at A, horizontal roller at B (H-R):

  $$\begin{array}{cccc}& {\mathbf{A}}_1={\mathbf{B}}_1:=\left[\begin{array}{ccc}0& & 0\\ -P{c}_{23}^{*}& & c_{33}-P{c}_{33}^{*}\end{array}\right],\hfill & \hfill & {\mathbf{A}}_0={\mathbf{B}}_0:=\left[\begin{array}{cc}1& 0\\ 0& P{c}_{23}^{*}\end{array}\right];\hfill \end{array}$$

  (9)
• clamp at A, free at B (C-F):

  $$\begin{array}{cccc}& {\mathbf{A}}_1:=\mathbf{0},\hfill & \hfill & {\mathbf{A}}_0:=\mathbf{I},\hfill \\ \hfill & {\mathbf{B}}_1:={\mathbf{K}}_2,\hfill & \hfill & {\mathbf{B}}_0:=\left[\begin{array}{cc}0& -c_{22}+P\left(c_{22}^{*}-1\right)\\ 0& P{c}_{23}^{*}\end{array}\right];\hfill \end{array}$$

  (10)
• clamp at A, horizontal roller at B (C-R):

  $$\begin{array}{cccc}& {\mathbf{A}}_1:=\mathbf{0},\hfill & \hfill & {\mathbf{A}}_0:=\mathbf{I},\hfill \\ \hfill & {\mathbf{B}}_1:=\left[\begin{array}{ccc}0& & 0\\ -P{c}_{23}^{*}& & c_{33}-P{c}_{33}^{*}\end{array}\right],\hfill & \hfill & {\mathbf{B}}_0:=\left[\begin{array}{cc}1& 0\\ 0& P{c}_{23}^{*}\end{array}\right];\hfill \end{array}$$

  (11)
• clamp at A, horizontal slider at B (C-S):

  $$\begin{array}{cccc}& {\mathbf{A}}_1={\mathbf{B}}_1:=\mathbf{0},\hfill & \hfill & {\mathbf{A}}_0={\mathbf{B}}_0:=\mathbf{I}.\hfill \end{array}$$

  (12)

The smallest load $P=P_c$ at which the boundary value problem (7) admits nontrivial solutions triggers incipient instability of the beam; this is referred to as the the global critical load. To evaluate the critical load, Equation (7)-a must be solved. This task can be easily accomplished in closed form by:

• finding the general solution of Equation (7)-a, namely the characteristic exponents ${\lambda }_{1,2}=0$, ${\lambda }_{3,4}=\pm \sqrt{f\left(P\right)}$, with f a function of P, and the associated eigenvectors, which depends on four arbitrary constants;
• enforcing the boundary conditions (7)-b,c, which provide a (homogeneous) system of linear algebraic equations in the four arbitrary constants;
• zeroing the determinant of the matrix of the coefficients of the linear system, which gives a transcendent equation in P, whose smallest root is $P_c$.

3. Elastic and Inertial Constant Identifications

The identification of the elastic constants follows the lines of [53], where it has been analytically carried out by imposing an energy equivalence between the fine, i.e., the frame, and coarse models. It is worth noticing that the analytical identification of the elastic constants, can be pursued when, as in the case at hand, the fibers are arranged in a regular pattern and/or the cross-section is assumed to be rigid (see [45,49,51]). In the remaining cases, a numerical identification procedure must be followed (see, e.g., [52]).

First a discrete map is defined that links at selected abscissas $x_i:=ih\left(i=0,1,\dots ,n\right)$, the displacements of the fine and coarse models. At the $x_i$, $x_{i+1}$ abscissas, the adjacent cross-sections ${\mathcal{D}}_i$, ${\mathcal{D}}_{i+1}$, respectively, bound a cell. It is assumed that the nodal points of the cross-section remain aligned, while the displacements inside the cells are unconstrained. The rigid motion of ${\mathcal{D}}_i$ is written as:

$$\begin{array}{cc}\hfill \tilde{u}\left(x_i,y\right)& =u\left(x_i\right)-\theta \left(x_i\right)y,\hfill \\ \hfill \tilde{v}\left(x_i,y\right)& =v\left(x_i\right),\hfill \end{array}\mathrm{in}{\mathcal{D}}_i,i=0,1,\dots ,n.$$

(13)

where the tilde denotes a quantity of the fine model.

Then the equivalence of the energies stored by a cell and a segment of equal length of the 1D-beam is imposed under the same assigned displacements at the ends. In this step, the averaged energy density is evaluated, namely $\tilde{\varphi }:=\tilde{U}/h$, $\tilde{U}$ being the elastic energy of a cell of the grid beam, which, by virtue of Equation (13), is a function of the configuration variables, i.e., $\tilde{\varphi }=\frac{1}{h}\tilde{U}\left(u\left(x_i\right),v\left(x_i\right),\theta \left(x_i\right),u\left(x_{i+1}\right),v\left(x_{i+1}\right),\theta \left(x_{i+1}\right)\right)$.

Finally, the energy densities $\tilde{\varphi }$ and $\varphi$ are rendered equal. This task requires the evaluation of a suitable strain test field. Here the simplest (constant) one is assumed, and, after integration of the strain–displacement relationships (1) in the $\left(x_i,x_{i+1}\right)$ interval, the end displacements of the beam segment reads:

$$\begin{array}{cc}\hfill u\left(x_i\right)=0,& u\left(x_{i+1}\right)=\epsilon h,\hfill \\ \hfill v\left(x_i\right)=0,& v\left(x_{i+1}\right)=\gamma h+\frac{1}{2}\kappa h^2,\hfill \\ \hfill \theta \left(x_i\right)=0,& \theta \left(x_{i+1}\right)=\kappa h.\hfill \end{array}$$

(14)

The elastic constants of matrix $\mathbf{C}$ follow from $\varphi =\tilde{\varphi }\forall \left(\epsilon ,\gamma ,\kappa \right)$, once use of Equation (14) is made.

It is important to remark that the end displacements (14) can be interpreted as the superposition of three independent deformation modes of the cell, namely:

• an extensional mode ($\mathrm{EX}$), in which ${\mathcal{D}}_{i+1}$ axially translates;
• a shear mode ($\mathrm{SH}$), in which ${\mathcal{D}}_{i+1}$ transversely translates;
• a flexural modes ($\mathrm{FL}$), in which ${\mathcal{D}}_{i+1}$ rotates around the z-axis and transversely translates.

Analytical Identification of the Elastic Constants

A cell made of m x-fibers and two y-fibers of half stiffnesses each, is considered. Only the nodal points on the y-fiber remain aligned (i.e., ${\mathcal{D}}_i$ coincides with the set of the isolated joints), but the same fiber is free to warp around the line.

In order to identify the elastic constants, the procedure previously resumed is applied (see also [53]). It calls for evaluating the elastic energy $\tilde{U}$ of the cell, when the three deformation modes, displayed in Figure 2, are assigned. Since the cell is prestressed, it is $\tilde{U}:={\tilde{U}}^0+{\tilde{U}}^{*}$, ${\tilde{U}}^0$ and ${\tilde{U}}^{*}$ being the elastic and the geometric energetic contributions, respectively.

To apply the procedure, it is also assumed that [51]: (i) in the shear mode $\mathrm{SH}$ all nodes rotate of the same and unknown angle $\phi$; (ii) in the flexural mode $\mathrm{FL}$, the effect of micro-warping is neglected. Then, the energetic contribution of the j-fiber, behaving as a Timoshenko micro-beam, to the cell’s energy, is:

$$\begin{array}{cc}& U_j={\int }_0^h\left(\frac{1}{2}EA{\left({\epsilon }^{\left(1\right)}\right)}^2+\frac{1}{2}G{A}^{*}{\gamma }^2+\frac{1}{2}EI{\kappa }^2+{\mathring{N}}_j{\epsilon }^{\left(2\right)}\right)\mathrm{d}x,\hfill \end{array}$$

(15)

where ${\mathring{N}}_j=-P/m$ for the (prestressed) x-fibers and ${\mathring{N}}_j=0$ for the (not prestressed) y-fibers, respectively. Moreover, for the assigned modes, the energy of each fiber, Equation (15), is evaluated by describing the strain field through the interpolating functions, which exactly satisfy the linear elastic problem of the Timoshenko beam (without prestress). The cell’s elastic and geometric energetic contributions are obtained by summing all the $U_j$; they read:

$$\begin{array}{cc}\hfill {\tilde{U}}^0& :=\frac{1}{2}h\left[mEA{\epsilon }^2+m\frac{12EI}{(1+\alpha )h^2}{\gamma }^2+\left(mEI+EA\sum _{j=1}^m{y}_j^2\right){\kappa }^2\right]\hfill \\ \hfill & +\frac{1}{2}h\left[\frac{12EI(2m-1)}{(1+\alpha )h^2}{\phi }^2-\frac{24EIm}{(1+\alpha )h^2}\gamma \phi \right],\hfill \\ \hfill {\tilde{U}}^{*}& :=-P\frac{1}{2}h\left[\frac{6+10\alpha }{5{(1+\alpha )}^2}{\gamma }^2+\frac{1-5{\alpha }^2}{5{(1+\alpha )}^2}{\phi }^2+\frac{h^2}{3}{\kappa }^2+h\gamma \kappa -\frac{2-10{\alpha }^2}{5{(1+\alpha )}^2}\gamma \phi \right],\hfill \end{array}$$

(16)

where $\alpha :=12EI/\left(G{A}^{*}{h}^2\right)$ is a non-dimensional stiffness ratio. The contributions in the square bracket terms (once multiplied by h) on the right-end sides of Equation (16) are, in the order: (i) in Equation (16)-a the extensional, shear and flexural elastic energies of the x-fibers, and the flexural-shear elastic energy of the y-fibers, respectively; (ii) the flexural-shear geometric contribution of the prestressed x-fibers.

By requiring $\frac{\partial \tilde{U}}{\partial \phi }=0$, the unknown $\phi$ is condensed, and $\tilde{\varphi }$ is known. Then, by imposing $\varphi =\tilde{\varphi },$$\forall \left(\epsilon ,\gamma ,\kappa \right)$ and enforcing the Green law, once $\varphi$ is expanded in Maclaurin series of P up to the first order, the elastic constants, Equation (5), follow:

$$\begin{array}{cc}& c_{11}=mEA,\hfill \\ \hfill & c_{22}=\chi m\frac{12EI}{\left(1+\alpha \right)h^2},\hfill \\ \hfill & c_{33}=\left(mEI+EA\sum _{j=1}^m{y}_j^2\right),\hfill \\ \hfill & c_{22}^{*}={\chi }^{*}\frac{6+10\alpha }{5{\left(1+\alpha \right)}^2},\hfill \\ \hfill & c_{23}^{*}=\frac{h}{2},\hfill \\ \hfill & c_{33}^{*}=\frac{h^2}{3},\hfill \end{array}$$

(17)

where:

$$\begin{array}{cc}& \chi :=\frac{m-1}{2m-1},\hfill \\ \hfill & {\chi }^{*}:=\frac{\left[5\alpha (3\alpha +8)+21\right]m^2-2\left[5\alpha (\alpha +4)+11\right]m+6+10\alpha }{2(5\alpha +3){(1-2m)}^2},\hfill \end{array}$$

(18)

are non-dimensional quantities named the shear factors. However, while $\chi$, accounting for the flexibility of the transverse fibers, is not new, since it has been already encountered in previous works [51,53], ${\chi }^{*}$ is a new corrective term, that takes into account for the flexibility of the transverse fibers also in the geometric part of the constitutive matrix ${\mathbf{C}}^{*}$. Finally, it can be seen that ${\chi }^{*}$ slightly varies with m. Indeed, in the particular case $\alpha =0$, i.e., when micro-beams behave as Euler–Bernoulli beams, it is ${\chi }_s^{*}=\left(21m^2-22m+6\right)/\left[6{(1-2m)}^2\right]$, which ranges from ${\chi }_s^{*}=23/27\simeq 0.85$, when $m=2$, to ${\chi }_s^{*}=7/8=0.875$, when $m\to \infty$.

4. Numerical Results

Numerical results concern two sample grid beams, made of a thermoplastic polymer, having the elastic modulus $E=2180$ $\mathrm{N}/{\mathrm{mm}}^2$ and the Poisson coefficient $\nu =0$. To evaluate the effect of the number of longitudinal fibers m on the buckling response, in the first case study $m=7$ is taken, while in the second it is $m=2$ and in both of them $h=10$ mm. The fibers’ cross-section is squared, of dimensions $1.5\mathrm{mm}\times 1.5\mathrm{mm}$ and its elasto-geometric characteristics are listed in

Table 1. Elasto-geometric properties of the fibers for the case studies I and II.

$\mathbf{EA}$ [$N$]	${\mathbf{GA}}^{*}$ [$N$]	$\mathbf{EI}$ [${N\times \mathrm{mm}}^2$]
4905	$2043.75$	$919.69$

(A is the area, $A^{*}$ the shear area, I the inertia moment). It is worth to notice that the above-mentioned case studies have been selected in order to be representative of grid beams which can be easily printed through additive manufacturing techniques, both for the chosen material and geometric characteristics.

Buckling analyses are carried out by considering different boundary conditions, already defined in Section 2 and corresponding to: hinge–roller (H-R), clamp–free (C-F), clamp–roller (C-R), clamp–slider (C-S). The analytical (exact) solution of the equivalent beam model has been determined according to the procedure sketched in Section 2. Accordingly, critical loads vs. cell number and critical deformed shapes are compared with numerical results obtained by FE models. Concerning FE analyses, the planar grid beam has been modeled as a planar elastic frame, made of Timoshenko micro-beams, eventually restrained at the centroid of the terminal cross-sections, after a rigid body constraints of the same ends is enforced. In the FE discretization procedure here developed: (i) a Timoshenko beam finite-element has been adopted to model the fibers, whose stiffness matrix accounts for both elastic and geometric effects; (ii) each fiber, comprised between two nodes of the planar frame, has been meshed into five finite-elements.

The elasto-geometric characteristics of the equivalent beam model, the cell’s critical loads and the shear factors are listed in

Table 2. Elasto-geometric coefficients of the equivalent beam model, cell’s critical loads and shear factors for the case studies I and II.

		Case Study I	Case Study II
$c_{11}$	$\left[\mathrm{N}\right]$	$3.43\times {10}^4$	$9.81\times {10}^3$
$c_{22}$	$\left[\mathrm{N}\right]$	$3.38\times {10}^2$	$69.81$
$c_{33}$	$\left[\mathrm{N}\times {\mathrm{mm}}^2\right]$	$1.37\times {10}^7$	$2.47\times {10}^5$
$c_{22}^{*}$	$\left[-\right]$	$1.04$	$1.02$
$c_{23}^{*}$	$\left[\mathrm{mm}\right]$	$5.00$	$5.00$
$c_{33}^{*}$	$\left[{\mathrm{mm}}^2\right]$	$33.33$	$33.33$
$P_{cell}$	$\left[\mathrm{N}\right]$	$325.78$	$67.99$
$\chi$	$\left[-\right]$	$0.46$	$0.33$
${\chi }^{*}$	$\left[-\right]$	$0.88$	$0.87$

for both the case studies. Throughout the next figures, the green dots and curves (labeled with A) represent the FE solution, while the continuous blue curves (label B) that of the Timoshenko beam model. Moreover, yellow (label C) and red (label D) curves are relevant to the results given by the Timoshenko and Euler–Bernoulli models, respectively, in which the geometric effect on the constitutive law is neglected, i.e., ${\mathbf{C}}^{*}=\mathbf{0}$.

Buckling analysis of case study I is addressed and the influence of the number of cells is first investigated. Results of this analysis are displayed in Figure 3 for the different boundary conditions described above, namely: H-R (Figure 3a), C-F (Figure 3b), C-R (Figure 3c), C-S (Figure 3d). It is seen that the critical buckling load multiplier $\mu :=P_c/P_{cell}$ decreases with the number of cells n, and, as it is expected, the more constrained the beam, the higher $\mu$ for a given n. Moreover, it is apparent that when $n>\overline{n}$, $\overline{n}$ being a ‘critical’ number of cell, the equivalent model gives excellent approximation of the FE results, i.e., the green curves (A) are almost superimposed to the blue (B) ones. This value $\overline{n}$ depends on how the boundary conditions are applied in the fine model (i.e., restrained end cross-sections modeled as restrained rigid bodies) and it is critical in the sense that below it the homogenization procedure may lead to wrong results. Indeed, the fine model, when $n<\overline{n}$ predicts $\mu >1$, which cannot be reproduced by the equivalent beam model, due to the fact that its constitutive law becomes singular when $P=P_{cell}$. It is also apparent in Figure 3 that, when the geometric effect of the prestress is not taken into account in the elastic law (5), both the Timoshenko and Euler–Bernoulli models, yellow (C) and red (D) curves, are ineffective in describing as the critical load varies with n, even if the former accounts for the flexibility of the transverse fibers via the corrective factor $\chi$, Equation (18). This sort of failure becomes much more evident the more constrained is the beam for a given n. As a final remark, the importance of the shear factors, defined in Equation (18), can be better stressed by looking at their numerical values in Table 2: indeed, it is seen that, not accounting for the transverse fibers flexibility would lead to (elastic and geometric) shear factors equal to 1, which are very different from those here found. Therefore, both the elastic and geometric corrective factors of the constitutive constants (17) play a crucial role in the constitutive law definition, the former producing the much more significant correction.

In order to compare the critical mode of the equivalent model with that furnished by FE discretization on the fine model, $n=50$ is taken. Comparison is shown in Figure 4, where the critical modes of the two models are superimposed, for different boundary conditions, namely: H-R (Figure 4a), C-F (Figure 4b), C-R (Figure 4c), C-S (Figure 4d). It is seen that the green curves, representing axis-line and cross-sections sampled at different abscissas of the Timoshenko equivalent model (remember Equation (13)), are almost coincident with the dots, relevant to the correspondent FE solutions. It is worth mentioning that the critical mode of the C-S grid beam (Figure 4d) resembles the one provided in [16] for the long-wavelength buckling of a square honeycomb with square cells, under uniaxial compression.

Same qualitative considerations are apparent for the buckling analysis of case study II, whose results, in terms of $\mu$ vs. n, are reported in Figure 5 for all the above-mentioned beam models and for the following boundary conditions: H-R (Figure 5a), C-F (Figure 5b), C-R (Figure 5c) and C-S (Figure 5d). Moreover, Figure 6 shows the superimposition of the analytical and numerical critical modes, namely those relevant to the equivalent beam and FE models, respectively, when $n=20$ and for the same boundary conditions. In addition, in this case, the accordance is very good both in terms of critical loads and modes, providing $n>\overline{n}$. Therefore, the effectiveness of the Timoshenko equivalent beam model, corrected by elastic and geometric shear factors, is confirmed also for a low number of longitudinal fibers.

Finally, in order to better investigate how the critical loads depend on m, a buckling analysis on a sample grid beam is carried out. The grid beam’s fibers have the same elasto-geometric characteristics described above and, moreover, $n=50$ and H-R boundary conditions are considered. Results of the buckling analysis are shown in

Table 3. Critical loads of a sample H-R grid beam with $n=50$, for different values of m.

m	${\mathit{P}}_{\mathit{c}}^{\mathbf{EQ}}$ [$N$]	${\mathit{P}}_{\mathit{c}}^{\mathbf{FE}}$ [$N$]	${\mathit{ϵ}}_{\%}$ [$\mathbf{-}$]
2	$8.55$	$8.60$	$0.60$
3	$29.59$	$29.53$	$-0.21$
4	$62.66$	$61.90$	$-1.22$
5	$104.91$	$102.71$	$-2.15$
6	$153.27$	$149.00$	$-2.86$
7	$205.30$	$198.60$	$-3.37$
8	$259.34$	$250.08$	$-3.71$
9	$314.37$	$302.55$	$-3.91$
10	$369.78$	$355.49$	$-4.02$
11	$425.20$	$408.61$	$-4.06$

, where the (dimensional) critical load values are reported for different values of m, namely $m\in \left[2,11\right]$; both the critical buckling loads of the coarse equivalent model, namely $P_c^{EQ}$, and of the refined finite-element model, $P_c^{FE}$, together with the percentage error ${ϵ}_{\%}:=100\left(P_c^{FE}-P_c^{EQ}\right)/P_c^{FE}$ are reported in the same Table 3. It is seen that: (i) the larger is m the larger is the critical load, coherently with the fact that the elastic stiffnesses of the equivalent beam model ($c_{11}$, $c_{22}$ and $c_{33}$ defined in Equation (17)) increase with m; (ii) the magnitude of the percentage error is small, thus confirming the effectiveness of the proposed Timoshenko equivalent beam model.

5. Conclusions and Perspectives

The buckling behavior of planar and uniformly compressed grid beams has been analyzed in this paper. A Timoshenko beam model, embedded in a 2D space, and accounting for the geometric effect of compression of longitudinal fibers, as well as for the micro-warping of the transverse fibers, has been formulated in the framework of a homogenization procedure. Accordingly, elastic and geometric corrective factors have been introduced in the constitutive coefficients, analytically derived through a suitable identification procedure. The buckling analysis of the coarse model, governed by a boundary value problem admitting an exact solution, has been illustrated. The effectiveness of the presented Timoshenko beam model has been discussed on sample grid beams, by comparing the results of exact buckling analyses with benchmark solutions obtained through finite-element discretization of the correspondent bi-dimensional frames.

The following conclusions have been drawn.

• An excellent agreement between the critical loads and modes of the equivalent beam model, with respect to finite-element analyses has been found, providing a sufficient number of cell is considered. This finding is independent from the boundary conditions.
• The geometric effect of the prestress has been introduced both in the equilibrium and constitutive equations of the coarse model. The second aspect, accounting for micro-effects, is shown to be crucial to correctly describe the buckling behavior of the here analyzed micro-structured beams.
• The corrective factors provide a significant correction of the elastic operator, revealing their essential contribution in the proper modeling of planar grid beams.

Finally, the main perspectives of this work are: (i) the extension of the presented model to a 3D space, also including for the presence of micro-structure’s imperfections, which may lead to a flexural-torsional buckling behavior; (ii) the development of buckling experimental tests on 3D-printed grid beams.