Asymptotic homogenization models for pantographic lattices with variable order rotational resistance at pivots

Journal of the Mechanics and Physics of Solids 134 (2020) 103718

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

Asymptotic homogenization models for pantographic lattices with variable order rotational resistance at pivots Nicole Coutris, Lonny L Thompson∗, Sai Kosaraju Department of Mechanical Engineering, Clemson University, United States

a r t i c l e

i n f o

Article history: Received 1 May 2019 Revised 26 August 2019 Accepted 7 September 2019 Available online 8 September 2019 Keywords: Asymptotic homogenization Pantographic lattice Pantographic sheet Meta-materials Second gradient continuum theory Strain gradient continuum theory

a b s t r a c t In this paper, different asymptotic orders of torsional resistance at the pivots for linear pantographic lattices are introduced within an asymptotic homogenization of the discrete micro-beam model associated with the summation of internal forces and moments at the pivots. The small parameter is defined by the ratio of the square root of a unit cell size for the lattice and the problem domain size. The objective is to systematically identify leading order continuum field equations and to define homogenized elasticity tensors and their properties for increasing rotational torque stiffness. To measure the effect of the torsional spring stiffness on the deformation of the lattice, the torsional resistance of a beam element between pivots which compares the torque stiffness to its bending moment stiffness is introduced as a power of the small parameter. Consistent with previous models for pantographic sheets with zero (perfect) or light torsional resistance, an asymptotic expansion with the small scaling parameter shows that the leading order continuum field equations are characterized by second gradient elasticity theory. The homogenized elasticity coefficients are determined in terms of the fiber characteristics and the torque stiffness. A shear modulus is calculated proportional to the torque stiffness. By increasing the torsional resistance to order one, the resulting field equations move to first-gradient classical elasticity with a shear modulus varying as a homographic function in the torsional resistance. By increasing torsional stiffness to higher orders, the shear modulus is obtained for the limiting case of rigid connections. Numerical results are presented for standard elongation and shear bias tests. Strain energies from the first gradient and second gradient terms present for the models with different orders of torsional resistance are computed to demonstrate the zero and infinite limiting cases, and the need to move from a second-gradient continuum model with low-order torsional resistance as scaled with the small parameter, to a first-gradient model for high-orders of torsional resistance. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Pantographic sheets are two-dimensional lattice structures assembled by extensible and continuous fibers that exhibit bending stiffness and interconnected by internal pivots (also called nodes). The nodes constrain the beam fibers to share the same displacement at their connections, however, when modeled as idealized pivots with zero torsional resistance, the intersecting fibers are allowed to rotate independently (Boutin et al., 2017; Rahali et al., 2015). In practice, all physical pivot ∗

Corresponding author. E-mail address: [email protected] (L.L. Thompson).

https://doi.org/10.1016/j.jmps.2019.103718 0022-5096/© 2019 Elsevier Ltd. All rights reserved.

2

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

joints will possess some level of torsional resistance that will store some deformation energy under load. As a result of this torsional resistance at the pivots, homogenized continuum models for the discrete pantographic beam lattice will exhibit shear strain energy; see e.g. Turco et al. (2016), Scerrato et al. (2016) and dell’Isola et al. (2016a). A topic of increasing interest for the design of meta-materials are microstructures with a high-contrast in microscopic properties; see e.g. Vescovo and Giorgio (2014), Hans and Boutin (2008), Boutin and Soubestre (2011), dell’Isola et al. (2012) and Kamotski and Smyshlyaev (2019). Highly heterogeneous microstructures, once homogenized, can produce generalized continuum models such as second-gradient theories (Germain, 1973; Mindlin, 1965); pantographic sheets fall in this category since the rotational resistance between beam fibers is of generally very small order; in the limit of zero rotational resistance, the homogenized shear stiffness is zero for the idealized pivot model; while still exhibiting large second-gradient stiffness behavior due to flexure of the beam fibers under general loading. Standard homogenization techniques used for classical Cauchy first gradient continuum modeling, e.g. classical analytical and computational Representative Elementary Volume (REV) methods when applied to multi-scale materials with high-contrast microstructures may not be adequate for identifying the homogenized macro-properties of the high-contrast meta-materials; see e.g. Biswas and Poh (2017), Kaczmarczyk et al. (2008), Vescovo and Giorgio (2014), Hans and Boutin (2008), Boutin and Soubestre (2011) and dell’Isola et al. (2012). In high-contrast materials, homogeneous deformations may not be sufficient to characterize the material behavior because they are sensitive to higher gradients of the displacement field (Forest, 2013; Forest and Trinh, 2011; Hütter, 2017). Generalized continua, such as second-order gradient models, incorporate a feature of the microstructure which may not be accounted for by standard homogenization methods, namely their size-dependent material response involving intrinsic lengths directly stemming from the microstructure of the material (Forest, 2013). Asymptotic homogenization for multi-scale analysis of finite-scale periodic discrete systems within the framework Cosserat (micro-polar) continuum models, in which an independent rotation field is included have also been studied where the macroscopic behavior depends on the ratio of the characteristic length of constituents and size of heterogeneity with respect to the size of the structure; see e.g. Forest et al. (2001) and Rezakhani and Cusatis (2016). Homogenization of meta-materials with large deformation and pattern transformation may also exhibit nonlocal effects (Rokoš et al., 2019). The homogenized field equilibrium equations and energy of planar pantographic sheets has been studied by several authors. In dell’Isola et al. (2016b), dell’Isola et al. (2016a) and Turco et al. (2016) and elsewhere, discrete models were introduced involving extensional and rotational springs valid for large, finite strains. In Rahali et al. (2015) and Boutin et al. (2017), the lattice has been modeled as a set of beam elements interconnected by internal pivots with zero rotational resistance, valid for small strains. In Rahali et al. (2015), the linearized homogenized model and formulas for effective mechanical properties of pantographic lattices are based on heuristic homogenization procedures described in Mindlin (1965) and dell’Isola et al. (2015). In this analytical method an evaluation of the strain energy density at the microlevel is modeled together with an asymptotic expansion up to the second order. In Boutin et al. (2017), an asymptotic expansion of the discrete beam lattice equilibrium equations, with a perfect pivot modeled with zero torsional resistance, is used to obtain a consistent homogenized continuum model that captures the underlying physics of the beam lattice at different orders. The equations for the continuum model are not obtained assuming an a priori continuum model such as classical first-gradient or micro-polar theories, which can lead to inconsistencies, but instead arise naturally and directly from the equations modeling the physical behavior of the underlying beam lattice. The strains and displacement fields are assumed to occur at a relatively large-scale behavior relative to a characteristic size that is large compared to the size of a unit cell in the lattice. Separation of scales is modeled by introducing a small parameter defined by the ratio of the unit cell area relative to the problem size. The macro-behavior is obtained from the leading asymptotic order as the small scaling parameter tends to zero. Through the asymptotic homogenization process, the continuum field equations and energy are well defined by second gradient continuum theory, as explained in Mindlin (1965) and Germain (1973). Material properties for the homogenized material are explicitly characterized in terms of the beam fibers Young modulus, the fiber cross section properties, (area and moment of inertia), and the distance between nearest pivots (Boutin et al., 2017). As recognized in Boutin et al. (2017) and Eremeyev et al. (2018), and elsewhere, the energy, assuming perfect pivots with zero rotational resistance, is not strongly elliptic in its dependence on second gradients. As a result, the determination of consistent boundary conditions that provide unique solutions is more difficult than found in classical first-gradient elasticity. As physical pivot joints possess some level of torsional resistance, it has been recognized that a shear energy arising from torsion resistance at the pivots should be included within the homogenized continuum modeling for pantographic lattices; see e.g. Turco et al. (2016), Scerrato et al. (2016) and dell’Isola et al. (2016a), and elsewhere. However, in prior work, the shear modulus used in these models has been fitted to data obtained from physical experiments and has not been explicitly defined in terms of the beam stiffness properties and the order of torsional resistance. There is a need to characterize the macro-behavior with different asymptotic orders of torsional resistance scaled with the bending stiffness of the beam fibers between nodes and the scale of the lattice unit cell with respect to the problem size. In this paper, different asymptotic orders of torsional resistance at the pivots for linear pantographic lattices, scaled with the bending stiffness of the beam fibers between nodes, are introduced within an asymptotic homogenization of the discrete micro-beam model associated with the summation of internal forces and moments at the pivots, in order to systematically identify corresponding continuum field equations and associated homogenized material properties. Generalizing the asymptotic homogenization in Boutin et al. (2017), from the discrete beam lattice force and moment equilibrium equations, we include rotational springs which couple the intersecting beam fiber rotations at the pivots, and include them in the asymptotic homogenization. A dimensionless torsional resistance is defined as the ratio of the torque coefficient and beam stiffness. Different cases are

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

3

determined considering the torsional resistance as expressed in terms of powers of a small parameter. The chosen small parameter is the square root of the ratio of the unit cell area and the domain area. Through the asymptotic process, consistent homogenized continuum models that capture the underlying physics of the beam lattice, including rotational resistance at different asymptotic orders, are obtained. The equations for the continuum models generated at the leading order, together with their homogenized material properties, arise naturally and directly from the equations modeling the physical behavior of the underlying beam lattice with the different asymptotic orders in the dimensionless torsional resistance at the pivots. Results from the asymptotic process show the homogenized shear modulus within these equations is defined in terms of different asymptotic orders of torsional resistance at the pivots, together with the beam bending stiffness between pivots. The rigorous asymptotic homogenization procedure shows how the effective behavior (at the macroscopic level) moves progressively from a classical (first gradient) model of elasticity to a second-gradient (strain gradient) model when the rotational stiffness of the pivots (at the microscopic level) decreases. By increasing torsional stiffness to high-orders, the continuum field equations and homogenized shear modulus is obtained for the limiting case of rigid connections. For the purpose of studying the homogenized models with different asymptotic orders of rotational resistance at the pivots, numerical results are presented for standard elongation and shear bias tests. Strain energies obtained with first gradient and second gradient contributions present in the model are computed to demonstrate the limiting cases governed by the torsional resistance and the need to move from a second-gradient continuum model with low-order torsional resistance to a first-gradient model for high-orders of torsional resistance. 2. Governing equations: equilibrium at each node Following Boutin et al. (2017), nodal equilibrium equations for forces and moments for the beam elements for the two fibers at a typical repeated node are written in terms of nodal displacement and rotation variables through stiffness relations. To homogenize these discrete beam lattice stiffness relations, equations are normalized and the displacement and rotation variables from the resulting normalized finite difference stiffness equations are expanded in a Taylor-series with respect to the node and then expanded in an asymptotic expansion with the small scaling parameter ε defined by the square-root of the unit cell area between pivots relative to the size of the problem domain; and after equating powers in ε , the leading order continuum equations corresponding to the discrete lattice equilibrium equations are obtained. 2.1. Nodal forces and moments The stiffness equations for a beam member relating nodal forces and moments to displacements and rotations can be obtained by integrating the static axial force, shear force, and moment equilibrium equations, and evaluating at the two ends. The dimensions of the beam cross-section areas are assumed small relative to the lengths between pivots; and the behavior of the beam segments between pivots are modeled using classical Bernoulli-Euler beam kinematics. For a typical beam member, denote the left end O and right end E. In the local coordinate system for a beam segment, the stiffness relations between (axial, shear) force vector components at O and E: (NO , TO ) and (NE , TE ), and bending moments at O and E: (MO , ME ), with displacement vector components at O and E: (u1,O , u2,O ), (u1,E , u2,E ), and rotations at O and E: (θ O , θ E ), can be expressed as,

NO = −ka (u1,E − u1,O )

 TO = −12k f (u2,E



(1a)

 ( θE + θO ) − u2,O ) − 2

MO = −6k f  (u2,E



 ( θE + θO ) − u2,O ) − 2

(1b)

 − k f  2 ( θE − θO )

NE = ka (u1,E − u1,O )

 TE = 12k f (u2,E



(2a)

 ( θE + θO ) − u2,O ) − 2

ME = −6k f  (u2,E

(1c)



 ( θE + θO ) − u2,O ) − 2

(2b)

 + k f  2 ( θE − θO )

(2c)

It can be shown that the end forces and moments are related by,

NE = −NO ,

TE = −TO ,

ME = TO  − MO .

(3)

In the above,

kf =

EI , 3

ka =

EA , 

(4)

4

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

Fig. 1. Equilibrium at pivot (m, n) at position (x1 , x2 ) = (m1 , n2 ). Two beam members for fiber 1: (B1− , B1+ ) and fiber 2: (B2− , B2+ ), respectively. Moment equilibrium includes torsional springs at nodes coupling fiber rotations Mθ = Kθ (θ2 − θ1 ).

are flexural and axial stiffness measures, where E √ is Young’s modulus, A is the cross-section area, I is the second-area moment of inertia, and  is the beam length. Let h = A be a measure of the thickness of the beam cross-sections, and define the slenderness ratio η = h/. The ratio of the flexural to axial stiffness k f /ka = I/(2 A ) = O(η2 ) is smaller than one, valid for the Bernoulli–Euler kinematics assumption. During the asymptotic process, a uniform or isotropic scaling is applied to the sheet leading to constant stiffness ratios as the small parameter tends to zero (this point is discussed further in Section 2.2). As discussed in the introduction, the pantographic sheet is a two-dimensional lattice modeled by orthogonal continuous fibers that exhibit extensional and bending stiffness, interconnected by internal pivots (also called nodes). The pivots constrain the beam fibers to share the same displacement at their connections, however, when modeled as idealized pivots with zero torsional resistance, the intersecting fibers are allowed to rotate independently. Each fiber is divided into elements located between two consecutive intersections with the other fiber. Each element will be considered as a beam member, the ends of which will be called nodes. Global displacement vector components (u1 , u2 ) at the intersections of the fibers or at the nodes are assumed equal to each other. For node (m, n), the end displacements for the two beam members for fiber 1: (B1− , B1+ ) and the two beam members for fiber 2: (B2− , B2+ ) satisfy the following constraint equations of equal displacements:

u1[m,n] := u1,E (B1− ) = u1,O (B1+ ) = −u2,E (B2− ) = −u2,O (B2+ )

(5a)

u2[m,n] := u2,E (B1− ) = u2,O (B1+ ) = u1,E (B2− ) = u1,O (B2+ )

(5b)

In the above constraint equations, u1[m,n] := u1 (m1 , n2 ), and u2[m,n] := u2 (m1 , n2 ) are defined as displacement component functions evaluated at fiber intersections for a typical repeating node (m, n); see Fig. 1. Rotations of the two fibers at a node are coupled by the resistance of a rotational spring exerting torque on the two fibers. Denote θ1 = θ1[m,n] and θ2 = θ2[m,n] to be the rotations of fiber 1 and fiber 2 at node (m, n). A rotational spring is applied between fiber 1 and fiber 2 at any node. The torque exerted by the rotational spring is:

Mθ = Kθ (θ2 − θ1 )

(6)

Due to the presence of the torque, the two fiber rotations are coupled. The limit Kθ → 0, corresponds to a node where fibers are free to rotate independently with no torsional resistance; this is the case considered in Boutin et al. (2017). The limit Kθ  1, corresponds to a rigid connection with the constraint θ1[m,n] = θ2[m,n] . At node (m, n), the equilibrium equations for the two beam members for fiber 1 (B1− , B1+ ), and two beam members for fiber 2 (B2− , B2+ ) are given below. The force equilibrium equations for the force reactions of the two beams exerted at node (m, n) in the x1 and x2 directions are:

NO (B1+ ) + NE (B1− ) − TO (B2+ ) − TE (B2− ) = 0

(7)

NO (B2+ ) + NE (B2− ) + TO (B1+ ) + TE (B1− ) = 0

(8)

The moment equilibrium equations for the moment reactions exerted from fiber 1 and from fiber 2, together with the moment reaction due to the rotational spring at node (m, n), are

MO (B1+ ) + ME (B1− ) − Kθ (θ2 − θ1 ) = 0

(9)

MO (B2+ ) + ME (B2− ) + Kθ (θ2 − θ1 ) = 0

(10)

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

5

2.2. Stiffness relations between kinematic variables at nodes Substitution of the stiffness equations for end forces and moments in terms of nodal displacements and rotations in the equilibrium equations and enforcing displacement compatibility gives the nodal equilibrium equations in terms of a repeating finite difference stencils for typical node (m, n). Sum of force reactions exerted on node (m, n) in the x1 and x2 directions is zero:



12k f2

 12k f1



δ22 u1[m,n] +

2 (θ2[m,n+1] − θ2[m,n−1] ) + ka1 δ12 u1[m,n] = 0 2

(11)

δ12 u2[m,n] −

1 (θ − θ1[m−1,n] ) + ka2 δ22 u2[m,n] = 0 2 1[m+1,n]

(12)



Sum of moment reactions exerted on node (m, n) from fiber 1 and the rotational spring is zero:





1 (4 + δ12 )θ1[m,n] − k f1 21 δ12 θ1[m,n] 2 − Kθ (θ2[m,n] − θ1[m,n] ) = 0

6k f1 1 −(u2[m+1,n] − u2[m−1,n] ) +

(13)

Sum of moment reactions exerted on node (m, n) from fiber 2 and the rotational spring





2 (4 + δ22 )θ2[m,n] − k f2 22 δ22 θ2[m,n] 2 + Kθ (θ2[m,n] − θ1[m,n] ) = 0

6k f2 2 (u1[m,n+1] − u1[m,n−1] ) +

(14)

In the above equations we have introduced notations for central difference operators defined by,

δ12 φm,n = φm−1,n − 2φm,n + φm+1,n

(15a)

δ22 φm,n = φm,n−1 − 2φm,n + φm,n+1

(15b)

The stiffness properties associated with bending and extension for the orthogonal fibers 1 and 2 are defined by,

k f1 =

E1 I1 , 31

ka1 =

E1 A1 , 1

k f2 =

E2 I2 , 32

ka2 =

E2 A2 , 2

Simplifying, the moment equilibrium Eqs. (13) and (14) for the two fibers at (m, n) can also be expressed as,

−6k f1 1 (u2[m+1,n] − u2[m−1,n] ) + 2k f1 21 (6 + δ12 )θ1[m,n] − Kθ (θ2[m,n] − θ1[m,n] ) = 0

(16)

6k f2 2 (u1[m,n+1] − u1[m,n−1] ) + 2k f2 22 (6 + δ22 )θ2[m,n] + Kθ (θ2[m,n] − θ1[m,n] ) = 0

(17)

When Kθ = 0, these equations match those in Boutin et al. (2017) which considered idealized pivots. Substituting a two-dimensional Taylor series expansion at beam end nodes with respect to a typical node (m, n), for functions, φ[m,n] = φ (m1 , n2 ) representing one of the kinematic variables into the finite difference equations gives the following expressions. For the central difference operators, (15), even powers of derivatives are obtained,

δ12 φ[m,n] = 21

∂ 2φ 2 4 ∂ 4φ 2 6 ∂ 6φ +  +  + ··· ∂ x21 4! 1 ∂ x41 6! 1 ∂ x61

(18)

δ22 φ[m,n] = 22

∂ 2φ 2 4 ∂ 4φ 2 6 ∂ 6φ +  +  + ··· ∂ x22 4! 2 ∂ x42 6! 2 ∂ x62

(19)

while odd derivatives are obtained for the other difference operators appearing in the equilibrium equations at node (m, n):

φ[m+1,n] − φ[m−1,n] = 21

∂φ 2 3 ∂ 3φ 2 5 ∂ 5φ + 1 3 +  + ··· ∂ x1 3! ∂ x1 5! 1 ∂ x51

(20)

φ[m,n+1] − φ[m,n−1] = 22

∂φ 2 3 ∂ 3φ 2 5 ∂ 5φ + 2 3 +  + ··· ∂ x2 3! ∂ x2 5! 2 ∂ x52

(21)

Define L1 and L2 as the lengths of variation of the displacements components u1 and u2 , large compared to 1 and 2 . This study is valid if the two lengths L1 and L2 have the same order of magnitude and if the two lengths 1 and 2 have the same order of magnitude. The small scale parameter is defined by the ratio of area of a unit cell of the lattice relative to the problem domain size;

ε=



( 1 2 )/ ( L1 L2 )

6

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

In the asymptotic homogenization process, as the cell size is decreased (number of lattice cells increased) for a finite size problem domain, ε → 0. The beam lengths scale with

1 = ∗1 ε , L

2 = ∗2 ε , L

(22)



where ∗1 and ∗2 are O(1), and L is chosen as L1 L2 . In order to clearly identify the orders of ε appearing in the asymptotic equations, non-dimensional equilibrium equations are written. Define dimensionless parameters of O(1).

x∗1 =

x1 , L

x∗2 =

x2 , L

u∗1 =

u1 , U1

u∗2 =

u2 , U2

θ1∗ =

θ1 , 2

θ2∗ =

θ2 , 1

(23)

from (16) and (17), it can be deduced that: U1 = U2 = U, and 1 = 2 = , with U = L . The four dimensionless equations ˜ 1 and R ˜ 2 , and K1 and K2 . The choice of the characteristic displacecan now be written with four dimensionless groups R ment U or the characteristic rotation depends on the loading or the boundary conditions. The moment equations lead to: U = L and the introduction of two dimensionless groups which are respectively the ratio of rotational torque stiffness to bending moment stiffness in beam 1 and the ratio of rotational torque stiffness to bending moment stiffness in beam 2 as

K1 = Kθ K2 = Kθ

   1 12E1 I1

   2 12E2 I2

,

(24)

,

(25)

Define the dimensionless ratio of the shear force in beam 2 relative to the normal force in beam 1:

˜ 1 = R1 R

22

, 2

1

R1 =



12k f2 12E2 I2 1 = ka1 E1 A1 32



(26)

Define the dimensionless ratio of the shear force in beam 1 relative to the normal force stiffness in beam 2:

˜ 2 = R2 R

21

, 2

2

R2 =



12k f1 12E1 I1 2 = ka2 E2 A2 31



(27)

√ For beam cross sections and lengths of the same order, let h = A be a measure of the thickness of the beam cross√ sections, with length  = 1 2 . If we consider that Ej , hj and j have the same order of magnitude, then the stiffness ratios ˜ j are O(η2 ). During the asymptotic process, a uniform or isotropic scaling is applied to the sheet leading to constant stiffR ness ratios as the small parameter tends to zero. The length j and thickness hj tend to zero, however, their ratio remains fixed contrary to Boutin et al. (2017) where it was assumed that the slenderness ratio is proportional to ε . If h/l is proportional to ε , it means that h would decrease more than . Working with the assumption in Boutin et al. (2017) for the case of ideal pivots (Case 1 with zero rotational resistance to be discussed in Section 2.3), the effect of the cross-section would be decreased and results in an inextensible macroscopic model that would need a smaller scaling to include microscopic beam extensibility at the macroscopic level. For the case of light rotational resistance (Case 2 to be discussed in Section 2.3), because the effect of the cross-section would be decreased, the scaling with ε on the order of torsional stiffness would also be decreased by O(ε 2 ) which is not physically justified or technically needed. Substituting the Taylor series for the difference operators acting on the kinematic variables in the force equilibrium equation at node (m, n) along x1 axis, gives



˜1 R

  ∗  ∂ 2 u∗1 2 2 ∗2 ∂ 4 u∗1 ∂θ2 1 2 ∗2 ∂ 3 θ2∗ + ε 2 ∗4 + · · · + + ε  + ··· 4! ∂ x∗2 3! 2 ∂ x∗3 ∂ x∗2 ∂ x2 2 2  2 ∗  ∂ u1 2 2 ∗2 ∂ 4 u∗1 + + ε 1 ∗4 + · · · = 0 4! ∂ x∗2 ∂ x1 1

(28)

Similarly, the Taylor series of force equilibrium equation at node (m, n) along x2 axis, for the finite difference stiffness stencil centered at node (m, n)



˜2 R

  ∗  ∂ 2 u∗2 2 2 ∗2 ∂ 4 u∗2 ∂θ1 1 2 ∗2 ∂ 3 θ1∗ + ε 1 ∗4 + · · · − + ε  + ··· 4! ∂ x∗1 3! 1 ∂ x∗3 ∂ x∗2 ∂ x1 1 1  2 ∗  ∂ u2 2 2 ∗2 ∂ 4 u∗2 + + ε 2 ∗4 + · · · = 0 4! ∂ x∗2 ∂ x2 2

(29)

Introducing the dimensionless coefficients K1 and K2 in the moment equations, the equilibrium equations become: moment reactions exerted at node (m, n) from fiber 1 and torsion spring:



   ∂ u∗2 1 2 ∗2 ∂ 3 u∗2 1 2 ∗2 ∂ 2 θ1∗ ∗ − + ε  + · · · + θ1 + ε 1 + · · · − K1 (θ2∗ − θ1∗ ) = 0 ∂ x∗1 3! 1 ∂ x∗3 3! ∂ x∗2 1 1

(30)

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

7

Equilibrium of moment reactions exerted at node (m, n) from fiber 2 and torsion spring:



   ∂ u∗1 1 2 ∗2 ∂ 3 u∗1 1 2 ∗2 ∂ 2 θ2∗ ∗ + ε  + · · · + θ2 + ε 2 + · · · + K2 (θ2∗ − θ1∗ ) = 0 ∂ x∗2 3! 2 ∂ x∗3 3! ∂ x∗2 2 2

(31)

2.3. Asymptotic expansion Following Boutin et al. (2017), by introducing asymptotic expansions and equating powers in the continuum field equations of leading order are obtained. Since the length scales in the Taylor-series and moment nodal equilibrium equations appear in even orders in ε , in order to determine the corresponding to the asymptotic order of torsional resistance and to consider the least degenerate parameters:

K1 = K1∗ ε 2 p , K1∗

K2 = K2∗ ε 2 p

small scaling parameter, expansions for the force leading order equations case, we introduce new

(32)

K2∗

where and are O(1) dimensionless parameters, and asymptotic expansions for the different variables with only even orders. In the above, p is a positive or negative integer used in the asymptotic expansion to define the order of magnitude for the torsion resistance of the rotational spring on the fibers. By changing the integer p, different field equations will be generated corresponding to physical models with different orders of torsional resistance relative to the bending moment stiffness in the beams. The asymptotic expansions for the displacement and rotation kinematic variables are (for α = 1, 2):

uα =

∞

ε 2k uα(2k) = uα(0) + ε 2 uα(2) + ε 4 uα(4) + · · ·

(33)

ε 2k θα(2k) = θα(0) + ε 2 θα(2) + ε 4 θα(4) + · · ·

(34)

k=0

θα =

∞ k=0

Substituting the asymptotic expansions into (28) and (29) and collecting leading order terms give the asymptotic force equations



 ∂ u1 ∗ ( 0 ) ∗ (0 ) + θ2 ∂ x2 ∗ ∂ x2 ∗   ∗ (2 )     ∗ (0 ) ∗ (0 ) 2 ∂ 4 u1 1 ∂ 3 θ2 ∂ ∂ u1 ∗ (2 ) 2 ∗2 4 +ε + θ2 + 2 + + ε [ ] + ... ∂ x∗2 ∂ x∗2 4! ∂ x∗4 3! ∂ x∗3 2 2

2 ∗ (0 )  2 ∗ (2 )   ∗ (0 ) ∂ u1 ∂ u1 2 ∂ 4 u1 2 + + ε + + . . . =0 4! ∂ x∗4 ∂ x∗2 ∂ x∗2 1 1 1

˜1 R

∂



  ∂ ∂ u∗2(0) ∗ (0 ) − θ 1 ∂ x∗1 ∂ x1 ∗   ∗ (2 )     ∗ (0 ) ∗ (0 ) 1 ∂ 3 θ1 2 ∂ 4 u2 ∂ ∂ u2 ∗ (2 ) 2 ∗2 4 +ε − θ1 + 1 − + ε [ ] + ... ∂ x∗1 ∂ x∗1 4! ∂ x∗4 3! ∂ x∗13 1

2 ∗ (0 )  2 ∗ (2 )   ∗ (0 ) ∂ u2 ∂ u2 2 ∂ 4 u2 2 + + ε + + ... = 0 4! ∂ x∗4 ∂ x∗2 ∂ x∗2 2 2 2

(35)



˜2 R

(36)

Similarly, substituting the asymptotic expansions (33) and (34) into (30) and (31) and collecting leading order terms give the asymptotic moment equations,





 ∗ (2 )  3 ∗ (0 )   ∂ u∗2(0) ∂ u2 ∂ u2 ∂ 2 θ1 ∗ ( 0 ) ∗ (0 ) ∗ (2 ) 2 ∗2 1 4 − θ1 + ε − θ1 +  1 − + ε [ ] + ... ∂ x∗1 ∂ x∗1 3! ∂ x∗2 ∂ x∗3 1 1 ∗(0) 

 ∗ 0 ∗ 2 = −K1 θ2 − θ1 ( ) + ε 2 θ2∗(2) − θ1 ( ) + . . .

(37)

 ∗ (2 )  3 ∗ (0 )   ∂ u∗1(0) ∂ u1 ∂ u1 ∂ 2 θ2 ∗ ( 0 ) ∗ (0 ) ∗ (2 ) 2 ∗2 1 4 + θ + ε + θ +  + + ε + . . . [ ] 2 2 2 ∂ x∗2 ∂ x∗2 3! ∂ x∗2 ∂ x∗3 2 2 ∗(0) 

 ∗ 0 ∗ 2 = −K2 θ2 − θ1 ( ) + ε 2 θ2∗(2) − θ1 ( ) + . . .

(38)

8

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

Case 1: Zero torsion resistance In this first case, p > 1 is chosen such that K1 and K2 , do not appear in first two leading order terms in the asymptotic expansion for the moment Eqs. (37) and (38). The resulting equations specialize to the case of zero torsional resistance for perfect pivots between beam fibers 1 and 2 at nodes, and derived by Boutin et al. (2017). To characterize the leading order terms, consider the equations necessary to obtain the displacements u∗1(0 ) and u∗2(0 ) , and rotations θ1∗(0 ) and θ2∗(0 ) . If the leading order equations couple displacements, they can be first-gradient or secondgradient homogenized continuum models; if the equations couple displacements and rotations, the model can be expressed in terms of micropolar theory. Equating the leading powers of ε in the asymptotic expansion for the moment Eqs. (37) and (38) gives the equations at order 0 and order 2:

∂ u∗2(0) − θ1∗(0) = 0, ∂ x∗1 ∂ u∗1(0) + θ2∗(0) = 0, ∂ x∗2

(39)

 ∗ (0 )  2 ∂ u2 ∂ u∗2(2) ∗ (2 ) ∗ (0 ) ∗2 1 ∂ − θ1 +  1 − θ1 = 0, ∂ x∗1 3! ∂ x∗2 ∂ x∗1 1  ∗ (0 )  2 ∂ u1 ∂ u∗1(2) ∗ (2 ) ∗ (0 ) ∗2 1 ∂ + θ +  + θ = 0, 2 2 2 ∂ x∗2 3! ∂ x∗2 ∂ x∗2 2

(40)

Taking into account the Eqs. (39) and (40) obtained at order 0 and 2, (35) and (36), can now be written as:



˜1 R



∗ (0 )

∂ 4 u1 2 − ε 2 ∗2 2 4! ∂ x∗4 2

˜2 − R

∗ (0 )

2 2 ∗2 ∂ 4 u2 ε 1 4! ∂ x∗4 1





 ∂ 2 u∗1(0) + ··· + + ··· = 0 ∂ x∗2 1   2 ∗ (0 )  ∂ u2 + ··· + + · · · =0 ∂ x∗2 2

(41)

(41)

It can be deduced that the tension in beams 1, respectively beams 2, is constant. The beams are inextensible at the order zero. This solution corresponds to a traction or a compression applied to the sheet in the direction of the fibers, which cannot be considered as the solution for general boundary conditions. It was assumed that L was the scale of variation of u1 and u2 in both directions. Suppose now that L remains the scale of variation of u1 in the direction orthogonal to fibers1 and u2 in the direction orthogonal to fibers 2, but there exist L1 and L2 , length scales of variation associated to u1 and u2 in the direction of fibers 1 and fibers 2. From (11) and (12), comparing the orders of magnitude of the different terms, it can be deduced: L1 = L2 = L/ε . With these different length scales, it is possible to capture the microscopic beam extensibility behavior at the macroscopic level (Boutin et al., 2017; dell’Isola and Steigmann, 2015; Rahali et al., 2015). For Case 1 and Case 2 of the pantographic sheet with zero or small rotational resistance at nodes, the change of variables needed is defined as:

u∗1(0) (x∗1 , x∗2 ) → u¯ ∗1(0) (x¯∗1 , x∗2 ),

x¯∗1 = ε x∗1

(42a)

u∗2(0) (x∗1 , x∗2 ) → u¯ ∗2(0) (x∗1 , x¯∗2 ),

x¯∗2 = ε x∗2

(42b)

∂ u∗1(0) ∂ u¯ ∗1(0) =ε ∗ ∂ x1 ∂ x¯∗1

(43a)

∂ u∗2(0) ∂ u¯ ∗2(0) = ε ∂ x∗2 ∂ x¯∗2

(43b)

The effect of this rescaling is to couple in the equations from the asymptotic process, the second-order derivatives (first gradient, axial strains) with the fourth-order derivatives (second gradient, orthogonal bending strains) in the leading order asymptotic force equations. Because of this coupling, the asymptotic equations capture the extensible behavior of the microscopic model. For the case of torsion resistance of order one (Case 3) and the limit of rigid connection (Case 4), only first gradient terms appear in asymptotic force equations and the change of variable is not needed. ˜ j are O(η2 ) and remain fixed in the asymptotic process as ε tends to zero. In As discussed earlier, the R Boutin et al. (2017), it was assumed that the slenderness ratio h/ is proportional to ε . With the effect of the beam crosssection being reduced in the asymptotic limit, the change of variables is now given by

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

9

u∗1(0) (x∗1 , x∗2 ) → u¯ ∗1(0) (x¯∗1 , x∗2 ),

x¯∗1 = ε 2 x∗1

(44a)

u∗2(0) (x∗1 , x∗2 ) → u¯ ∗2(0) (x∗1 , x¯∗2 ),

x¯∗2 = ε 2 x∗2

(44b)

˜ j is considered fixed, in order to include microscopic beam extensibility at the macroscopic level. In the following, since R this new scaling is not needed. Substituting the equations of lower asymptotic order into the equations for higher-order in the moment equations and substituting these equations into the leading order asymptotic expansions for the force Eqs. (35) and (36), with (42) leads to the following dimensionless equations

˜ 1 ∗2 R 2

∗ (0 ) 2 ∂ 4 u1 = 4! ∂ x∗4 2

˜ 2 ∗2 R 1

2 ∂ 4 u2 4! ∂ x∗4 1

∗ (0 )

=

∂ 2 u¯ ∗1(0) , ∂ x¯∗2 1

(45a)

∂ 2 u¯ ∗2(0) ∂ x¯∗2 2

(45b)

From this result, the corresponding dimensional field equations corresponding to extensible beams derived in Boutin et al. (2017) are obtained

E2 I2 1

∂ 4 u1(0) E1 A1 ∂ 2 u1(0) = , 2 ∂ x42 ∂ x21

(46a)

E1 I1 2

∂ 4 u2(0) E2 A2 ∂ 2 u2(0) = 1 ∂ x41 ∂ x22

(46b)

For this model, secondary field variables of rotation and shear strain are defined by,

∂ u2(0) ∂ u (0 ) , θ2(0) = − 1 ∂ x1 ∂ x2  ∗ (0 )  ∂ u1 ∂ u∗2(0) (0 ) (0 ) (0 ) γ = − ( θ2 − θ1 ) = + ∂ x∗2 ∂ x∗1 θ1(0) =

(47)

(48)

The presence of 4th-order derivatives in (46) shows that the second gradient elasticity theory is needed to describe the homogenized beam system. The energy expression associated with these equilibrium equations are given in Boutin et al. (2017). With idealized pivots with zero torsional resistance, these force equilibrium equations do not couple shear strain when the material coordinates are aligned with the orthogonal beam fiber orientations. Case 2: Light torsion resistance In this case, we set p = 1 such the nondimensional torsional resistances K1 = K1∗ ε 2 , K2 = K2∗ ε 2 , is small and now appears in the second leading order equations in the asymptotic expansion for the moment Eqs. (37) and (38). The resulting equations model light torsional resistance between beam fibers 1 and 2 at nodes of the order ε 2 . Equating the leading powers of ε in the asymptotic expansion for the moment equations and substituting the equations of lower asymptotic order into the equations for higher-order lead to

ε0 :

∂ u∗2(0) − θ1∗(0) = 0, ∂ x∗1 ∂ u∗1(0) + θ2∗(0) = 0, ∂ x∗2

ε2 :

∂ u∗2(2) − θ1∗(2) = −K1∗ (θ2∗(0) − θ1∗(0) ), ∂ x∗1 ∂ u∗1(2) + θ2∗(2) = −K2∗ (θ2∗(0) − θ1∗(0) ), ∂ x∗2

Substituting this result into the leading order asymptotic expansions for the force Eqs. (35) and (36) leads to the existence of different length scales in different directions as in Case 1. Since the torsional resistance is small, to model both beams as extensible at the macrolevel, and to match orders of first and second gradients in the asymptotic force equations, the change of scale in (42) is used. The dimensionless equations read now:

10

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

 ˜ 1 K2∗ R

 ˜ 2 K1∗ R

∂ 4 u∗1(0) ∂ ∗ (0 ) ∗ (0 ) ∗2 2 ( θ − θ ) +  2 1 ∂ x∗2 2 4! ∂ x∗4 2 ∗ (0 ) ∂ 2 ∂ 4 u2 (θ2∗(0) − θ1∗(0) ) + ∗2 1 ∗ ∂ x1 4! ∂ x∗4 1

Substituting,

θ

∗ (0 ) 2

−θ

∗ (0 ) 1



∂ u∗1(0) ∂ u∗2(0) =− + ∂ x∗2 ∂ x∗1

 =

∂ 2 u¯ ∗1(0) , ∂ x¯∗2 1

=

∂ 2 u¯ ∗2(0) ∂ x¯∗2 2



 (49)

the dimensionless field equations for displacements of zero-order are:



   ∂ 4 u∗1(0) ∂ ∂ u∗1(0) ∂ u∗2(0) ∗2 2 + + 2 = ∂ x∗2 ∂ x∗2 ∂ x∗1 4! ∂ x∗4 2   ∗ (0 )   ∂ u∗2(0) ∂ 4 u∗2(0) ∂ ∂ u1 ∗2 2 ˜ 2 −K1∗ ∗ R + +  = 1 ∂ x1 ∂ x∗2 ∂ x∗1 4! ∂ x∗4 1 ˜ 1 −K2∗ R

∂ 2 u¯ ∗1(0) , ∂ x¯∗2 1 ∂ 2 u¯ ∗2(0) ∂ x¯∗2 2

The corresponding dimensional equations for displacements are then,

  ∂ 2 u1(0) ∂ ∂ u1(0) ∂ u2(0) + β + = ∂ x2 ∂ x2 ∂ x1 ∂ x21   (0 ) E2 A2 ∂ 2 u2 ∂ ∂ u1(0) ∂ u2(0) +β + = 1 ∂ x1 ∂ x2 ∂ x1 ∂ x22 E1 A1 2

E2 I2 1

∂ 4 u1(0) , ∂ x42

(50a)

E1 I1 2

∂ 4 u2(0) ∂ x41

(50b)

with shear modulus

β=

Kθ 1 2

(51)

for the homogenized sheet. For this model, secondary field variables of rotation and shear strain are defined by (47) and (48). With the addition of torsion resistance, shear strains are coupled with normal strains in the equilibrium Eq. (50). Due to the presence of 4th-order derivatives, the force equilibrium equations can also be expressed in terms of a second-gradient continuum elasticity theory. By introducing the non-dimensional parameter of torque resistance with bending moment resistance in the asymptotic expansions, the relationship between the order of magnitude of torsion resistance Kθ , the homogenized shear modulus β and the other material parameters is established. This result is consistent with the shear modulus obtained from linearization of finite strain shear energy models of torsion resistance given in Turco et al. (2016), Scerrato et al. (2016) and dell’Isola et al. (2016a). Shear modulus was not explicitly defined in terms of the beam stiffness properties and the order of torsional resistance. With K1∗ = O(1 ) and K2∗ = O(1 ) both with the same order of magnitude K1∗ /K2∗ = O(1 ), then the magnitude for torsion resistance for the equations of Case 2 is O(ε 2 ):

Kθ = max

 12E I 12E I  1 1 2 2 1

,

2

K∗ ε 2 ,

(52)

where K∗ of O(1) is specified, ε = max(1 , 2 )/L is a small parameter, and L is the smaller dimension of the problem domain size. As a result, the shear modulus is O(ε 2 ) and can be expressed as:

β=





Kθ 12E1 I1 12E2 I2 K∗ 2 = max , ε 1 2 1 2 1 2

(53)

Alternatively, if the different quantities (1 , 2 ) and/or (I1 , I2 ), (E1 , E2 ), have the same order of magnitude, i.e, 1 /2 = O(1 ), I1 /I2 = O(1 ), E1 /E2 = O(1 ), then the magnitude for torsion resistance for the equations of Case 2 can be chosen as,

Kθ =



12E1 I1 1

 12E I  2 2 2

K∗ ε 2

(54)



where the small parameter is defined by the ratio of areas is ε = (1 2 )/(L1 L2 ). When E = E1 = E2 and I = I1 = I2 , and  = 1 = 2 , the shear modulus normalized with bending moment stiffness is

β=

12EI ∗ 2 K ε 3

(55)

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

11

Case 3: Torsion resistance of order 1 Continuing with the progression of increasing torsional resistance, we move to p = 0 such the size of the nondimensional torsional resistance K1 = K1∗ = O(1 ), and K2 = K2∗ = O(1 ), such that the torsional resistances appear in the first leading order asymptotic equations. The resulting equations model torsional resistance between beam fibers 1 and 2 at nodes of order 1. Equating the leading power in the asymptotic expansion for the moment Eqs. (37) and (38) gives

∂ u∗1(0) + θ2∗(0) = −K2 (θ2∗(0) − θ1∗(0) ), ∂ x∗2

(56a)

∂ u∗2(0) − θ1∗(0) = −K1 (θ2∗(0) − θ1∗(0) ), ∂ x∗1

(56b)

Substituting this result into the asymptotic expansions for the leading terms in the force Eqs. (35) and (36) leads to the following dimensionless equations



 ∂ 2 u∗1(0) ∂ ∗ (0 ) ∗ (0 ) ( θ − θ ) = , 1 ∂ x∗2 2 ∂ x∗2 1   ∂ 2 u∗2(0) ∂ ˜ 2 K1∗ ∗ (θ2∗(0) − θ1∗(0) ) = R ∂ x1 ∂ x∗2 2 ˜ 1 K2∗ R

(57a)

(57b)

From the leading order moment equations, (56), solving for the zero-order rotations

θ2∗(0) − θ1∗(0) = −

1 (1 + K1∗ + K2∗ )



∂ u∗1(0) ∂ u∗2(0) + ∂ x∗2 ∂ x∗1



(58)

and substituting into (57) leads to the following dimensionless equations for displacements of zero-order:

˜ 1 K2∗ R ∂ − (1 + K1∗ + K2∗ ) ∂ x∗2 ˜ 2 K1∗ R ∂ − (1 + K1∗ + K2∗ ) ∂ x∗1

 

∂ u∗1(0) ∂ u∗2(0) + ∂ x∗2 ∂ x∗1 ∂ u∗1(0) ∂ u∗2(0) + ∂ x∗2 ∂ x∗1



=

∂ 2 u∗1(0) , ∂ x∗2 1

(59a)

=

∂ 2 u∗2(0) ∂ x∗2 2

(59b)



The corresponding dimensional equations are,

  ∂ 2 u1(0) ∂ ∂ u1(0) ∂ u2(0) + β + = 0, ∂ x2 ∂ x2 ∂ x1 ∂ x21   (0 ) E2 A2 ∂ 2 u2 ∂ ∂ u1(0) ∂ u2(0) + β + =0 1 ∂ x1 ∂ x2 ∂ x1 ∂ x22 E1 A1 2

(60a)

(60b)

with homogenized shear modulus

β=

Kθ (1 + K1 + K2 )−1 1 2

(61)

or simplifying

β=

1 2

  1 12E1 I1

1 +

2 1 + 12E2 I2 Kθ



(62)

In this case, the shear modulus varies as a homographic function with linear fractional behavior in Kθ > 0. For this model, secondary field variables of rotation and shear strain are defined by

θ1(0) =

∂ u2(0) K1 − γ (0 ) ∂ x1 ( 1 + K1 + K2 )

θ2(0) = − γ

(0 )

(63a)

∂ u1(0) K2 + γ (0 ) ∂ x2 ( 1 + K1 + K2 ) (0 )

= −(1 + K1 + K2 )(θ2

(0 )

− θ1 ) =

(63b)



∂ u1(0) ∂ u2(0) + ∂ x2 ∂ x1

 (64)

12

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

In this case, the force equilibrium equations for displacement field variables are expressed in terms of second-order derivatives and are compatible with classical elasticity theory. With K1∗ /K2∗ = O(1 ), the magnitude for torsion resistance for the equations of Case 3 is of order 1,

Kθ = max

 12E I 12E I  1 1 2 2 ,

1

2

K∗

(65)

where K∗ of O(1) is specified. As a result, the shear modulus in (62) is O(1) and varies as a homographic function with respect to torsional resistance. Alternatively, 1 /2 = O(1 ), I1 /I2 = O(1 ), E1 /E2 = O(1 ), then the magnitude for torsion resistance for the equations of Case 3 can be chosen as,

Kθ =



 12E I  2 2

12E1 I1 1

2

K∗

(66)

Case 4: Rigid Connection In this case, we set p = −1 such the nondimensional torsional resistance K1 = K1∗ ε −2 and K2 = K2∗ ε −2 are large. Equating the leading order terms in the asymptotic expansion for the moment Eqs. (37) and (38) gives

0 = −K1∗ (θ2∗(0) − θ1∗(0) ),

(67a)

0 = −K2∗ (θ2∗(0) − θ1∗(0) ),

(67b)

∂ u∗2(0) − θ1∗(0) = −K1∗ (θ2∗(2) − θ1∗(2) ) ∂ x∗1

(68a)

∂ u∗1(0) + θ2∗(0) = −K2∗ (θ2∗(2) − θ1∗(2) ) ∂ x∗2

(68b)

and

From (67), the fiber rotations of order-zero for the two beams are equal corresponding to a rigid connection,

θ ∗(0) = θ2∗(0) = θ1∗(0) Solving (68) for

θ ∗ (0 ) = gives,

1 1 + rb∗

 rb∗

∂ u∗2(0) ∂ u∗1(0) − ∂ x∗1 ∂ x∗2

 (69)

 ∗ (0 )  ∂ u1 ∂ u∗2(0) ∂ u∗2(0) 1 ∗ (0 ) − θ = + 1 ∂ x∗1 (1 + rb∗ ) ∂ x∗2 ∂ x∗1  ∗ (0 )  rb∗ ∂ u1 ∂ u∗1(0) ∂ u∗2(0) ∗ (0 ) + θ = + 2 ∂ x∗2 (1 + rb∗ ) ∂ x∗2 ∂ x∗1

(70a)

(70b)

where

rb∗ =

K2∗ K2 E1 I1 2 = = K1∗ K1 E2 I2 1

(71)

is a non-dimensional ratio of bending moment stiffness in beam 1 relative to beam 2. In the special case where the two beams have the same material and geometric properties, E1 = E2 , I1 = I2 , and 1 = 2 , then rb∗ = 1. Substituting (70) into the leading order terms for the asymptotic expansion of the force Eqs. (35) and (36) leads to the dimensionless equations for displacements of order-zero,

  ∂ ∂ u∗1(0) ∂ u∗2(0) + + (1 + rb∗ ) ∂ x∗2 ∂ x∗2 ∂ x∗1   1 ∂ ∂ u∗1(0) ∂ u∗2(0) ˜ R2 + + (1 + rb∗ ) ∂ x∗1 ∂ x∗2 ∂ x∗1 ˜1 R

rb∗

∂ 2 u∗1(0) =0 ∂ x∗2 1

(72a)

∂ 2 u∗2(0) =0 ∂ x∗2 2

(72b)

The corresponding dimensional equations are,

E1 A1 2

  ∂ 2 u1(0) ∂ ∂ u1(0) ∂ u2(0) +β + = 0, ∂ x2 ∂ x2 ∂ x1 ∂ x21

(73a)

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

  ∂ 2 u2(0) ∂ ∂ u1(0) ∂ u2(0) + β + =0 ∂ x1 ∂ x2 ∂ x1 ∂ x22

E2 A2 1

13

(73b)

with homogenized shear modulus

β=

1 2

  1

1

12E1 I1

+

2 12E2 I2



(74)

or simplifying

β=

12(E1 I1 )(E2 I2 ) 1 2 (E1 I1 2 + E2 I2 1 )

(75)

This form shows that the homogenized shear modulus is a function of the bending stiffness E1 I1 /1 and E2 I2 /2 for the beams. The coefficient also follows from (62) when considering Kθ infinitely large. The limit does not depend on Kθ . In the special case where the two beams have the same material and geometric properties, E = E1 = E2 , I = I1 = I2 , and  = 1 = 2 , then, β = 6k f = (12EI/3 )/2, or β = 6EI/3 . This result shows that the homogenized shear is directly proportional to beam bending. From this result, we find that in the case where the normalized torsional resistance is large of order ε −2 , where ε = max(1 , 2 )/L is small,

K1 = K1∗ ε −2 ,

K2 = K2∗ ε −2

the torsional resistance Kθ , disappears from the equations. In addition, a new dimensionless number introduced in (71) defined by the ratio of bending moment stiffness in beam 1 relative to beam 2, can be used to study the behavior of rigidly connected frames. Similar to the previous case, the force equilibrium equations for the displacement field variables are expressed in terms of second-order derivatives and are compatible with classical elasticity theory. In the case of rigid connection the first gradient model (69) arises provided that the stiffness of the beams are of the same order, while if the orthogonal beams have very high contrast, the combination of terms of different orders can lead to a second gradient model. For the rigid case, the rotation field variable obtained in (69) can be expressed as

θ

(0 )

=

(0 )

where,

(0 ) = −

1 2

1 (1 − rb∗ ) + 2 (1 + rb∗ )





∂ u1(0) ∂ u2(0) + ∂ x2 ∂ x1

∂ u1(0) ∂ u2(0) − ∂ x2 ∂ x1



(76)

 (77)

is the infinitesimal rotation of a differential element. When the two orthogonal beams have the same material and geometric properties, rb∗ = 1, and θ (0 ) = (0 ) . We remark that it can also be shown that for p < −1, corresponding to normalized torsional stiffness of even larger order, carrying out the asymptotic expansion leads to the same equations. For the predominant terms the homogenized field equations for displacements (73) with homogenized material properties derived using the asymptotic expansion of the moment equations with torsional resistance of order ε −2 or larger, and without assuming a priori θ1 = θ2 , match the field equations obtained with the a priori assumption that the beam fibers are connected rigidly such that beam rotations are constrained with θ = θ1 = θ2 . To show this, assume the beam fibers 1 and 2 are connected rigidly a priori, at a node with θ[m,n] = θ1[m,n] = θ2[m,n] . With this constraint, the corresponding moment equilibrium equation at node ([m, n]), is obtained from the sum of (13) and (14),





1 (4 + δ12 )θ[m,n] − k f1 21 δ12 θ[m,n] 2   2 + 6k f2 2 (u1[m,n+1] − u1[m,n−1] ) + (4 + δ22 )θ[m,n] − k f2 22 δ22 θ[m,n] = 0 2

6k f1 1 −(u2[m+1,n] − u2[m−1,n] ) +

(78)

Dividing (78) with 12k f2 22 = 12E2 I2 /2 , and writing the Taylor-series for the resulting stiffness difference stencil, the nondimensional equilibrium equation summing moment reactions from both beam fiber 1 and 2, exerted at node (m, n) can be expressed as,



   ∂ u∗2 1 2 ∗2 ∂ 3 u∗2 6 2 ∗2 ∂ 2 θ ∗ ∗ + ε  + · · · + + ε  + · · · θ 1 ∂ x∗1 3! 1 ∂ x∗3 4! ∂ x∗2 1 1

 ∗    2 ∗ ∂ u1 1 2 ∗2 ∂ 3 u∗1 6 ∂ θ ∗ 2 ∗2 + + ε  + ··· + θ + ε 2 + ··· =0 ∂ x∗2 3! 2 ∂ x∗3 4! ∂ x∗2 2 2

rb∗ −

(79)

14

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

where rb∗ is the bending stiffness ratio defined earlier in (71). Applying the asymptotic expansions for the field variables to (79) gives,



 ∗ (2 )    ∗ (0 ) ∂ u∗2(0) ∂ u2 1 ∂ 2 θ ∗ (0 ) 1 ∂ 3 u2 ∗ (0 ) 2 ∗ (2 ) ∗2 4 − θ + ε − θ +  − + ε + . . . [ ] 1 ∂ x∗1 ∂ x∗1 3! ∂ x∗3 3! ∂ x∗2 1 1

∗ (0 )  ∗ (2 )    ∗ (0 ) ∂ u1 ∂ u1 1 ∂ 3 u1 1 ∂ 2 θ ∗ (0 ) ∗ (0 ) 2 ∗ (2 ) ∗2 4 + +θ +ε +θ + 2 + + ε [ ] + .... = 0 ∂ x∗2 ∂ x∗2 3! ∂ x∗3 3! ∂ x∗2 2 2

−rb∗

(80)

Equating the leading order term in the asymptotic expansion for this moment equation gives



−rb∗

  ∗ (0 )  ∂ u∗2(0) ∂ u1 ∗ (0 ) ∗ (0 ) − θ + + θ =0 ∂ x∗1 ∂ x∗2

Solving this equation we find,

θ ∗ (0 ) =

1 1 + rb∗



rb∗

∂ u∗2(0) ∂ u∗1(0) − ∂ x∗1 ∂ x∗2

(81)

 (82)

With θ ∗(0 ) = θ1∗(0 ) = θ2∗(0 ) , then

 ∗ (0 )  ∂ u1 ∂ u∗2(0) ∂ u∗2(0) 1 ∗ (0 ) − θ1 = + ∂ x∗1 (1 + rb∗ ) ∂ x∗2 ∂ x∗1  ∗ (0 )  rb∗ ∂ u1 ∂ u∗1(0) ∂ u∗2(0) ∗ (0 ) + θ2 = + ∂ x∗2 (1 + rb∗ ) ∂ x∗2 ∂ x∗1

(83a)

(83b)

which matches Eq (70) obtained earlier when considering a nondimensional torsional resistance of order ε −2 or larger. Substituting (83) into the leading order terms in the asymptotic expansion of the force Eqs. (35) and (36) with θ ∗(0 ) = θ1∗(0) = θ2∗(0) leads to the same leading order homogenized Eq. (73) for displacements of order-zero. 3. Behavior of the homogenized shear modulus for different orders of torsional resistance To further examine the behavior of the different orders of torsional resistance with ε for the different cases in the asymptotic homogenization, consider the case where the two beams have the same properties E = E1 = E2 and I = I1 = I2 , and  = 1 = 2 . With this simplification, define a non-dimensional shear modulus by the ratio of the shear modulus with the bending moment stiffness,

β ∗ = β 2

   12EI

Moving from increased torsional resistance orders for the different cases, the shear modulus formulas derived from the asymptotic expansions specialize to,

β ∗ = f α ( K ∗ )ε 2 p (K ∗ )

where fα ( f 3 = (2 +

(84)

is a function of the dimensionless torsional resistance of order 1. For Cases α = 2, 3: ( f2 = p = 0 ). For Case 1 β ∗ = 0, and for the rigid case 4, β ∗ = 1/2. K∗

1 −1 K∗ ) ,

K∗ ,

p = 1 ),

• Case 2 ( p = 1) The shear modulus normalized with bending moment stiffness and the small parameter ε 2 is a linear function in K∗ = O(1 ).

β ∗ ε −2 = K∗

(85) K∗ .

As a result, the shear strain energy will increase proportionally to increasing • Case 3 ( p = 0) The shear modulus normalized with bending moment stiffness is a homographic function in K∗ = O(1 ) > 0.

β∗ =

K∗ 1 = 1 + 2K ∗ 2 + K1∗

(86)

In the limit of K∗ → ∞,

lim

K∗ →∞

β∗ =

1 2

(87)

asymptotes to a constant value which matches the rigid case. As a result, the shear strain energy will asymptote to the rigid case.

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

15

Fig. 2. Normalized Shear Modulus vs. Non-dimensional Torsion Resistance K∗ .

• Case 4 (Rigid) The shear modulus normalized with bending moment stiffness is constant.

β∗ =

1 2

(88)

Fig. 2 graphs the nondimensional shear modulus for the different cases in Eqs. (85) and (86) normalized with order ε 2p vs. K∗ . Case 4 rigid is obtained from Case 3 with K∗ → ∞. Case 1 is obtained from Case 2 with K∗ = 0. The Case 2 shear expression, scaled with ε 2 , included in the second-gradient equations, is the linear approximation of the shear expression for Case 3 at its origin and the Case 2 shear expression degenerates from the shear expression in the classical first-order equations for Case 3. However, increasing torsional resistance in the second-order gradient equations for Case 2 (p = 1) will not asymptote to the rigid case. To asymptote to the rigid case, we must move from Case 2 equations to the first-order gradient equations for Case 3 (p = 0) and increase the torsional resistance to infinity. 4. Homogenized field equations As discussed earlier, the presence of 4th-order derivatives for Case 1 and Case 2 shows that the second gradient elasticity theory is needed to describe the homogenized beam system. When moving to torsional resistance of order 1 and rigid connections at nodes, for Case 3 and Case 4, the force equilibrium equations for displacement variables are expressed in terms of second-order derivatives and are compatible with classical elasticity theory. The complete set of field equations of equilibrium for second gradient theory can be traced back to Mindlin and Eshel (1968). In Germain (1973), it is shown that second gradient theory may be obtained from the theory of micromorphic media of order 1 when a material point is considered as a classical continuum of small extent, with the gradients of the relative displacements (velocities) of particles, assumed to have the same deformation as the general macroscopic continuum. When reference coordinates are oriented with material coordinates aligned with the orthogonal fibers with basis vectors e1 and e2 , the field equations for displacement vector components u1 (x1 , x2 ) and u2 (x1 , x2 ) derived from the asymptotic expansions are:

  ∂ 2 u1 ∂ ∂ u1 ∂ u2 ∂ 4 u1 + C + − D1 = 0, 12 2 ∂ x2 ∂ x2 ∂ x1 ∂ x1 ∂ x42   ∂ 2 u2 ∂ ∂ u1 ∂ u2 ∂ 4 u2 C2 + C + − D2 =0 12 2 ∂ x1 ∂ x2 ∂ x1 ∂ x2 ∂ x41 C1

C1 =

E1 A1 , 2

C2 =

E2 A2 , 1

C12 = β ,

D1 =

E2 I2 , 1

(89a)

(89b)

D2 =

E1 I1 2

where C12 = β is the shear modulus defined for the different orders of torsional resistance as discussed earlier in the asymptotic expansions. For the two-dimensional model, the dimensions of the homogenized material properties C1 , C2 and C12 are force per length, while D1 , D2 are units force times length (units of moment). As discussed earlier in the development of the homogenized equations from the asymptotic expansions: For Case 1: the shear modulus C12 = 0 and corresponds to the equations derived in Boutin et al. (2017); For Case 2: C12 is positive and defined in (51) with Kθ > 0. For Case 3: C12 is positive and defined in (62) with Kθ > 0. For Case 4: C12 is positive and defined in (75). In all cases, C1 > 0, and C2 > 0. For cases 1 and 2, D1 > 0 and D2 > 0. For cases 3 and 4, D1 = D2 = 0. We introduce the following notations: Symmetric Cauchy stress (rank-2) tensor,

σ = σi j ei  e j This classical stress tensor is related through the constitutive relations to first-order gradients.

(90)

16

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

When reference coordinates are oriented with material coordinates aligned with the orthogonal fibers, the Cauchy stress components are defined in terms of first-derivatives of displacements,

σ11 = C1 u1,1 , σ22 = C2 u2,2 , σ12 = C12 γ12 = 2 C12 (u2,1 + u1,2 )

(91)

Here commas are used to denote partial derivatives with respect to coordinates. A hyperstress (rank-3) tensor is defined by the second gradient of the displacement.

τ˜ = τ˜pqm e p  eq  em When reference coordinates are oriented with material coordinates aligned with the orthogonal fibers, the only nonzero hyperstress tensor components are

τ˜221 = D1 u1,22 , τ˜112 = D2 u2,11 , other τ˜i jk = 0

(92)

With these definitions, (89) can be written

∂σ11 ∂σ12 ∂ 2 τ˜221 + − =0 ∂ x1 ∂ x2 ∂ x22

(93a)

∂σ12 ∂σ22 ∂ 2 τ˜112 + − =0 ∂ x1 ∂ x2 ∂ x21

(93b)

Following Boutin et al. (2017), the equilibrium equations in term of displacement (89) and (93) can be expressed in tensor form consistent with the field equilibrium equations for second-order gradient theory (Germain, 1973):

∇ · (σ − ∇ · τ˜ ) = 0

(94)

or ∇ · s = 0, where the intrinsic stress s = σ − ∇ · τ˜ is defined in terms of the symmetric tensor σ and the divergence of the hyperstress tensor τ˜ , symmetric with respect to the first two indices. The weak-form and boundary conditions associated with the equilibrium equations for the displacement field vector derived from the asymptotic expansions, can be determined by weighting the residual with the vector function δ u and integrating over the domain D



D

{∇ · (σ − ∇ · τ˜ )} · δ u da = 0

(95)

Applying the divergence (Gauss) theorem twice to move the gradients to δ u gives



D

(σ · ·∇δ u + τ˜ · · · ∇∇δ u ) da =



∂D

[n · σ − n · (∇ · τ˜ )] · δ u ds +



∂D

(n · τ˜ ) · ·∇δ u ds

(96)

In the above, we have introduced the notations:

σ · ·∇ u = σ pm

∂ um ∂ xp

τ˜ · · · ∇ ∇ δ u = τ˜pqm

(97)

∂ 2 um = τ˜pqm κ pqm = τ˜ · · · κ ∂ x p ∂ xq

(n · τ˜ ) · ·∇ u = (n · τ˜ ) pm

(98)

∂ um ∂ xp

(99)

where the convention of repeated indices implies summation. The boundary conditions are typically prescribed in local (n, t) coordinates, where n = (n1 , n2 ) and t = (t1 , t2 ) = (−n2 , n1 ) are the unit normal and tangent vectors to the boundary ∂ D, respectively. Cartesian and normal-tangent derivatives are related by the transformation rule for the partial derivatives,

∂ um ∂ um ∂ um = ni + ti , ∂ xi ∂n ∂s

∂ um ∂ um = ni , ∂n ∂ xi

∂ um ∂ um = ti ∂s ∂ xi

Working with the last boundary integral in (96)

 ∂D

(n · τ˜ ) · ·∇δ u ds =

 ∂D

[n · (n · τ˜ )] ·

∂ δ u ds + ∂n

 ∂D

[t · (n · τ˜ )] ·

∂ δ u ds ∂s

(100)

The boundary is composed of a number of smooth curves with possible discontinuity at vertex points. Integrating-by-parts the tangent derivative on the boundary

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

 ∂D

[t · (n · τ˜ )] ·

∂ δ u ds = − ∂s

 ∂D

nv ∂ [[t · (n · τ˜ ) · δ u(sv )]] [t · (n · τ˜ )] · δ u ds + ∂s v=1

17

(101)

where sv are boundary vertex points. Combining gives

 D

(σ · · · δε + τ˜ · · · δκ ) da =

 ∂D

Q n · δ u ds +

 ∂D

Mn ·

nv ∂ δ u ds + [[M nt (sv ) · δ u(sv )]] ∂n v=1

(102)

where δε is the symmetric part of ∇δ u. The traction on the smooth segments of boundary ∂ D, in terms of stress tensors are

Qn = Tn − Vn −

∂ M on ∂ D ∂ s nt

M n = n · (n · τ˜ )

on

(103)

∂D

(104)

where Qn is a traction force vector per unit length, T n − V n is the intrinsic stress tensor projected on the normal to the boundary,

T n = n · σ,

V n = n · (∇ · τ˜ )

(105)

and M n = n · (n · τ˜ ) is defined on the boundary ∂ D excluding any boundary vertex points. The jump in Mnt occurs at boundary vertices sv on ∂ D (if any) where ∂ D is discontinuous,

[[M nt (sv )]] = [[t · (n · τ˜ )]],

for sv ∈ v

(106)

4.1. Boundary conditions Following Germain (1973), from the boundary integrals of the weak form, essential boundary conditions are identified as

u = u¯ on

∂D

(107)

∂ u = u¯ n on ∂ D ∂n

(108)

where the normal derivative on smooth parts of the boundary ∂ D excluding any vertex points sv , of the vector field u(x) is,

∂ ∂ um u = ( n · ∇ )u = ni e ∂n ∂ xi m

(109)

Natural (traction) or dual boundary conditions are:

Q n = Q¯ n ¯n Mn = M

on on

∂D

(110)

∂D

¯ nt (sv )]] = h¯ (sv ), [[M nt (sv )]] = [[M

(111) for sv ∈ v

(112)

where the bar denotes prescribed values for the external forces and moments applied on ∂ D. For second-gradient theory, ¯ n , and M ¯ nt (sv ) will not appear septhe hyper-stress tensor couples both Qn and Mn ; as a result, the external loads Q¯ n , M arately; instead they are interrelated (Germain, 1973). In the above, Q¯ n is a boundary traction force per unit length, and ¯ n is a moment per unit length (moment intensity) referred in the literature to as a prescribed normal double traction M ¯ nt at boundary vertices with discontinu(Germain, 1973), and h¯ (sv ) are any concentrated force vectors due to a jump in M ous tangents. For Case 1 and Case 2, to obtain unique solutions to the weak-form, Dirichlet (essential) boundary conditions must be specified at points on ∂ D in a manner that will ensure zero energy modes are removed from the solution space. For boundary edges that are parallel with the orthogonal fibers, restrictions on the boundary conditions for unique solutions are discussed in Andreaus et al. (2016) and Eremeyev et al. (2018). Cases 3 and 4 with β > 0, and τ˜ = 0, specialize to the classical first-gradient elasticity theory, and only the displacement boundary condition (107) and/or dual (110) is specified. Existence and uniqueness results for Case 3 and 4 within classical elasticity are proven using standard techniques (Brenner and Scott, 1994; Ciarlet and Lions, 1991).

18

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

4.2. Internal strain energy The internal strain (deformation) energy is

U (u ) =

1 2



D

(σ · · · ε + τ˜ · · · κ ) da

(113)

The deformation energy can be divided into the sum of contributions from first-order gradient axial and shear strains, and bending (second-gradient) strains. When reference coordinates are oriented with material coordinates aligned with the orthogonal fibers,

σ11 = C1 ε11 , σ22 = C2 ε22 , σ12 = 2 C12 ε12 = C12 γ12

(114)

τ˜221 = D1 κ221 , τ˜112 = D2 κ112 , other τ˜i jk = 0

(115)

and the different contributions to strain energy can be expressed as the sum of first-order and second-order gradient terms,

U (u ) = (Ua + Us ) + Ub

(116)

where the contributions from the first-order gradient strains are

Ua =

1 2

Us =

1 2



D

 D



2 2 C1 ε11 + C2 ε22 da =

2 C12 γ12 da =

1 2

 D

1 2



 D

1 2



D

 da

β (u1,2 + u2,1 )2 da

and second-gradient bending strains

Ub =

E1 A1 2 E2 A2 2 u + u 2 1,1 1 2,2



2 2 D1 κ221 + D2 κ112 da =

1 2

 D

E I 2 2 1

(117) (118)

u21,22 +



E1 I1 2 u da 2 2,11

(119)

For Case 1, C12 = β = 0, and Us = 0 and the energy matches the expressions given in Boutin et al. (2017) and Placidi et al. (2017). For Case 3 and Case 4 (rigid), D1 = D2 = 0, and Ub = 0 reverts to classical first-gradient elasticity. For Case 3 and Case 4 with Ub = 0, a positive shear modulus β > 0 is required for the strain energy to be coercive and positive definite in the energy space for classical first-gradient elasticity. As recognized by Boutin et al. (2017) and Eremeyev et al. (2018), and others, energy for Case 1 with zero shear modulus β = 0, is neither elliptic (semi-definite) nor strongly elliptic (positive-definite) in its dependence on second gradients. However, as pointed out in Eremeyev et al. (2018), the energy forms a semi-norm in anisotropic Sobolev spaces (Besov et al., 1978a; 1978b) but is elliptic in the sense of Agmon–Douglis–Nirenberg (Agmon et al., 1959; 1964). Using these spaces, it is shown that the strain energy is coercive and positive definite in the natural energy space for second-gradient theory. When considering positive shear modulus in Case 2, the corresponding energy is elliptic and similar to classical elasticity (Case 3 and Case 4 with positive shear modulus), considering Dirichlet boundary conditions sufficient to eliminate rigid body motions (zero energy states), the Korn inequality can be used to prove coercivity, and using the Lax–Milgram theorem, existence and uniqueness of the weak solutions can be proven (Eremeyev et al., 2018). 5. Numerical results Similar to Boutin et al. (2017) and Placidi et al. (2017), and others, we model standard axial and shear bias tests for a rectangular domain with a 45◦ orientation with respect to the orthogonal fibers; i.e., in the undeformed configuration, the orthogonal beam fibers are aligned in the horizontal and vertical directions, not parallel to the rectangle edges. The homogenized material properties are orthotropic and exhibit highly anisotropic behavior with respect to different coordinate rotations. This test problem has also been used in Placidi et al. (2016) and dell’Isola et al. (2016a) and others for large strain simulations and comparison to physical experiments. For the purpose of studying the homogenized models with different asymptotic orders of rotational resistance at the pivots, the numerical results are obtained by a discretization with C1 -continuous elements. In this work, we minimize the energy using complete bicubic Hermitian polynomial rectangle elements (F.K. Bogner and Schmit, 1965) together with compatible triangle elements (Gileva et al., 2013) to fit the polygonal domain boundary and to ensure C1 continuity of the field variable and first-derivatives across inter-element boundaries. Use of these elements provides numerical results which correlate well with others (Boutin et al., 2017; Scerrato et al., 2016), and allows for a comparison of the homogenized models with different orders of rotational resistance at the pivots. Numerical results are also consistent with the nonlinear discrete beam lattice (micro) and homogenized continuum (macro) models in Placidi et al. (2016), dell’Isola et al. (2016a), dell’Isola et al. (2016b) and Turco et al. (2016). For the numerical studies, we model the rectangular domain with aspect ratio L2 /L1 = 3; where L1 = 101 and L2 = 302 with beam lengths between nodes  = 1 = 2 = 4.95 mm. The beam section properties for the two beam fibers are

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

19

equal and assumed to be rectangular with a slenderness ratio of η = h/ = 1/10 and b = 1.6 mm. The material properties for the beam fibers are E = E1 = E2 = 1600 MPa. With this geometry, A = bh, I = bh3 /12, 12k f = 12E I/3 = E bη3 , and the homogenized shear modulus for the different orders of torsional resistance specializes to • Case 1: β = 0, • Case 2: β /(Eb) = η3 K∗ ε 2 , K∗ • Case 3: β /(Eb) = η3 1+2 K∗ • Case 4: β /(Eb) =

1 3 2η

Recall that Case 2 is obtained from the asymptotic expansion with torsional resistance Kθ = 12k f K∗ ε 2 and due to the 4th-order derivatives in the homogenized field equations is modeled with second-gradient theory. Case 1 corresponding to zero torsional resistance at the pivot nodes is equivalent to Case 2 with K∗ = 0 and has zero shear strain energy. For the numerical results, the value of the small parameter ε from the asymptotic expansion is taken as ε = /L1 = 1/10. Moving to Case 3, Kθ = 12k f K∗ , the homogenized shear modulus is order one with K∗ > 0 and modeled with first-gradient classical elasticity and zero bending energy. Case 4 modeling the beam lattice with rigid connections is equivalent to Case 3 with K∗ = ∞; in this limit, the shear modulus does not depend on a torsional resistance value and is constant. 5.1. Bias elongation test Fig. 3 graphs deformation energies for Case 1 and Case 2 models over a short and longer range of nondimensional torsional resistance√K∗ obtained from the bias test at 45◦ orientation with essential boundary conditions defining elongation with magnitude 1/ 2 on the short edge boundary with zero normal derivatives. As K∗ increases from zero (corresponding to Case 1 with zero shear modulus β ), the shear modulus increases linearly, and we observe both the shear energy Us and total energy increases. In contrast, the bending second-gradient energy Ub and first-order gradient axial strain energies Ua increase at a much lower rate. For this problem, when the dimensionless torsional resistance K∗ approaches about 0.1, the shear energy begins to exceed the second-gradient bending energy. The results also show the first-gradient axial energy is very small compared with the second-gradient bending energy. As expected, the Case 2 second-gradient model with K∗ → ∞, and β → ∞ does not match Case 4 corresponding to rigid beam connections, modeled with classical first-gradient equations which do not include any second-gradient bending energy. Fig. 4 graphs deformation energies for Case 3 and Case 4 models obtained from the bias elongation test. Recall that these cases model order one and infinite torsional resistance at pivots of the beam lattice, respectively, and are solved with

Fig. 3. Bias elongation test strain energies for Case 2. Case 1 is equivalent to Case 2 with K∗ = 0.

Fig. 4. Bias Elongation Test Strain Energies for Case 3. Case 4 (Rigid) is equivalent to Case 3 with K∗ = ∞.

20

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

Fig. 5. Bias Shear Test Strain Energies. Case 1 is equivalent to Case 2 with K∗ = 0. Case 4 (Rigid) is equivalent to Case 3 with K∗ = ∞.

the homogenized equations and boundary conditions within classical first-order gradient theory. For Case 3, as K∗ becomes increasingly large, the shear modulus and energies asymptote to Case 4 (Rigid). Conversely, as K∗ approaches zero, the shear modulus within the classical first-gradient theory approaches zero, resulting a singularity with zero energy (loss of ellipticity) in the limit as expected. For Case 3, the axial energy is significantly smaller than the shear energy for this bias test. Since the bending second-gradient energy is zero for Case 3, there is no torsional resistance value for this model that will match a solution from Case 2 or Case 1. Fig. 5 graphs the energies for Case 2 and Case 3 for the bias test at 45◦ orientation with a shear boundary condition. In this shear loading, the deformation energies are an order of magnitude smaller than the energies from the axial loading, but the same trends with increasing torsional resistance for both Case 2 and moving to Case 3 are similar to the trends observed for the elongation bias test. Fig. 6 shows contours representing the distribution of total strain energy for the bias elongation test problem superimposed on the classical deformed geometry with necking shape. The strain energy distribution is computed from the area weighted average of integrated element energies at shared vertices. Fig. 6a, shows Case 1 with K∗ and β equal to zero. Consistent with Boutin et al. (2017), the energy is concentrated within distinct bands of (internal) boundary layers contrasting with triangular zones with much smaller energy. Fig. 6b shows how these energy bands spread within the central region of the domain for Case 2 modeled with the second-gradient equations with a small shear modulus, corresponding to a small increase in torsional resistance at the pivots. As discussed earlier from the energy graphs in Fig. 3, when the dimensionless torsional resistance increases to K∗ = 1 and K∗ = 3; see Fig. 6c and 6d, the shear energy becomes larger than the bending and axial energies and manifests as a homogenous zone in the center region of the problem domain. Moving from the Case 2 model to the Case 3 model with different equations and boundary conditions result in different solutions. For Case 3, as K∗ becomes increasing large, the magnitudes and distribution of energies asymptote to Case 4 (Rigid). By trial and error, a K∗ value can be found such that the total energy for Case 2 and 3 matches, however, the solutions obtained from the different equations and boundary conditions for these cases are not the same. This is illustrated when comparing Case 3 with a small K∗ in Fig. 6g with Case 2 with K∗ = 3 in Fig. 6d; in both these cases, the total energies computed over the domain are the same, but the energy distribution is very different. Similarly, a K∗ for Case 2 can be defined such that the total energy matches the energy obtained for Case 4 rigid, however, the solutions are not the same due to the second-gradient shear energy present in Case 2 which is not present in the solution for Case 4. This is illustrated when comparing Case 2 with a large K∗ selected in Fig. 6e to match the total energy for Case 4 in Fig. 6h; while the energy distribution in the contours for total energy look similar, they are not the same since Case 4 does not have any secondgradient bending energy and the shear modulus is obtained in the limit of K∗ → ∞ for Case 3. Indeed, if K∗ is increased further for Case 2, the shear and total energy completely overshoot the energies for Case 4 with the distribution of energy dominated by a large homogeneous shear energy zone; this is illustrated in Fig. 6f for Case 2 with a large dimensionless torsional resistance of K∗ = 10 0 0. These results illustrate that for a large shear modulus in the Case 2 second-gradient model, not only is the strain energy much larger that the rigid Case 4 model, the overall deformed geometry is changed significantly. These results demonstrate the zero and infinite limiting cases, and the need to move from a second-gradient continuum model with low-order torsional resistance as scaled with the small asymptotic parameter, to a first-gradient model for highorders of torsional resistance. 5.2. Bias shear test In Fig. 7 contours representing the distribution of total strain energy for the shear bias test problem are shown for the different cases. Consistent with Boutin et al. (2017), K∗ and β equal to zero, the energy is concentrated within distinct

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

Fig. 6. Total Strain Energy Distribution for the Bias Extension Test.

21

22

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

Fig. 7. Total Strain Energy Distribution for Bias Shear Test.

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

23

Fig. 8. Strain Energy Distributions for Bias Axial Test for Case 2 with K∗ = 0.05.

bands of (internal) boundary layers. Similar to the elongation bias test, the energy bands spread for Case 2 modeled with the second-gradient equations with a small shear modulus, from a small increase in torsional resistance at the pivots of order ε 2 . When the dimensionless torsional resistance increases to K∗ = 1 and K∗ = 3, the shear energy becomes larger than the bending and axial energies and manifests as a homogenous zone. Moving from Case 2 to Case 3 models with different equations and boundary conditions result in different solutions. For Case 3, as K∗ becomes increasing large, the magnitudes and distribution of energies asymptote to Case 4 (Rigid). Similar to the elongation bias test, by trial and error, a K∗ value can be found such that the total energies for Case 2 and Case 3 match, however, the solutions obtained from the different equations and boundary conditions for these cases are not the same. In Figs. 8 and 9, contours are shown for the different strain energies for Case 2 with small torsional resistance of order ε2 , for the elongation and shear bias tests, respectively. For these results, the dimensionless K∗ values are small such that shear strain energy is smaller than the second-gradient strain energy. For both the elongation and shear tests, the firstorder (axial) energy is relatively small compared to the shear and bending energy, consistent with physical observation of the micro-beams where the axial deformation is significantly less than bending deformation for pantographic lattices. The second-gradient (bending) energy is concentrated within distinct bands of (internal) boundary layers, while the shear energy distribution is characterized with homogenous zones in the middle for the elongation bias test, and at the long edges for the shear test.

24

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

Fig. 9. Strain Energy Distributions for Bias Shear Test for Case 2 with K∗ = 0.2.

6. Conclusions Asymptotic homogenization of pantographic lattices constructed from orthogonal beam fibers interconnected by pivots with internal moments defined by a torsional stiffness coupling the two beam fibers has been performed to obtain the leading order macro continuum equations characterizing the micro behavior. In this this paper, the asymptotic homogenization developed in Boutin et al. (2017) for Case 1 (with idealized pivots modeled with zero torsional resistance), has been generalized to include different asymptotic orders of dimensionless torsional resistance in terms of a small geometric scaling parameter ε , defined by the ratio of a unit cell size in the lattice relative to the problem domain size. Through the asymptotic homogenization process, including a torsional resistance, the leading order continuum field equations and associated homogenized shear modulus properties were systematically identified. Consistent with previous models for pantographic lattices with zero (Case 1), or light torsional resistance (Case 2), the asymptotic expansion shows that for torsional resistance of low-order ε 2 , scaled with the bending stiffness of the beam segments between pivots, the corresponding continuum field equations show second gradient effects involving a second gradient elasticity theory, and the shear modulus for the homogeneous sheet. By increasing the torsional resistance to order one, scaled with the bending stiffness of the beams within the asymptotic expansions, the resulting field equations (Case 3) move to first-gradient classical elasticity with homogenized shear modulus varying as a homographic function in the torsional resistance. By increasing torsional stiffness to higher-orders (Case 4), the continuum field equations and shear modulus for the homogeneous sheet is obtained for the limiting case of rigid beam connections. To examine the

N. Coutris, L.L. Thompson and S. Kosaraju / Journal of the Mechanics and Physics of Solids 134 (2020) 103718

25

behavior of the homogenized models with different asymptotic orders of rotational resistance, numerical results were presented for standard elongation and shear bias tests. Strain energies from the first-order gradient (axial and shear) and second-order gradient terms present for the models with different orders of torsional resistance were computed which illustrate the behavior at the zero and infinite limiting cases, and the need to move from a second-gradient continuum model with low-order torsional resistance as scaled with the small asymptotic parameter, to a first-gradient model for high-orders of torsional resistance. References Agmon, S., Douglis, A., Nirenberg, L., 1959. Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions, I. Communications on Pure and Applied Mathematics 12 (4), 623–727. Agmon, S., Douglis, A., Nirenberg, L., 1964. Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions, II. Communications on Pure and Applied Mathematics 17 (1), 35–92. Andreaus, U., dell’Isola, F., Giorgio, I., Placidi, L., Lekszycki, T., Rizzi, N.L., 2016. Numerical simulations of classical problems in two-dimensional (non) linear second gradient elasticity. International Journal of Engineering Science 108, 34–50. doi:10.1016/j.ijengsci.2016.08.003. Besov, O., II’in, V., Nilol’skii, S., 1978a. Integral Representations of Functions and Imbedding Theorems, 1. Wiley. Besov, O., II’in, V., Nilol’skii, S., 1978b. Integral Representations of Functions and Imbedding Theorems, 2. Wiley. Biswas, R., Poh, L., 2017. A micromorphic computational homogenization framework for heterogeneous materials. Journal of the Mechanics and Physics of Solids 102, 187–208. Boutin, C., dell’Isola, F., Giorgio, I., Placidi, L., 2017. Linear pantographic sheets: Asymptotic micro-macro models identification. Mathematics and Mechanics of Complex Systems 5 (2), 127–162. doi:10.2140/memocs.2017.5.127. Boutin, C., Soubestre, J., 2011. Generalized inner bending continua for linear fiber reinforced materials. International Journal of Solids and Structures 48 (3-4), 517–534. Brenner, S.C., Scott, L.R., 1994. The Mathematical Theory of Finite Element Methods. Springer-Verlag. Ciarlet, P., Lions, J., 1991. Handbook of Numerical Analysis, v. II. Finite Element Methods (Part 1). North-Holland, Amsterdam, New York, Oxford. dell’Isola, F., Andreaus, U., Placidi, L., 2015. At the origins and in the vanguard of peridynamics, non-local and higher-gradient continuum mechanics: An underestimated and still topical contribution of Gabrio Piola. Mathematics and Mechanics of Solids 20 (8), 887–928. doi:10.1177/1081286513509811. dell’Isola, F., Corte, A.D., Giorgio, I., Scerrato, D., 2016a. Pantographic 2D sheets: Discussion of some numerical investigations and potential applications. International Journal of Non-Linear Mechanics 80, 200–208. doi:10.1016/j.ijnonlinmec.2015.10.010. dell’Isola, F., Giorgio, I., Pawlikowski, M., Rizzi, N.L., 2016b. Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472 (2185), 20150790. doi:10.1098/rspa.2015.0790. dell’Isola, F., Seppecher, P., Madeo, A., 2012. How contact interactions may depend on the shape of cauchy cuts in nth gradient continua: approach ‘à la d’alembert”. Journal of Applied Mathematics and Physics / Journal de Mathématiques et de Physique appliquées 63 (6), 1119–1141. dell’Isola, F., Steigmann, D., 2015. A two-dimensional gradient-elasticity theory for woven fabrics. Journal of Elasticity, 1 (118), 113–125. Eremeyev, V.A., dell’Isola, F., Boutin, C., Steigmann, D., 2018. Linear pantographic sheets: Existence and uniqueness of weak solutions. Journal of Elasticity 132 (2), 175–196. doi:10.1007/s10659- 017- 9660- 3. F.K. Bogner, R.F., Schmit, L., 1965. The generation of interelement compatible stiffness and mass matrices by the use of interpolation formulas. In: Proceedings of the Conference on Matrix Methods in Structural Mechanics. Wright-Patterson Air Force Base, Ohio, pp. 397–444. Forest, S., 2013. Generalized Continua from the Theory to Engineering Applications. CISM International Centre for Mechanical Sciences (Courses and Lectures), H Altenbach and V A Eremeyev (Eds), 541. Springer, Vienna doi:10.1007/978- 3- 7091- 1371- 4_5. Micromorphic Media Forest, S., Pradel, F., Sab, K., 2001. Asymptotic analyis of heterogeneous cosserat media. International Journal of Solids and Structures 38, 4585–4608. Forest, S., Trinh, D., 2011. Generalized continua and non-homogeneous boundary conditions in homogenisation methods. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 91 (2), 90–109. doi:10.10 02/zamm.2010 0 0109. Germain, P., 1973. The method of virtual power in continuum mechanics. Part 2: Microstructure. SIAM Journal on Applied Mathematics 25, 556–575. Gileva, L., Shaydurov, V., Dobronets, B., 2013. The triangular Hermite finite element complementing the Bogner-Fox-Schmit rectangle. Applied Mathematics 4 (No. 12A), 50–56. doi:10.4236/am.2013.412A006. Hans, S., Boutin, C., 2008. Dynamics of discrete framed structures: a unified homogenized description. Journal of Mechanics of Materials and Structures 3 (9), 1709–1739. doi:10.2140/jomms.2008.3.1709. Hütter, G., 2017. Homogenization of a cauchy continuum towards a micromorphic continuum. Journal of the Mechanics and Physics of Solids 99, 394–408. doi:10.1016/j.jmps.2016.09.010. Kaczmarczyk, L., Pearce, C.J., Bicanic, N., 2008. Scale transition and enforcement of rve boundary conditions in second-order computational homogenization. International Journal for Numerical Methods in Engineering 74 (3), 506–522. Kamotski, I.V., Smyshlyaev, V.P., 2019. Bandgaps in two-dimensional high-contrast periodic elastic beam lattice materials. Journal of the Mechanics and Physics of Solids 123, 292–304. Mindlin, R., 1965. Second gradient of strain and surface-tension in linear elasticity. International Journal of Solids and Structures 1 (4), 417–438. doi:10. 1016/0 020-7683(65)90 0 06-5. Mindlin, R., Eshel, N., 1968. On first strain gradient theories in linear elasticity. International Journal of Solids and Structures 4, 109–124. Placidi, L., Andreaus, U., Giorgio, I., 2017. Identification of two-dimensional pantographic structure via a linear D4 orthotropic second gradient elastic model. Journal of Engineering Mathematics 103 (1), 1–21. doi:10.1007/s10665- 016- 9856- 8. Placidi, L., Barchiesi, E., Turco, E., Rizzi, N.L., 2016. A review on 2D models for the description of pantographic fabrics. ZAMP - Zeitschrift für angewandte Mathematik und Physik 67 (5), 121. doi:10.10 07/s0 0 033- 016- 0716- 1. Rahali, Y., Giorgio, I., Ganghoffer, J., dell’Isola, F., 2015. Homogenization à la Piola produces second gradient continuum models for linear pantographic lattices. International Journal of Engineering Science 97, 148–172. doi:10.1016/j.ijengsci.2015.10.003. Rezakhani, R., Cusatis, G., 2016. Asymptotic expansion homogenization of discrete fine-scale models with rotational degrees of freedom for the simulation of quasi-brittle materials. Journal of the Mechanics and Physics of Solids 88, 320–345. doi:10.1016/j.jmps.2016.01.001. Rokoš, O., Ameen, M., Peerlings, R., Geers, M., 2019. Micromorphic computational homogenization for mechanical metamaterials with patterning fluctuation fields. Journal of the Mechanics and Physics of Solids 123, 119–137. doi:10.1016/j.jmps.2018.08.019. The N.A. Fleck 60th Anniversary Volume Scerrato, D., Zhurba Eremeeva, I.A., Lekszycki, T., Rizzi, N.L., 2016. On the effect of shear stiffness on the plane deformation of linear second gradient pantographic sheets. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 96 (11), 1268– 1279. doi:10.10 02/zamm.20160 0 066. Turco, E., dell’Isola, F., Cazzani, A., Rizzi, N.L., 2016. Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models. ZAMP - Zeitschrift für angewandte Mathematik und Physik 67 (4), 85. doi:10.10 07/s0 0 033- 016- 0681- 8. Vescovo, D.D., Giorgio, I., 2014. Dynamic problems for metamaterials: Review of existing models and ideas for further research. International Journal of Engineering Science 80, 153–172. doi:10.1016/j.ijengsci.2014.02.022.