The symmetric buckling mode in laminated elastoplastic micro-structures under plane strain

The present work considers lamellar (micro) structures of thin, elastic lamellae embedded in a yielding matrix as a stability problem in the context of the theory of stability and uniqueness of path-dependent systems. The volume ratio of the stiff lamellae to the relatively soft matrix is assumed low enough to initiate a symmetric buckling mode, which is investigated by analytical and numerical means. Using a highly abstracted, incompatible model, a first approach is made, and the principal features of the problem are highlighted. Assuming plane strain deformation, an analytic expression for the bifurcation load of a refined, compatible model is derived for the special case of ideal plasticity and verified by numerical results. The effect of lamella spacing and matrix hardening on the bifurcation load is studied by a finite element unit cell model. Some of the findings for the ideal plastic matrix are shown to also apply for a mildly hardening matrix material. Furthermore, the postbuckling behaviour and the limit load are investigated by simulating a bulk lamella array.


Introduction
Failure modes of composites are very different for tension or compression, and effects related to (in)stability play an important role in the latter case. For relatively high volume ratios of the reinforcements, i.e. stiff fibres or lamellae, the unsymmetric kink mode is generally regarded as the prevalent failure mode. Under simplified conditions, i.e. inextensible reinforcements without bending stiffness, several different estimates for the critical effective stress were given in the literature. The mechanism leading to kinking is either assumed to be elastic shear buckling [1], matrix yielding due to the shear stress state resulting from misaligned reinforcements [2], or a synthesis of both [3,4].
The short-wavelength symmetric (micro) buckling mode is usually only observed at very low reinforcement volume fractions in purely elastic composites [5, p. 93]. Plasticity in the matrix, however, amplifies the stiffness mismatch and seems to promote the symmetric mode. Plastic micro-buckling has been observed in detailed numerical simulations [6] for layered titanium aluminides as an intermediate mode that is overshadowed by the eventual kink failure that it precedes.
It is the foremost concern of the present article to investigate this elastoplastic symmetric micro-buckling mode and how it relates to the kink mode that is usually observed. The presence of inelasticity substantially alters the stability behaviour and requires a treatment of the problem by the stability theory of path-dependent systems. This aspect of the problem is particularly emphasized in the following.

Preparations
Unless explicitly stated otherwise, deformation processes are considered to occur in a quasistatic manner. In this context, time is interpreted as a load proportionality factor, and time derivatives are rates with respect to this factor. The matrix yield locus is assumed isotropic and smooth with an associated flow potential so that it can be modelled by standard J 2 -plasticity. The problem is reduced to two dimensions by assuming plane strain deformation. Out-of-plane deformations are in general possible in an actual structure however, and this constitutes a simplification. Possible implications of this assumption are discussed further down in the text. Furthermore, it is assumed that there is either ideal plasticity, or only moderate, isotropic hardening. Then, the total deformations are either completely or at least predominantly plastic, so the total strain can be considered isochoric while yielding. Coordinate directions are assigned so that the first direction is perpendicular to the unbuckled lamella plane, the second direction coincides with the direction of loading, and displacements in the third direction are suppressed by the plane strain assumption. The components of the rate of deformation tensor, d , are then related by (1) in terms of displacement rate components, v i , and current coordinates, x j , On the fundamental, i.e. unbuckled, deformation path, the direction of homogeneous loading falls within the lamella plane, and the principal stress axes coincide with the adopted coordinate axes. The standard J 2 -yield function, (2), is used with the stress deviator,σ , and the radius of the yield surface, κ, It has been assumed that elastic deformations are negligible while yielding. Therefore, the plastic strain rate in the third direction must vanish to comply with the plane strain assumption, i.e. ∂F /∂σ 33 = 0. With the requirement that the stress state causes yielding, i.e. F = 0, this allows the elimination of two principal stress components. Evaluation of the stress deviator and taking into account that σ 22 is compressive reveals that the remaining principal stress component cancels out and the deviator is given in terms of κ alone, regardless of the magnitude of the actual stress components, cf. (3), This further implies a constant expression for the Eulerian plastic material tangent tensor, L p , given in (4). There, the operator ⊗ stands for the dyadic product, and the colon operator indicates contraction over the last two indices. The symbol E represents the standard linear elastic material tangent tensor, and θ is the isotropic hardening parameter related to the tangent modulus of a uniaxial test, E t , as stated below, Evaluation of (4) with (3) yields the components of the Eulerian plastic material tangent tensor in terms of the Lamé parameters, λ and μ, and the hardening parameter. All nonzero components are listed in (5), and the abbreviations L 1 and L 2 are defined for reference further down. For the special case of ideal plasticity, there is θ = 0 and L 1 = L 2 , The actual material tangent tensor in the matrix depends on whether the velocity solution triggers plastic yielding or elastic unloading. For this reason, the Eulerian constitutive relation (6) between the objective Zaremba-Jaumann rate of Kirchhoff stress, • τ , and the rate of deformation tensor is incrementally nonlinear. Calculation of the plastic consistency factor for the conditions discussed above shows that it is proportional to d 11 and, therefore, permits the dependence on the velocity field, v, of the matrix material tangent tensor, L mat , to be expressed by (7) where E mat stands for the standard linear elastic material tangent in the matrix. The lamellae are modelled by beam theory and assumed to remain elastic at all times so that the Eulerian material tangent tensor of the lamella consists of only one nonzero entry E lam The mixed Eulerian-Lagrangian formulation (8) in terms of the material time derivatives of the first Piola-Kirchhoff stress tensor, . S , and the deformation gradient, .
F , is preferential to (6), and the appropriate mixed material tangent tensor, C , is obtained from the transformation (9). The reference configuration is set to coincide with the current, prestressed configuration so that Eulerian and Lagrangian coordinates are the same and the deformation gradient (but not its rate) is equal to identity. The initial stresses in the reference configuration affects the mixed material tangent tensor via the bracketed term in (9) where δ denotes the Kronecker symbol [7,8], Only discretized problems of a finite degree of freedom n are considered. Using the notation from [9], the velocity field is approximated by a field, v, assembled from shape functions, φ j , and the rates of the generalized displacements q j (GDs), i.e. . q j , in (10). Some GDs might be controlled, and it shall be required that the numbering is in such a manner that the GDs 1 to m are free, while the GDs m + 1 to n are controlled. Then, a perturbation of the discretized velocity compatible with the BCs at constant load shall be labelled by the symbol w and is given by (11), Because the reference configuration has been identified with the current configuration, the velocity gradient and the material time derivative of the deformation gradient coincide. The gradient of a field f is abbreviated as ∇( f ). Inserting the discretized velocity fields into the rate form of the virtual work principle (12) and cancelling the arbitrary perturbation δ . q i yield the rate of the internal forces .
Q i (v) on the left-hand side of (13). There, the argument v indicates the dependence on the actual velocity field because of the nonlinearity in (6). The right-hand side of (12) represents the virtual work of the external force rates due to changes in the tractions, T , which are represented in the following by concentrated forces . P i , and body forces, B, which are absent here. The external loading is assumed to be independent of the body's deformation. Equation (12) is equivalent to the condition of continued equilibrium (14), . .
The resulting tangent stiffness matrix of the system depends on v on account of the incremental nonlinearity of the material tangent tensor, (15),

Stability of path-dependent systems
Hill [10,11] showed that the key aspect of path-dependent stability problems stems from the fact that the concepts of stability and uniqueness do not coincide as they do for path-independent systems. In [9,12], an extended theory of inelastic bifurcation problems for materials that allow a rate potential was introduced on the basis of an energy interpretation. For the associative plasticity model under consideration here, this theory is applicable. The energy functional, E , is defined as the sum of the energy contained in the deformable body and the potential of the loading device. External loading is thought to be controlled independently of the body's deformation. The energy functional is related to the velocity functional appearing in Hill's theory and differs only by a factor two and additive constants [13]. It allows to assess both stability of an equilibrium state and stability of a deformation path.

Stability of an equilibrium state
For an equilibrium state under constant loading, the first time derivative of the energy functional, .
E (w), vanishes. The additional energy required to move the system from the equilibrium state in an arbitrary, but kinematically admissible, direction, at constant loading, is then given by .. (16). It is argued in [9,12] that if .. E (w) is strictly positive for all w, then spontaneous departure from an equilibrium state along a direct path is prevented by an energy barrier, and the equilibrium state is stable. If (16) does not hold for at least one nonzero w, the equilibrium is unstable, .. ..

Stability of the loading process
Uniqueness of the solution to the rate problem and stability of equilibrium are separate concepts for pathdependent systems [10,11]. A deformation path may be unstable and thus not physically realizable, even though it entirely consists of stable equilibrium states. Moreover, bifurcation points are not isolated, but form a continuous interval. The classical example to demonstrate this behaviour is known as Shanley's column [14][15][16]. Criteria to detect nonuniqueness and reject unphysical solutions are given in [8,9,12] for materials that allow for a rate potential, and meet certain constitutive inequalities (relative convexity property). For the case of associated J 2 -plasticity, those requirements are met, and the definition of a stable loading path given in the references applies: in a stable process of quasistatic deformation, the increment of E calculated with accuracy to second-order terms is minimized within the class of all kinematically admissible deformation increments. For all velocity solutions v, the first derivative .
E is the same, and the minimum is determined by the second-order terms of the energy functional ..
Hence, a velocity solution v 0 following a stable loading path is characterized by (18), This leads to a condition for the first bifurcation point. Let v 0 denote a velocity solution following the fundamental deformation path in the direction of further loading and K (v 0 ) be the associated tangent stiffness matrix. Then, positive definiteness of K (v 0 ) is necessary and sufficient for uniqueness of the velocity solution, provided the equilibrium is stable according to (17) [9]. When K (v 0 ) turns singular, bifurcation occurs, and (18) provides a criterion to choose which path to follow in the postbifurcation regime. Beyond this first bifurcation, there usually exists an interval of continuously arranged bifurcation points on the fundamental deformation path where the stability condition for equilibrium (17) still holds, but the stability condition for the loading path (18) does not. The first bifurcation point can only be exceeded on the fundamental path by temporarily imposing additional constraints. The bifurcation point interval terminates when the undeflected equilibrium itself becomes unstable and (17) fails [9]. Departure from the fundamental deformation path initiated by failure of (17) takes place suddenly at constant load, but deflection due to failure of (18) with (17) still holding occurs only gradually with increasing load.  Fig. 1 The introductory model. a Undeformed, b buckled

Introductory model
In order to relate the problem under consideration to the theory of inelastic buckling reviewed in the previous Section, a first approach is made by using a highly simplified, introductory model. For this purpose, suppose an infinite periodic formation of lamellae and matrix layers under plane strain. The lamellae are thought to be connected to the matrix only at equidistantly distributed individual pinning points, rather than by a continuous bedding. It is further assumed that the buckling mode of these lamellae is such that there is symmetry with respect to the midplane between two lamellae and the planes orthogonal to the lamellae passing through the connection points. A justification for this symmetry assumption is given in Sect. 5. In this regular assembly, the spacing of the connection points defines the half-wavelength of the lamella's buckling mode.
A rectangular cut-out extending between two neighbouring connection points and matrix-midplanes tiles the entire cross section by suitable mirror transformations, and one such cut-out is, therefore, representative for the entire formation. Two cut-outs are depicted in Fig. 1 side by side, one of them outlined in Fig. 1a. The lamella is assumed to be thin enough to be handled by beam theory modified for plane strain, while for the matrix deformation a crude kinematic assumption is made: the matrix in the cut-out is thought to be divided into four subdomains that remain rectangular during deformation and slide without friction relative to each other. The model is loaded by the external force P, which is in balance with the internal forces, resulting from matrix and lamella stresses. When deflection of the lamella occurs, a void forms at the centre of the model, cf. Fig. 1b. The lacking accuracy of this introductory model in representing the actual situation makes it unsuitable for quantitative predictions, but it allows a description simple enough to be analysed in the context of the theory and still captures the lamella-matrix interaction at least in principle. The so-defined model has only four GDs related to the displacement fields in the lamella and the matrix by the shape functions given in (19),

Internal force rates and stiffness
By adapting (13) to the present model, the quasistatic internal force rates per unit depth are given in (20). Note the nonlinear dependence on the velocity field, and, eventually, its gradient, via (7), The rate of deformation tensor is constant in each of the four matrix subdomains and depends on the current coordinates in a fixed frame only via the signum function. The abbreviationsd andd are introduced for convenience, It has been found in (7) that yielding takes place when d 11 > 0 and (21) assigns each subdomain a yielding or unloading state, depending on the value of the GD rates. Only the four cases listed in Table 1 need to be considered, and (20) is resolved into the expression (22), .
The superscript L ∈ {ll, ul, lu, uu} indicates in which subdomains yielding and unloading are active, and C L i and D L i are given by (23), In several terms in (22), a direct subtraction of stress components from the material constants C L i or the lamella's Young's modulus is present in the expression for the internal forces. Since C L i is of the same order of magnitude as the matrix Young's modulus, it seems pertinent to neglect the stress components in these expressions. For the matrix, this has the same effect as completely neglecting its geometric stiffness and (20). Then the tangent stiffness matrix is given by (25), where the superscript L is representative of the dependence on v.
The tangent stiffness matrix depends on the loading via σ lam 22 in K L 11 . The usual sign convention, i.e. a negative sign for compressive stress, applies.

Bifurcation and stability
With the stiffness for each load case at hand, it can be shown that the bifurcation and stability behaviour of the introductory model are very similar to the classic example of Shanley's column.
Starting from an undeformed state, the loading path first traverses the elastic regime. There, the tangent stiffness matrix is given by K uu and does not depend on the the direction of the velocity field. The lowest eigenvalue is initially positive but decreases with progressive loading. If it turns zero in the elastic regime, the criteria for stability of equilibrium, (17), and for the stability of the loading path, (18), fail simultaneously. In this case, elastic buckling occurs in the usual manner, and deflection can occur at constant load because the system is in indifferent equilibrium, according to second-order theory. Elastic buckling shall be excluded in the following.
The symbol σ lam y is used in the following to denote the absolute value of the lamella stress when the matrix yields for the first time. Note that in the model the lamella remains always elastic. Upon yielding the incremental, i.e. tangential, stiffness changes discontinuously to K ll , and this requires some consideration. The discontinuity might render the tangent stiffness matrix indefinite without the lowest eigenvalue turning zero first. This is a somewhat special case that is discussed further down, but the case when the lowest eigenvalue of K ll is still positive upon initial yielding shall be considered before. In this event, the instant, when (18) fails for the first time, is distinct from the instant of initial yielding, and the system continues on the fundamental deformation path in loading direction with all subdomains in a yielding state.
Positive semidefiniteness of the incremental stiffness is equivalent to first failure of the loading path stability criterion (18). Because all eigenvalues have been previously positive, this load state is characterized by (26), From the criterion (26), the lamella stress at the first bifurcation point, σ * , follows as stated in (27). This stress marks the first of an interval of bifurcation points, cf. the shaded area in Fig. 2. The fundamental equilibrium path has no physical relevance beyond this point, unless additional constraints are present. The right-hand side of (27) consists of two terms with clearly recognizable interpretations. The first term represents the Euler buckling load of the laterally unsupported lamella, and the second term expresses the effect of the lateral support exerted by the matrix on the lamella. In case of ideal plastic matrix behaviour, the second term vanishes, and the matrix does not affect the buckling stress. The bifurcation stress is assigned a positive value by definition, even though it is compressive, The surfaces representing the energy functional at different lamella stresses are shown in Fig. 3a-c as contour plots. See also [13] for plots of the energy functional of Shanley's column. For the purpose of generating these plots, the GD rates . q 2 and . q 3 were chosen so that the energy functional is minimal with respect to them. Numerical values are calculated for a lamella spacing of 2 mm, a fixed half-wavelength of 1 mm, and ideal plasticity in the matrix. Other parameters are taken from Table 2 where a lam stands for the lamella thickness.
The plot for a load state with the lamella loaded to 3/4 of σ * and . P = 10 3 N/s is shown in Fig. 3a. At this state, a single minimum of E .. exists, which is obtained by the velocity solution following the fundamental deformation path, v 0 . Bifurcation is, therefore, excluded. The situation changes when the lamella stress is increased to σ * , and it can be seen in Fig. 3b that the condition (26) manifests itself in the plots as a degeneration of the ellipsoidal isolines to parallels in the range of the GD rates that activate yielding in all subdomains, i.e., the region characterized by 0 <d <d, which is shaded in the plots. There is no longer a single minimum point for the functional, but the functional is minimized on an entire line segment spanning the yielding-everywhere region, and this admits departure from the fundamental deformation path. An ever so slight load increase results in three stationary points, two minima at both endpoints of the line segment and a saddle point at . Fig. 3b, the saddle point is indicated by v 0 and the right minimum point by v * . According to the requirements for a stable loading path discussed in Sect. 3, deflection must occur if . P > 0. With further increasing load in the postbuckling regime, unloading in two subdomains takes place as the solution now falls in the region 0 <d < |d|, cf. Table 1.
Stability of equilibrium, on the other hand, is determined by (17), and the energy functional is to be evaluated at constant loading. The velocity dependence of the internal force rate is indicated by the superscript L in (22), and the superscript has to be paired with a given direction of w as specified by (21) and Table 1. Consequently, the velocity dependence of the internal force rates, i.e. , the activation of different constitutive branches in different directions, results in a positive definite expression for the left-hand side of (17), and, thus, stability of equilibrium at constant load. As an illustration, refer to Fig. 3d. Even though there is zero curvature with respect to . q 1 in the yielding-everywhere region (shaded in the plot), departure in direction . q 1 is prevented by the curvature of the adjacent surface patch for L = ul. At a certain higher stress σ * * (see Fig. 2), the energy barrier in the lu region will vanish in a particular direction, and, thus, equilibrium becomes unstable, if the load P is controlled. This state is characterized by det K lu = 0, and lamella deflection can occur at constant load.
Even though the first bifurcation point can generally not be exceeded on the fundamental bifurcation path by a continuous loading process, additional constraints can keep the lamella undeflected beyond σ * . Typically, this occurs when σ * is less than the stress required to cause yielding in the matrix. While σ * < −σ lam 22 < σ lam y < σ * * , the entire assembly is still in an elastic state, and plastic stability considerations do not apply. When σ lam y is reached, condition (18) immediately fails, and deflection occurs at an elevated stress, but still with increasing load.
In the actual lamella-matrix compound, no constraint on the buckling wavelength is apparent. In the simplified model under consideration in this Section, the model height a 2 was defined to be equal to the halfwavelength of the buckling mode and must hence be considered as unknown. It can be obtained from (27) by specifying that only the half-wavelength a 2 , that minimizes σ * , can occur, cf. (28), It is noteworthy that a 2 tends to infinity for ideal plasticity (θ → 0 in (24)). This is a consequence of the kinematical assumptions, and a refined model for ideal plastic matrix behaviour, that is discussed in Sect. 5, reveals a constraint for a 2 arising from the equations governing the matrix deformation.

FEM Model
The analytical findings were verified by numerical simulations, using the commercial finite element method (FEM) software Abaqus [17]. The FEM model is very simple and consists of a single, first-order interpolation, plane strain element for each matrix subdomain, and several beam elements for the lamella. The sliding conditions are implemented by linear constraint equations for nodal displacements. Horizontal expansion is free, and σ 11 is zero on the fundamental deformation path. Due to the incremental nonlinearity (7), the bifurcation load cannot be obtained via a linear perturbation step and has to be derived from the equilibrium path of a quasistatic analysis. Stable postbuckling behaviour is expected, but in order to extend the calculations beyond a limit load that might be reached in the postbuckling regime, the modified Riks algorithm is employed. The standard solution algorithm of the software does not have any measures in place to ensure convergence to the solution that minimizes E at the bifurcation point. As a remedy, a system with a tiny geometric imperfection of 1×10 −5 mm lateral offset of the topmost lamella node is simulated. The equilibrium path of the imperfect system closely approximates the path of the perfect system, but does not cross the first bifurcation point. Figure 4 shows load-displacement diagrams for undeformed lamella spacings of 2 mm and 4 mm, and an initial model height of 1 mm. The tangent hardening modulus is set to zero in both cases, and other parameters are given in Table 2, so the bifurcation loads calculated from (27) coincide (dash-dotted line). The start of the lateral deflection closely matches the expected value with a certain deviation due to shortening of the model in the prebuckling regime, cf. Fig. 4a, c. The dashed curve in the Figures stands for the ratio, r , of currently yielding elements, i.e. the 'yield fraction'. It can be seen in Fig. 4b, d that half of the subdomains are unloaded after bifurcation. Although the axial stiffness decreases only gradually after bifurcation, a limit load is reached eventually. The mechanism determining the limit load appears to be very different for the narrow and the wide lamella spacing. In the first case, the postbuckling stiffness hardly decreases until the unloading has proceeded to such an extent that yielding occurs in opposite direction. This sudden change in stiffness then triggers unstable collapse with a snap-back behaviour in the quasistatic case. For the wide lamella spacing, the axial stiffness softens faster, and a snap-through limit point is reached before the stress state in the unloaded matrix subdomains reaches the yield surface on the opposite side. It is shown in Sect. 5.2 further down that this behaviour is present for a compatible model as well.

Compatible unit cell model for ideal plasticity
Biot formulated the governing equations for the incremental deformation of a prestressed, possibly orthotropic, incompressible, elastic material without body forces, (29.1) to (29.4), in terms of a fourth-order partial differential equation (PDE) [18]: 12 (29.1) The incompressibility condition is identically fulfilled by defining a stream function, ϕ, in the manner of (30.1). Likewise, the equilibrium can be fulfilled by a stress function, ψ, (30.2), Combining the remaining two equations and eliminating the stress function yield the fourth-order PDE (31) [18, p.193]. Elimination of the stream function leads to the same equation, but with ψ replaced by ϕ [19]. This equation was used in [20] to calculate the critical value of p for internal buckling of homogenized, elastic laminates, An approach different from Biot's is pursued here, and lamella and matrix are kept as separate domains, but (31) is used to describe the displacement field in the matrix. For this purpose, it shall be presumed that the lamella stress at matrix yielding, σ lam y , is less then σ * , i.e. the lamella stress at the first bifurcation point. This means that elastoplastic bifurcation is not suppressed by a very high yield stress. Then, as has been argued in Sect. 3, the first bifurcation point can be expected with the entire matrix in a yielding state. Limiting further considerations to ideal plastic materials, it is found from (6), the left equation in (29.1), and the identification of the reference configuration with the current configuration that the material constant μ * is equal to zero, (32), Further assuming that the geometric stiffness of the problem is dominated by rotations of the thin lamella, as it has been found to be the case for the incompatible model, and neglecting the geometrical stiffness in the matrix by setting p = 0 leads to the simplified PDE (33), If the matrix consisted of an isotropic, elastic material, then μ would equal μ * , and the biharmonic equation ∇ ∇ ϕ = 0 would have been obtained. In (33), however, the middle term has a negative sign due to the ideal plastic behaviour. The equation can be solved by a function ϕ( In the present case, two deformation modes exist in the matrix: homogeneous shortening and the strain fluctuations resulting from lamella deflection. In analogy to the incompatible model, these modes are represented by their amplitudes, the GDs q 1 and q 4 , and the respective shape functions, φ mat 1 and φ mat 4 , which are to be determined from the PDE. The lamella deflection is assumed sinusoidal, and the function ϕ is thus identified with a trigonometrical function. The solution given in (34) fits the deformation sketched in Fig. 5 with b 1 and b 2 as real constants. The last term is added to account for homogeneous strain, and it is clearly also a solution of (33), To determine the constants b i and c i , the assumption that the lamellae are very thin is invoked, and, therefore, vertical displacements at the interface due to rotation of the lamellae's cross sections are negligible. This suppresses modes that might be present for laminates where the fractions of both phases are similar. Such a mode has been observed in FEM simulations for a lam /(2a 1 ) close to 1. The conditions to match the lamella displacement at the left boundary are given in (35) leading to (36).
The sign of b 2 is undetermined and does not affect the velocity field. The matrix velocity field is: With all constants determined, the BC on the right-hand side interface still needs to be accounted for. Beforehand the model as shown in Fig. 5 shall be given an interpretation as a (minimum) periodic unit cell (UC). For the assumed case of evenly spaced lamellae, the micro-structured material constitutes an inhomogeneous periodic medium. It can be efficiently described via a cut-out representing a whole period of the repeating arrangement, the UC [21]. The shape of the UC is taken as rectangular for the considerations here. UCs are often utilized to calculate homogenized or 'smeared' properties for materials that are heterogeneous on a length scale small enough to meet the following homogenization conditions. The macroscopic properties can be derived from the structural properties at a smaller length scale [22]. This approach gives meaningful results only if the length scale at the micro-level and the length scale of the macroscale are sufficiently different, so that the gradient of the macro scale field variables can be neglected on the micro-scale. On the micro-level, i.e. in the UC, the local strain at a particular point x is given by (38), whered stands for the macroscopic strain andd for the local fluctuation of strain. The fluctuations do not contribute to the macroscopically apparent strain and cancel out when averaging the strain field over the UC, Likewise, the velocity field can be decomposed in mean and fluctuating parts, cf. (39.1) where the mean velocity is expressed in terms of the local coordinates and the mean strain. The velocity fluctuations are periodic in two directions, and, therefore, take the same value when an integer multiple of two linear independent vectors, r a and r b , is added to the argument x. These vectors are called the periodicity vectors [23], 2) t = n a r a + n b r b n a , n b ∈ N. (39.3) Returning to the case at hand, the decomposition of v mat is immediately apparent from (37) withṽ mat = .
q 4 φ mat 4 . The periodicity vectors are given in (40). There, an allowance for a phase shift angle, β, between the mode of two neighbouring lamellae is made, rather than prejudicing the deformation mode as symmetric by setting β = π/4. At this point, it must be mentioned that the plane strain assumption introduced in Sect. 2 suppresses both the macroscopic and the fluctuating out-of-plane strain. This is a simplification, and in a more realistic model one would assume (generalized) plane strain for macroscopic strain only, while still allowing for out-of-plane fluctuations. Such an extended model would require a 3-dimensional formulation of the PDE (31) and is not considered here, The fluctuating part of the velocity field in the UC adjacent to the right of the considered UC is the same for two corresponding points x and x + r a , (41). If x is located on the left boundary of the UC, then this relation ensures compatibility on the right boundary, Explicit calculation of (41) for the periodicity vector specified in (40.1) and application of trigonometric identities show that the conditions stated in (42) must be met for compatibility, This ultimately results in the requirement that the buckling wavelength, a 2 , must be an integer multiple of the lamella spacing, 2a 1 , (43). Obviously, the lowest buckling load corresponds to the longest possible wavelength, and, therefore, k takes the value of 1, With the shape functions at hand, the internal force rates are calculated by (13) for the yielding everywhere case which is indicated by the superscript l. The matrix stress is capped by the yield stress, and consistently with setting p = 0 in (31), the matrix geometric stiffness is neglected by replacing C mat (v) with L p . Although the expression L 1 − L 2 evaluates to zero in the ideal plastic case, it is retained here to allow for an approximation of the mildly hardening case that is discussed in Sect. 5.1 further down, .
Bifurcation occurs when det K (v 0 ) = 0, which reduces here to ∂ . Q 1 /∂ . q 1 = 0. For ideal plasticity, the bifurcation stress coincides with the critical stress of a pin-ended Euler beam with a length equal to the lamella spacing. The second term in (45) vanishes, and the matrix does not provide any direct support against lateral deflection. The matrix geometry does, however, affect the buckling load via a 1 ,

FEM UC model
To verify the analytical findings regarding the bifurcation load, finite element method (FEM) simulations were performed. For this purpose, the bifurcation load of the numerical model is defined as the load at which unloading occurs for the first time anywhere in the matrix. The UC used for the numerical calculations differs from the minimum UC in Fig. 5. To enforce periodicity, equations for nodal displacements in analogy to the rate equations (39) are required, with x referring to current coordinates. The finite element program, however, uses a total Lagrangian formulation, and equations relating nodal displacements are to be specified with respect to their undeformed position. Setting t equal to r a in the undeformed configuration, subjects the equations to the homogeneous, but often substantial, prebuckling deformation distorting t and thereby affecting the phase angle. To nullify this effect, the numerical UC is enlarged horizontally to span 4a 1 , and t is set to 2r a − r b . In a regular mesh, this allows to specify the periodicity equations between nodes at the same vertical position, and, therefore, cancels the effect of the homogeneous deformation on the phase angle. In vertical direction, the UC needs to be enlarged as well. Because the number of full periods within the UC must be integer, and the relation a 2 = 2a 1 is to be verified, the UC should be as high as possible to avoid undue interference of UC height and wavelength. To ensure this, the numeric UCs generally have a height of 40 times the lamella spacing. Again, these proportions refer to the undeformed configuration. Furthermore, the origin of the FEM UC is translated by −a 1 in horizontal direction relative to the origin in Fig. 5. As in Sect. 4.3, the mean horizontal expansion is unconstrained in the simulations, and a geometric imperfection is introduced by a random variation of the nodal positions with an amplitude of 10 −4 mm. Geometric and material parameters are reused from Table 2.
Numerically determined bifurcation stress values are compared to values predicted by (45) in Fig. 6. A series of simulations for different lamella spacings, Fig. 6a, shows good agreement in the range 0.5 mm ≤ a 1 ≤ 1.5 mm and a noticeable discrepancy for 1.5 mm < a 1 ≤ 2.0 mm. In this range, the initiation of bifurcation in the FEM model appears to be sensible to solver parameters like load increment size. For the largest lamella spacings, the predicted bifurcation stress approaches the lamella stress at matrix yielding, σ lam y . The bifurcation stress is inversely proportional to the square of the lamella spacing. Therefore, the symmetric mode is likely superseded by other modes for small values of a 1 . For the intermediate value of a 1 = 1 mm, the horizontal component of the normal strain at the increment of buckling is depicted in Fig. 6b. Lamella beam elements are not visible in the Figure, but their approximate horizontal position is indicated by the vertical double lines added to the plot. The strain field is caused by prebuckling deformation of the imperfect system and is affine to the strain rate field triggered by the buckling mode. Comparison to the strain rate field derived from (37) shown in Fig. 6c shows good agreement.
The velocity field (37) was derived under the assumption of ideal plastic matrix behaviour, but the internal load rates, (44), were still calculated for the general case. Under the assumption that matrix hardening has little effect on the velocity distribution, this allows to use (45) as an estimate for the general case as well. To investigate the suitability of this assumption and to assess the sensibility of the bifurcation stress to matrix hardening, the predicted values are compared to FEM results in Fig. 7a. The approximation appears to hold up well for the rather small range of values for the tangent moduli considered. Even so, the strain field qualitatively changes at higher hardening moduli, as can be seen in Fig. 7b, invalidating the assumption for larger hardening moduli. It becomes also apparent that the sensitivity of the bifurcation stress to matrix hardening is very pronounced.

FEM array model
To further corroborate the previous results, models comprising an entire array of ten lamellae and free from any periodicity BCs were simulated. It is shown further down that in these models the strain pattern of the symmetric mode arises naturally and interference by periodicity constraints can be ruled out. Furthermore, the array models allow an extension of the simulations beyond the limit point into the regime of unstable collapse. The UC models are not suitable for this purpose, because the collapse mode cannot be expected to follow the same periodicity requirements as the bifurcation mode. A minor drawback of the models is, however, that they introduce a slight free edge effect.
In the array models, the horizontal displacements of all but one node are unconstrained, and vertical displacement is controlled at the top and the bottom of the model to prescribe the mean engineering strain component,ε 22 . To simulate the collapse, a dynamic analysis was performed, but the effect of inertial forces is kept to a minimum in the stable regime. To obtain convergence in the collapse regime, an artificially increased  mass density was necessary. The Riks method fails to converge after the limit point, and following the snapback segment of the equilibrium path would lead to unphysical results on account of the path dependence of the material. The undeformed model height is 20-times the lamella spacing, and the material is ideally plastic in all simulations. Here the lamellae are modelled by continuum elements. They are not shown in the plots, but their location becomes apparent from the white gaps between the matrix sections. Again, an imperfect model is simulated with the same random initial imperfections as in Sect. 5.1. Unspecified parameters are taken from Table 2.
The normal strain component in horizontal direction at an increment before buckling is shown in Fig. 8. The figure is cropped, and only the top section of the model is shown. The strain pattern is essentially the same that has been observed in the UC-model Fig. 6b. In Fig. 9, the total load, i.e. the reaction force induced by the controlled displacement, is graphed for a narrow and a wide lamella spacing. The general behaviour is analog to the results in Fig. 4. For the narrow lamella spacing, re-yielding of elastically unloaded regions takes place almost simultaneously with reaching the limit point, while for the wide spacing the ratio of yielding integration points is strictly decreasing. During the collapse, the plastic deformation localizes in a band at an angle of roughly 45 • to the lamella, Fig. 10. The maximal postbuckling amplitude remains relatively small, and a typical deflection outside the localization zone is marked by the arrow in Fig. 10b. This circumstance and the initially marginal effect of bifurcation on the axial stiffness might make it difficult to detect symmetric buckling in experiments.

Conclusions
Elastoplastic buckling in materials with lamellar structures has been discussed in the context of the theory of stability of path-dependent systems. It has been shown to feature some of the essential properties predicted by the theory as, for instance, stable postbuckling behaviour. From the analytical considerations and the results of the simulations, it appears that the symmetric buckling mode can conceivably occur at moderate lamella stress levels for a matrix to lamella volume fraction of about 10 or larger with very low matrix hardening moduli. However, elastoplastic buckling might be difficult to detect in experiment as the apparent axial stiffness of the material is affected only gradually by the bifurcation itself, and lamella deflection amplitudes remain small. Even so, it appears to be an important mechanism for certain configurations as it eventually leads to unstable collapse in the postbuckling regime. The limit loads are typically much larger than the bifurcation load. The collapse seems to be caused by a sudden change in stiffness due to re-yielding of formerly unloaded matrix areas for narrow lamella spacings or by softening due to lateral lamella deflection for wider lamella distances. During collapse, plastic deformations localize in a band. This appears to be in accordance with the localized mode observed in [6].