M Ahmer Wadee

Stylianos Yiatros

Marios Theofanous

Comparative studies of localized buckling in sandwich struts with different core bending models

Keywords: Structural stability, bifurcation, interactive buckling, localization, shear deformation, composite structures

Analytical models with geometric nonlinearities accounting for interactions between local and global instability modes leading to localized buckling in sandwich struts are formulated. For the core material response, two increasingly sophisticated bending models are compared against each other: Timoshenko beam theory (TBT) and Reddy-Bickford beam theory (RBT). Numerical solutions of the analytical models are validated with the commercial finite element code Abaqus. It is found that there is a small but significant difference in the critical load between the two models and that the previously obtained solution slightly underestimates the linear buckling strength. More importantly, it is found that the RBT model predicts the onset of interactive buckling before the TBT model and, according to the results from the finite element study, matches the actual behaviour of a strut in both its initial and advanced post-buckling states with excellent correlation.

Introduction

Sandwich construction, comprising two stiff face plates separated by a softer core material, is popular as a provider of structural strength combined with

A c c e p t e d m a n u s c r i p t

weight efficiency and is used extensively in astronautic [START_REF] Kodiyalam | Composite sandwich structure optimization with application to satellite components[END_REF], aeronautic [START_REF] Bannink | Advanced design composite material aircraft study[END_REF][START_REF] Hsiao | Core crush problem in manufacturing of composite sandwich structures: mechanisms and solutions[END_REF] and marine applications [START_REF] Knox | Adhesively bonded steel corrugated core sandwich construction for marine applications[END_REF]. However, precisely because they are both specialized and efficient, the responses of sandwich struts are liable to exhibit complicated collapse mechanisms [START_REF] Hunt | Interactive buckling in sandwich structures[END_REF][6][START_REF] Sokolinsky | Nonlinear behavior of sandwich panels with a transversely flexible core[END_REF][START_REF] Fleck | End compression of sandwich columns[END_REF]. It is well known from classical work [START_REF] Allen | Analysis and design of structural sandwich panels[END_REF] that compressed sandwich panels, or more specifically struts, sometimes fail by a combination of overall (Euler-type) buckling and local buckling (wrinkling) of the face plates. Figure 1 shows a sequence of test photographs [START_REF] Wadee | Experimental evaluation of interactive buckle localization in compression sandwich panels[END_REF] along with a characteristic equilibrium diagram of compression sandwich panels representing the different phases of the loading history. The phases in sequence are: (1) pre-buckling followed by (2) Euler buckling and then proceeding to (3) interactive buckling where the structure becomes unstable and localization is clearly observed. This type of structural response has been previously modelled using a combination of nonlinear structural stability theory and a Timoshenko beam approach; the appearance of shearing strains within the core material being vital in introducing a nonlinear interaction between the overall buckling wavelength scale and the local buckling "strut on an elastic foundation" wavelength scale [START_REF] Hunt | Comparative lagrangian formulations for localized buckling[END_REF][START_REF] Hunt | Localization and mode interaction in sandwich structures[END_REF]. More advanced work on the sensitivity due to imperfections [START_REF] Wadee | Effects of periodic and localized imperfections on struts on nonlinear foundations and compression sandwich panels[END_REF], the possibility of face-core delamination [START_REF] Wadee | Delamination from localized instabilities in compression sandwich panels[END_REF][START_REF] Wadee | Localized buckling in sandwich struts with pre-existing delaminations and geometrical imperfections[END_REF] and the buckling of panels with differing face plate thicknesses [START_REF] Wadee | Asymmetric secondary buckling in monosymmetric sandwich struts[END_REF] has also been presented. In spite of significant progress, these earlier works have focused only on the more simplified case of plane sections remaining plane within the core material once buckling occurs, although not necessarily normal to the neutral axis of bending. This, of course, introduces shear strains into the model [START_REF] Attard | Sandwich column buckling -a hyperelastic formulation[END_REF], but for the more practical cross-section geometries, principally with deeper and rela-

A c c e p t e d m a n u s c r i p t

tively softer cores, the local shear strains-particularly in the neighbourhood of the face-core interface-may distort the core nonlinearly such that the assumptions used previously may need updating. Higher order bending models have been developed [START_REF] Frostig | A high order theory for the bending of sandwich beams with flexible core[END_REF] and been applied in the literature to account for stress concentrations in the vicinity of lateral point loads [START_REF] Frostig | Localized load effects in high-order bending of sandwich panels with transversely flexible core[END_REF][START_REF] Frostig | Localized effects in the nonlinear behavior of sandwich panels with a transversely flexible core[END_REF] and thermal buckling studies [START_REF] Frostig | Thermal buckling and postbuckling of sandwich panels with a transversely flexible core[END_REF]. In the current work, the core is also represented with a higher-order bending model, so-called Reddy-Bickford Theory (RBT). This distributes the distorting displacement of the previously plane section as a cubic polynomial function [START_REF] Wang | Shear deformable beams and plates: Relationships with classical solutions[END_REF], leading to a quadratic shear strain distribution in the cross-section that diminishes to zero at the top and bottom surface which maintains compatibility with the adjacent face plates. The aim of the current work is to evaluate the mechanical accuracy of the elastic behaviour from models formulated with the RBT approach and with Timoshenko Beam Theory (TBT) approach, the latter being used in the earlier work outlined above.

The current paper begins with the development of the analytical model for the sandwich strut with RBT bending in the core material. A system of nonlinear ordinary differential equations and integral constraints is derived from minimizing the total potential energy using variational principles. This system of equations is solved within the numerical continuation software package Auto97 [START_REF] Doedel | AUTO97: Continuation and bifurcation software for ordinary differential equations[END_REF]. The results from the analytical models developed with RBT currently and TBT previously are compared against each other and validated in conjunction with a purely numerical model formulated within the commercial nonlinear finite element software Abaqus [START_REF] Abaqus | ABAQUS/Standard: User's manual version 6[END_REF]; the latter being used as a benchmark. The practically significant results are presented with discussions on the effects of changing various physical parameters that significantly affect the behaviour; particular emphasis is placed on the relative validity of each analytical model and the parameter ranges where the principal differences lie. Conclusions are then drawn.

Definitions

The fundamental interactive buckling model previously developed [START_REF] Hunt | Localization and mode interaction in sandwich structures[END_REF]6] adopts a Timoshenko Beam Theory (TBT) approach, which deviates from Euler-Bernoulli Theory (EBT: also known as "Engineer's Bending Theory") by relaxing the constraint of orthogonality between the deflected neutral axis of bending and the plane section. This introduces shear strains in the core material during bending that are neglected in EBT and have been shown to be vital in allowing the prediction of interactive buckling in sandwich struts that leads to localization from an analytical perspective. A higher order bending theory such as the Reddy-Bickford Theory (RBT) takes this a step further by relaxing the assumption that forces plane sections to remain plane [START_REF] Wang | Shear deformable beams and plates: Relationships with classical solutions[END_REF].

A c c e p t e d m a n u s c r i p t

This makes the theory more suitable for structural elements with deep crosssections or softer cores, allowing the axial displacements to vary nonlinearly over the depth of the section as shown in Figure 2. More importantly, a higher order theory such as RBT satisfies the vanishing shear strain condition at the top and bottom flange while having a quadratic distribution through the depth of the cross-section [START_REF] Reddy | A simple higher order theory for laminated composite plates[END_REF][START_REF] Reddy | A general non-linear third-order theory of plates moderate thickness[END_REF]. Poisson's ratio ν. For the core material, provision for orthotropic behaviour is included by allowing different Young's moduli and Poisson's ratios in the x and y directions-E x and ν x with E y and ν y respectively. The overall (Euler-type) buckling mode is represented by two generalized coordinates that comprise the amplitudes of a sway and a tilting mode. These functions are taken as trigonometric in x with dimensionless amplitudes of q s and q t respectively for the sway and tilt components, thus:

dW dx θ(x) dW dx = θ(x) dW dx θ(x) dW dx = θ(x) dW dx θ(x) dW dx = θ(x) (
W (x) = q s L sin πx L , θ(x) = q t π cos πx L , (1)
and the overall mode components are shown in Figure 4. The expressions for W and θ are given for the buckling mode of a simply-supported strut, found by solving the linear Euler strut equation [START_REF] Timoshenko | Theory of elastic stability[END_REF], which is the most common support case. Alternative support conditions can easily be implemented by changing W and θ to reflect these and the same procedure can henceforth be followed. The corresponding shear strain, assumed to be zero in EBT, is the Local-bottom: u(x) and w(x) becoming non-zero beyond the second bifurcation.

The generalized coordinate Δ gives the purely compressive strain arising from the axial load P .

difference between the slope of the deformed neutral axis due to the lateral deflection W (x) and the angle of tilt θ(x) such that:

γ xy = dW dx -θ (2)
The localized mode, w(x), is an initially unknown function that is maximum at the bottom face plate and is assumed to decrease linearly to zero at the top face. The associated local in-plane displacement u(x), which arises from the formation of w, follows the same distribution, thus:

u c (x, y) = b -2y 2b u(x), w c (x, y) = b -2y 2b w(x), (3)
where u c and w c represent the distribution of the local modes through the core and are assumed to be linear in y; this is a pragmatic choice as it had been found in preliminary studies that a cubic distribution in y, which would be similar to RBT, yielded similar results to the linear case. Hence, for simplicity the linear distribution is kept as in earlier work [START_REF] Hunt | Localization and mode interaction in sandwich structures[END_REF]. The localized mode components are shown in the final diagram of Figure 4.

It should also be noted that displacement functions accounting for localized buckling in the top face plate are not considered in the modelling. The reason is that for the panels considered presently, overall buckling is assumed to occur first since the critical load for local buckling for the face plates is usually greater. However, even in the cases where local buckling or wrinkling

A c c e p t e d m a n u s c r i p t

occurs first, previous studies [START_REF] Wadee | Effects of periodic and localized imperfections on struts on nonlinear foundations and compression sandwich panels[END_REF][START_REF] Wadee | Asymmetric secondary buckling in monosymmetric sandwich struts[END_REF] have clearly shown that the overall buckling is triggered almost instantaneously leading to interactive buckling and localization in the way described by the current model.

Reddy-Bickford Theory and Model Formulation

The formulation of the model follows the procedure outlined in [6], but it now incorporates a new expression for the in-plane tilt displacement during overall buckling, u R :

u R (x, y) = -y θ(x) + 4y 2 3b 2 ∂W ∂x -θ(x) . (4)
Note that the first component of the tilt displacement is the same as for the TBT model whilst the second component provides the nonlinear variation that releases the common "plane sections remain plane" assumption. Differentiating this expression with respect to y gives the new angle of tilt θ R that has a quadratic distribution in y and depends on both q s and q t :

θ R (x, y) = - ∂u R ∂y = θ(x) + 4y 2 b 2 ∂W ∂x -θ(x) . (5)
It is worth noting that to avoid further complication to the model the distributions of u c and w c remain linear through the core depth, as given in Equation (3).

For the nonlinear buckling problem, the total potential energy V of the strut under axial loading is formulated. It is composed of the strain energy U minus the work done P E, where P is the axial load and E is the total end-shortening at the point of load application. The strain energy is integrated over the volume of the strut and is derived from bending and membrane strains in the face plates together with axial, transverse and shear strains in the core. In the current model it is assumed that the core behaves linearly elastically upon loading; nonlinearities arising from the core constitutive behaviour under axial compression [START_REF] Hunt | Localization and mode interaction in sandwich structures[END_REF] and from face-core delamination [START_REF] Wadee | Delamination from localized instabilities in compression sandwich panels[END_REF] are left for further work.

The total potential energy V of the current system is accumulated in a hybrid formulation using continuous functions and generalized coordinates. This is minimized using variational principles [START_REF] Thompson | A general theory of elastic stability[END_REF][START_REF] Fox | An introduction to the calculus of variations[END_REF] to produce a system of equilibrium equations comprising ordinary differential and integral equations. The initial bifurcation point can be found through classical linear eigenvalue analysis. Moreover, solving the nonlinear equilibrium equations simultaneously allows the position of the secondary bifurcation to be found together with the post-buckling path and its associated localized buckling mode.

A c c e p t e d m a n u s c r i p t

3.1 Strain energy

Face plates

There are two components of strain energy stored within the face plates: from bending and membrane action. The former component arises from the local and global curvature of the face plates along the strut:

U b = EI 2 L 0 2 Ẅ 2 + ẅ2 dx. ( 6
)
where EI is the flexural rigidity of one face plate about its local minor axis, hence:

EI = Ect 3 12 (1 -ν 2 ) . (7)
Note that the localized mode w is confined only to one face and dots denote differentiation with respect to x. The membrane strain energy component arises from the axial strain in the two face plates due to squashing and subsequent buckling:

U m = D L 0 ε 2 xb + ε 2 xt dx, (8)
where D = Etc/2 with ε xt and ε xb being the axial strains in the top and bottom face plates respectively. The axial strain in the top face is given by:

ε xt = - b 2 θ - b 6 Ẅ -θ -Δ, (9)
with terms accounting for pure squash strain Δ and overall buckling, while for the bottom face plate, which is susceptible to local buckling:

ε xb = b 2 θ + b 6 Ẅ -θ -Δ + u + 1 2 ẇ2 , (10)
there are two additional von Kármán strain terms that are familiar from large deflection plate theory [START_REF] Timoshenko | Theory of elastic stability[END_REF].

Core material

The strain energy stored in the core has three sources: axial, transverse and shear strains. Firstly, both the axial and transverse strain contributions are considered. These arise from the direct strains where conditions of plane stress (σ z = τ xz = τ yz = 0) and orthotropy are assumed; the latter assumption giving the following reciprocal relationship between the Young's moduli E x and E y and the Poisson's ratios ν x and ν y :

E x ν y = E y ν x . (11)

A c c e p t e d m a n u s c r i p t

Therefore, the core energy from axial strains is:

U ca = c 2 (1 -ν x ν y ) L 0 b/2 -b/2 E x ε 2 x + E y ε 2 y + 2ν x E y ε x ε y dy dx. ( 12
)
where:

ε x = -y θ - 4y 3 3b 2 Ẅ -θ -Δ + uc + 1 2 ẇ2 c , (13)
and

ε y = -ν x Δ + ∂w c ∂y . ( 14
)
It should be noted that the first term in Eq. ( 14) is necessary to remove spurious terms in the total potential energy that would arise if the Poisson's ratio effect under pure compression is not considered [6]. The other contribution to the strain energy in the core U cs comes from the shear strain γ xy present, where:

U cs = c 2 L 0 b/2 -b/2 G c γ 2 xy dy dx. ( 15
)
The quantity G c is the shear modulus of the core material and the shear strain γ xy contains both overall and local components:

γ xy = ∂W ∂x -θ R + ∂w c ∂x + ∂u c ∂y . (16)
Note that there are no von Kármán shear strain terms present in this last expression, this is because these are always zero for cylindrical bending cases.

Work done by load

The work done comprises the axial load P multiplied by the total end-shortening E. The contributions come from pure compression, overall buckling and the in-plane displacement due to local buckling:

P E = P L 0 Δ + 1 2 Ẇ 2 - 1 2 u dx. ( 17
)

Total potential energy

Integrating over the depth and assembling all the strain energy and work done contributions, the total potential energy V is given as an integral over

A c c e p t e d m a n u s c r i p t

the length of the strut.

V = U b + U m + U ca + U cs -P E = L 0 D (2q t + q s ) 2 b 2 π 4 9L 2 sin 2 πx L + (2q t + q s ) bπ 2 6L sin πx L ẇ2 + 2 u + 2Δ 2 + 1 4 ẇ4 + u ẇ2 -2 uΔ -ẇ2 Δ + u2 + EI 2 q 2 s 2π 4 L 2 sin 2 πx L + ẅ2 + C x 1260 5q 2 s + 32q s q t + 68q 2 t b 2 π 4 L 2 sin 2 πx L + 315 u ẇ2 + 420 u2 -ẇ2 Δ -21 (q s + 4q t ) bπ 2 L sin πx L ẇ2 + 2 u + 1260(Δ 2 -u) + 63 ẇ4 + C y ν x ν x Δ u + 1 3 ẇ2 -Δ - uw b Δ - 1 3b w ẇ2 + 1 2 kw 2 + G (q s -q t ) 2 8π 2 15 cos 2 πx L + (q s -q t ) 2π 3 ẇ - 2u b cos πx L + 1 3 ẇ2 - u ẇ b + u 2 b 2 -P Δ + q 2 s π 2 2 cos 2 πx L - 1 2 u dx. (18)
For brevity, some of the material and geometric property terms have been combined into the quantities given below:

C x = E x bc 2 (1 -ν x ν y ) , C y = E y bc 2 (1 -ν x ν y ) , k = 2C y b 2 , G = G c bc 2 . ( 19
)
The total potential energy is minimized using the calculus of variations [START_REF] Fox | An introduction to the calculus of variations[END_REF] by employing the procedure detailed in Hunt and Wadee [START_REF] Hunt | Localization and mode interaction in sandwich structures[END_REF]. The first variation of V vanishes by assuming simply supported boundary conditions for w where w(0) = ẅ(0) = w(L) = ẅ(L) = 0, while for u more complicated conditions need to be invoked by matching the applied stress at the ends; for instance, the condition for x = 0 is:

u(0) 2D + 2C x 3 + ẇ2 (0) D + C x 4 -Δ 2D + C x -C y ν 2 x + P 2 = 0, ( 20
)
and a similar condition exists for x = L.

More importantly, the minimization yields a system of two non-autonomous coupled ordinary differential equations (ODEs)-a fourth-order ODE in w and a second-order ODE in u-along with three integral equations that are obtained by differentiating V with respect to the generalized coordinates q s , A c c e p t e d m a n u s c r i p t q t and Δ respectively and setting those expressions to zero:

EI .... w + G (q s -q t ) 2π 2 3L sin πx L + u b - 2 3 ẅ + 2C y 3 ν x ẇ2 b -ẅ ν x Δ - w b + D (2q t + q s ) bπ 2 3L sin πx L ẅ + π L cos πx L ẇ + 2 ẅΔ -2 ( ẇü + ẅ u) -3 ẅ ẇ2 + C x 2 3 ẅ - 3 5 ẇ2 ẅ - Δ 2 ( ẇü + ẅ u) + (q s + 4q t ) bπ 2 30L sin πx L ẅ + ẇ π L cos πx L - C y ν x b u + ẇ2 3 + kw = 0, (21)
2D + 2 3 C x ü + D 2 ẇ ẅ -(q s + 2q t ) bπ 3 3L 2 cos πx L - C y b ν x ẇ + C x 1 2 ẇ ẅ - bπ 3 30L 2 (q s + 4q t ) cos πx L + G b ẇ + 4 3π (q s -q t ) cos πx L - 2u b = 0, ( 22
) L 0 2Gπ 3 cos πx L ẇ - 2u b - D 6 + C x 60 bπ 2 L ẇ2 + 2 u sin πx L dx + EIπ 4 L + Db 2 π 4 18L (2q t + q s ) + C x b 2 π 4 1260L (16q t + 5q s ) + 8GLπ 2 15 (q s -q t ) - P Lπ 2 2 q s = 0, (23)
L 0 D 3 + C x 15 bπ 2 L sin πx L 2 u + ẇ2 dx + Db 2 π 4 9L (2q t + q s ) + C x b 2 π 4 315L (17q t + 4q s ) + 8GLπ 2 15 (q t -q s ) = 0, (24)
P = Δ 4D + 2C x -2C y ν 2 x + L 0 C y ν 2 x -C x u + ẇ2 3 -2D u + ẇ2 2 dx. ( 25
)
Assuming that the overall mode occurs first, by setting w(x) and u(x) to zero, a linear eigenvalue analysis can be performed that yields the critical load.

The following expressions from the RBT model differ significantly from the equivalent expressions from the TBT model owing to the presence of the sway mode W in the membrane energy:

qs 1 - 5 Dπ 2 φ 24 G - Cx π 2 φ 42 G = qt 1 + 5 Dπ 2 φ 12 G + 17 Cx π 2 φ 168 G , (26)
P C = π 2 8 Gφ + Dπ 2 φ 18 G + Cx π 2 φ 252 G + 8 15 + A 1 A 2 Dπ 2 φ 9 G + 4 Cx π 2 φ 315 G - 8 15 , (27)

A c c e p t e d m a n u s c r i p t

where:

A 1 = 1 - 5 Dπ 2 φ 24 G - Cx π 2 φ 42 G , A 2 = 1 + 5 Dπ 2 φ 12 G + 17 Cx π 2 φ 168 G . ( 28
)
The terms with tildes are dimensionless quantities of their respective terms, thus:

qs = bπ 2 L q s , qt = bπ 2 L q t , D = DL 2 8EI , G = GL 3 8EIb , P = P 2G , Cx = C x L 2 8EI , Cy = C y L 2 8EI , φ = b L . (29)
The system of nonlinear equations for the RBT model are solved using the same software as for the TBT model namely the numerical continuation package Auto97 [START_REF] Doedel | AUTO97: Continuation and bifurcation software for ordinary differential equations[END_REF], which is well-known for its capability to pinpoint bifurcation points and plot out multiple branching paths as model parameters are varied. To reduce the computational cost, symmetry in the physical system is exploited in solving the equations for only half the length of the strut. This is achieved by employing the following conditions at midspan:

ẇ(L/2) = ... w(L/2) = u(L/2) = 0. ( 30
)
Obviously as a result of the above conditions any possible asymmetric solutions of the nonlinear ordinary differential equations are automatically neglected. This is not alarming as previous studies [START_REF] Wadee | Effects of periodic and localized imperfections on struts on nonlinear foundations and compression sandwich panels[END_REF] have shown that the symmetric solution has the lowest load at which interactive buckling occurs. Applying continuation principles from the critical point (P = P C ) by introducing a small perturbation to the amplitude of tilt, q t , a secondary bifurcation is identified which leads to buckle localization. Once identified, continuation is applied again, this time by varying the axial load P to trace out the new equilibrium path of secondary buckling. Care needs to be exercised when searching for the secondary bifurcation since if the step-size is too large, there is a possibility that the solver might jump over the most severe secondary buckling mode and as a result a higher energy localized buckling mode may be captured. This is an important feature of this class of nonlinear model; it is well known that the equations of a post-buckled axially-loaded cylindrical shell [START_REF] Lord | Computation of localized postbuckling in long axially-compressed cylindrical shells[END_REF] and a strut on a softening foundation [START_REF] Hunt | Comparative lagrangian formulations for localized buckling[END_REF][START_REF] Champneys | Bifurcation of a plethora of large amplitude homoclinic orbits for Hamiltonian systems[END_REF] have a multiplicity of possible solutions, only one of which is the energy minimizer and hence the practically observed solution. Once the practical localized mode is found, homotopy-the use of numerical continuation in varying physical parameters-is used extensively to find post-buckling equilibrium solutions for different strut configurations.

A c c e p t e d m a n u s c r i p t 4 Results and Discussion

Material and geometrical properties

The set of properties for the sandwich strut is given below and is in accordance with Hunt and Wadee [START_REF] Hunt | Localization and mode interaction in sandwich structures[END_REF] as a practical configuration for a strut in order to draw comparisons with the previous models reported in the literature. The strut width c can be cancelled out from the governing equations without any loss in accuracy, hence the width c has been taken to be unity and the dimensional axial forces are given as a force per unit width in the results below. The length of the strut, L, and the core depth b are left as the principal parameters to vary.

Equilibrium paths

The equilibrium load versus end-shortening paths of the two analytical models for the same strut are shown in Figure 5(a), plotting the normalized load P/P C against the normalized end-shortening, E/L. From these results it is evident that the two models are well correlated with minor differences in the initial stiffness. More importantly, the RBT model predicts the secondary bifurcation before the TBT model, which can be attributed to a more flexible representation of the cross-sectional deformation. This greater flexibility in shear leads indirectly to larger axial stresses in the face plates, which exceeds the compressive stress needed to cause local buckling at midspan significantly earlier than the equivalent TBT model. As a result of the increase in the axial strains in the face plates, there is a redistribution of shear strain through the depth of the strut as shown in Figure 5(b). Having zero shear strain at the edges of the core might not be entirely accurate but a good approximation for the case of thin face plates, as presented in the current formulation, especially where the local centroid of the face plates is taken at the edge of the core. The higher order theory for core bending is more suitable for sandwich panels with deeper cores as it better accounts for the local nonlinear deformations of the cross-sectional planes in the neighbourhood of the extreme fibres, which can be seen in physical experiments [START_REF] Wadee | Experimental evaluation of interactive buckle localization in compression sandwich panels[END_REF] and from finite element simulations presented below (see Figure 11).

Interactive buckling

Beyond the secondary bifurcation, a localized buckle akin to wrinkling appears purely on the bottom face plate at the location of maximum compressive stress. This effect manifests itself as w(x) and u(x) become non-zero and grow as the axial load carrying capacity decreases from P C , as shown in Figure 6. The evolution of the two local modes for the RBT model is similar to the one observed in the results of the TBT model, which is triggered at a slightly lower critical load but at a larger amplitude of sway q s . The modes of the two analytical models are compared at different stages in the post-buckling path (Figure 7). During the early stages of post-buckling, both the amplitude and the so-called "wavelength of localization" λ-the distance between the first extrema on either side of midspan, see Figure 8 [10]-are smaller for TBT compared to the equivalent values for RBT. Nevertheless, as the path evolves the difference in the modes from the two models is reduced as both the wavelength and the maximum amplitudes of w converge; a similar trend is followed by the u(x) mode. An explanation for the results of the two models converging at the advanced post-buckling state is probably due to the effect of the large amount of lateral deflection, attributed to overall buckling, which then dominates the behaviour. The growth of the sway amplitude q s reduces the relative effect of the nonlinear in-plane deformation field-the main difference in the modelling-which grows at a much smaller rate. This progressively reduces the relative difference in the response between the two models.

Parametric study

Length variation

Increasing the length L of the strut configuration given in §4.1 increases the space between the critical and secondary bifurcations as overall buckling becomes dominant. For lengths ranging from 100 mm to 508 mm (φ = 0.0508 → 0.0100) localized buckling is triggered at a larger normalized end-shortening with increasing levels of snap-back being observed as the aspect ratio φ is reduced. The wavelength of localization λ remains relatively constant. As the total length increases the buckle appears more localized; an extra wave peak or trough is observed giving rise to a quasi-linear relationship between the length and the number of peaks and troughs clearly visible in the localized mode, see Figure 9(a).

Core depth variation

The core depth is probably the most informative parameter to vary in the current study, as the objective of developing a higher order theory model such as RBT is to account for nonlinear cross-sectional deformations. As the depth increases the effect of the higher order deformation, the assumption of plane sections remaining plane is less likely to hold even approximately. This is applicable for softer cores where the axial stresses are principally concentrated at the much stiffer extreme fibres and comparatively less strain energy is stored in the core.

Using the previous strut configuration for three different depths (b = 5.6 mm, 6.6 mm and 10.2 mm), the interaction between the overall and local mode is examined. The pattern shown in Figure 9(b) clearly indicates the increase of the sway mode at the secondary bifurcation, with increasing depth. In all cases, comparing the RBT model with the TBT model, the amount of the sway mode q s before the secondary bifurcation is triggered is again always smaller and the subsequent evolution of the localized mode occurs earlier. The graph also shows that the comparison between RBT and TBT becomes worse as the core depth becomes relatively larger; this is to be expected since the RBT model should be superior for deeper cores owing to its improved modelling of the cross-section distortion.

Validation

The results of the current research work are based on analytical formulations of perfectly elastic struts. However, no provision for delamination and hence no discontinuity in the stress distribution at the face-core interface is considered. The analytical models have been validated by means of the Finite Element (FE) method using the general purpose commercial software Abaqus [START_REF] Abaqus | ABAQUS/Standard: User's manual version 6[END_REF]. The sandwich strut with the geometric and material properties given in §4.1 was modelled as a 2D continuum with the plane stress assumption. The core material was discretized with the two-dimensional solid element CPS4R, which is a 4-noded bilinear element with reduced integration and hourglass control. Mesh convergence studies revealed that ten elements through the core thickness with an aspect ratio of unity accurately capture the structural response

A c c e p t e d m a n u s c r i p t

of the strut. Regarding the discretization of the face plates, two alternative modelling approaches were considered. The first approach involved modelling the face plates with plane stress solid elements, while the second approach consisted of modelling the face plates as stringers bonded to the edges of the existing core while specifying suitable engineering properties. The 2D linear Timoshenko beam element B22 was used to discretize the stringers.

Although the first approach is conceptually simpler and easier to implement, it is computationally expensive as several layers of plane stress solid elements are needed to discretize each face plate in order to avoid shear locking. Owing to the small thickness of a face plate (t = 0.5 mm), in conjunction with the large length and the need to maintain a reasonable element aspect ratio, the first approach results in an unduly large number of elements. The second approach is less computationally demanding as only one beam element through thickness of the face plate is used without compromising accuracy and has therefore been implemented currently.

As in the analytical model, symmetry has been utilized and only half the length of the strut has been modelled in the FE simulations by applying suitable boundary conditions along the axis of symmetry at the midspan of the strut. Simple pinned conditions were simulated by restraining the degree of freedom perpendicular to the strut's axis at the mid-height of the loaded edge. Displacement control was utilized to apply loading by prescribing a displacement at all the nodes of the loaded edge of the strut; the corresponding axial force was derived from the reaction forces at the strut's mid-height, where symmetry boundary conditions were applied. No further constraints were applied to the model; the lateral expansion of the cross-sections due to the Poisson's ratio effect and the non-planar deformation of the strut's crosssections were free to develop.

It is well known that FE packages cannot analyse perfect systems in the postbuckling range without a small perturbation [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]; hence very small imperfections in the form of the lowest buckling mode shape were incorporated in order to give the closest possible comparison to the perfect case from the analytical model. A linear eigenvalue analysis was initially conducted to extract the lowest buckling mode shape. This was thereafter used as an initial perturbation in the model geometry in the subsequent geometrically nonlinear analysis, employing the modified Riks method [START_REF] Abaqus | ABAQUS/Standard: User's manual version 6[END_REF]. An imperfection amplitude of L/10000 was selected as it was sufficient to trigger instability, yet small enough to give results close to the perfect case. It should be noted that although no local imperfection was explicitly incorporated in the model, a secondary instability and buckle pattern localization was subsequently observed from a self-generated local imperfection that evolved during the post-critical analysis.

A c c e p t e d m a n u s c r i p t

Comparison of the finite element results with the Auto97 results of the RBT model reveals excellent agreement at the full range of the elastic response including initial stiffness, the critical load, the triggering of localized buckling and the post-buckling path as shown in Figure 10. The TBT result is also

Critical loads

A common feature of all the analytical models is the underestimation of the critical load at the critical bifurcation as established by the method described in Allen [START_REF] Allen | Analysis and design of structural sandwich panels[END_REF] for thin face plates and weak cores.

P C Allen = P E 1 + (P E /P S ) , (31)
where:

P E = Dπ 2 L 2 (b + t) 2 , P S = G c c b (b + t) 2 . ( 32
)
The reason for this is that the lever arm that was used to calculate the axial strains in the faces was taken to be b/2, which simplifies the formulation considerably but effectively calculates the strains at the face-core interface rather than the local face plate neutral surface. For a strut of the above configuration (b/t = 10) and L = 200 mm, the difference in P C is approximately 14% while for a configuration approaching to b/t = 20 the difference drops to less than 4% (see Table 1). The model could be improved to evaluate the critical load more accurately by altering the in-plane tilt displacement field, u R (x, y) to: where y ranges from:

u R (x, y) = -y θ(x) + 4y 2 3 (b + t) 2 ∂W ∂x -θ(x) , (33)
-(b + t)/2 y (b + t)/2. ( 34
)
Changes would also need to be made to the vertical distribution of the core in-plane deformation functions, w c and u c , leading to longer and more complicated coefficients in the ordinary differential equations. However, as the objective of the analytical approach was to capture the localized post-buckling response, this change is deemed not to be necessary even though it is appreciated that the applicability of the model depends to some extent on the level of acceptable percentage difference in the estimation of the critical load. Another option is to consider that the critical load results provided by the model are more suitable for a strut with a core of effective thickness b , where b = b -t, since the axial stresses in the core are actually very small and the reduction in the total strain energy would be minimal. Such an approach would reduce the error further, effectively matching the Allen's critical load. It should be noted that the second FE modelling approach described in the previous subsectionmodelling the face plates with stringers-is actually in accordance with the assumptions regarding the effect of the face thickness made in the analytical model; this partly explains the excellence in the comparisons between the FE and TBT models.

Interactive buckling

In addition, the comparison between the RBT and the FE models for the localized buckling modes of the bottom face plate with respect to the localized wavelength and the lateral displacement at midspan is excellent. This correlation improves for deeper beams as observed in Figure 12, highlighting the strength of the RBT model for deeper sandwich struts in particular. The deflected shapes of the more compressed face plate in Figure 12 were compared at different stages of post-buckling. The effective match in the amplitude of localization between the RBT and FE models at early post-buckling can be attributed to the close proximity of the secondary bifurcation points. The secondary bifurcation for the TBT model is found to occur at a lower load and therefore the localized buckling mode initially has a much smaller amplitude, see Figures 12(b) and 12(e), when compared to both the RBT and FE models. Further along the post-buckling path, the results for the analytical models have more significant numbers of peaks and troughs away from midspan. This introduction of more peaks and troughs that spread towards the ends of the strut produces a slightly stiffer equilibrium path when compared to the FE model. Conversely, the FE model shows the localization being confined to a limited number of peaks and troughs that only increase in amplitude, not in number. The difference between the post-buckling paths of the FE and RBT models is evaluated by comparing the load levels of each path and quantifying their difference as a percentage of the critical load. It has been found that the overestimation, or difference, of the load level at various stages in post-buckling is close to 3% in early post-buckling (0.9 P/P C 1.0), whilst further down the post-buckling path (P/P C ≈ 0.6), where the strut has lost a significant proportion of its original strength, the error increases to approximately 5%, which remains an excellent comparison.

A c c e p t e d m a n u s c r i p t

A rational explanation for the slightly softer response of the FE model in the advanced post-buckling state is connected with the fact that the change in the local cross-sections is accounted for in the mechanical response, whereas the analytical models assume that the core depth remains at b. To account for this in the analytical model would be quite cumbersome as the limits of integration in y for the core energy and the lever arm for the membrane energy would have to be adjusted to account for the displacement w and the commensurate change in the position of the cross-section neutral axis. This would introduce a considerable number of new nonlinear terms that may improve the comparison even more, but it is arguable whether any gains in accuracy would be only marginal given that there are probably more important effects to consider such as core material nonlinearities and plasticity in the faces or the core.

Concluding Remarks

An analytical model with higher order bending theories for the interactive buckling of sandwich struts was developed. This was achieved by relaxing the constraint of plane sections remaining plane under bending and allowing the presence of a nonlinear shear strain distribution over the depth of the strut. The model was compared against an existing analytical model that uses the Timoshenko beam approach and validated against a nonlinear finite element formulation. The validation shows excellent agreement between the new analytical formulation and the purely numerical model for the full range of elastic behaviour from critical buckling through to the far post-buckling range.

The results of the comparative study between the previous and current analytical formulations have shown that the current model significantly improves the estimation of the critical load and, more importantly, provides a better prediction within the nonlinear buckling range. The location of the secondary instability that leads to interactive buckling, the resulting post-buckling path and the mode profile that localizes on the more compressed face plate after overall buckling all show that the current model surpasses the previous one in accuracy. It is worth noting though that the difference in the responses of the two models is more important in the early stages of post-buckling; the superiority of the current model is greater for struts with larger core depths and those with a larger difference in the elastic moduli between the face plates and the core.

Work is now continuing to develop the higher-order analytical model further to incorporate combinations of axial force and moment, the effect of imperfections and the effects of using of functionally graded core materials.

Fig. 1 .

 1 Fig. 1. A sandwich panel under axial compression shown in photographs (1)-(3): left to right and response graph of equilibrium paths of load P versus end-shortening E. Photo (1) pre-buckling-path (a) on graph; photo (2) overall buckling-point C and path (b) on graph; photo (3) interaction between overall and local buckling leading to localization-path(c) on graph. Path (d) on the graph represents the imperfect response, S and l mark the location of the secondary bifurcation for perfect geometries and the limit point for imperfect geometries respectively.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Outline of the three bending theories. (a) EBT, (b) TBT, and (c) RBT

Fig. 4 .

 4 Fig. 4. Modal descriptions. Overall-top: sway mode W (x); middle: tilt mode θ(x).Local-bottom: u(x) and w(x) becoming non-zero beyond the second bifurcation. The generalized coordinate Δ gives the purely compressive strain arising from the axial load P .

Face

  plate Young's Modulus: E = 68900 N/mm 2 , Face plate Poisson's ratio: ν = 0.3, Face plate thickness: t = 0.5 mm, Core Young's modulus: (in x and y) E x = E y = 199 N/mm 2 , Core Poisson's ratio: (in x and y) ν x = ν y = 0.2, Core Shear modulus: G c = 83 N/mm 2 , Core depth (range): b = 5.1 mm to 10.2 mm, Strut length (range): L = 100 mm to 508 mm.

Fig. 5 .

 5 Fig. 5. (a) Equilibrium paths for a strut of length 100 mm and depth 5.08 mm for the two models. (b) Shear strain distribution through the core depth of the two models at the same load in the post-buckling range (P = 0.9P C ), as expected the RBT model shows a parabolic distribution of γ xy with depth as opposed to the constant shear strain given by the TBT model. Note that the dashed lines represent the TBT model and the solid lines represent the RBT model.

A c c e p t e d m a n u s c r i p tFig. 6 .Fig. 7 .

 67 Fig. 6. Mode evolution versus the falling axial load for (a) w(x) and (b) u(x). In this case: L = 100 mm and b = 5.08 mm.

Fig. 8 .

 8 Fig. 8. Definition of localized buckle wavelength λ and maximum wave height H.

Fig. 9 .

 9 Fig. 9. (a) The quasi-linear relationship between the number of half-sine waves in localized buckling and the length for the sandwich strut with b = 5.08 mm. (b)Local versus overall mode for sandwich struts (L = 100 mm) with varying depth showing a comparison in the response between TBT (dashed line) and RBT (solid line) models. The maximum wave height H is defined in Figure8and the secondary bifurcation occurs where the curves intersect the abscissa.

Fig. 10 .

 10 Fig. 10. The equilibrium paths for a strut of length 100 mm for the two models and FE simulation: (a)-(b) strut depth 5.1 mm; (c)-(d) strut depth 10.2 mm. Note that W max is given by q s L, see Equation (1).

A c c e p t e d m a n u s c r i p tFig. 11 .

 11 Fig. 11. A snapshot and close-up of the Abaqus model of a sandwich panel of length 100 mm and depth 5.1 mm during unloading exhibiting interactive buckling (P = 0.6P C ). Note the non-planar profile of the originally plane sections in the localized buckling region.

A c c e p t e d m a n u s c r i p t

  Percentage difference in P C L/b = 39.4 L/b = 28.1 L/b = 19.7 compared to Allen (1969) b/t = 10 b/t = 14 b/t = 20

Fig. 12 .

 12 Fig. 12. The local mode at different stages in post-buckling for a strut of length 100 mm and depth of 5.1 mm for (a)-(c) and 10.2 mm for (d)-(f). (a)P/P C = 0.95, (b) P/P C = 0.9, (c) P/P C = 0.7, (d) P/P C = 0.9, (e) P/P C = 0.85, (f) P/P C = 0.7. TBT is absent in (a) and (d) as its localized mode is triggered at a lower load level. The results for the RBT model is represented by the thick solid line, the TBT model by the dashed line and the FE model by the thin solid line.

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

Preprint submitted to Int J. Non-Linear Mech.

Acknowledgement

Financial support for S.Y. was provided by a UK Engineering and Physical Sciences Research Council (EPSRC) Doctoral Training Award.


P Grammenoudis

Ch Tsakmakis
email: tsakmakis@mechanik.tu-darmstadt.de

D Hofer

Micromorphic continuum. Part III: Small deformation plasticity coupled with damage

Keywords: micromorphic plasticity, finite deformation, isotropic and kinematic hardening rules, isotropic continuum damage mechanics, finite element implementation

Properties of the micromorphic theory proposed in Part II are discussed for the case of small deformations. Model responses for beam specimens under bending loading and plates with circular holes under tension loading are calculated by employing the finite element method. The results reported are concerned with the capabilities of the theory to predict size effects.

A c c e p t e d m a n u s c r i p t 1 Introduction

A finite deformation micromorphic plasticity theory, exhibiting isotropic and kinematic hardening, and incorporating damage effects, has been proposed in Part II. The theory is consistent with the second law of thermodynamics and deals with a plastic micromorphic curvature tensor, which is not required to fulfill some compatibility conditions, i.e. it is not related to some gradient terms. Furthermore, a measure of smallness ε has been introduced in Part II, and the theory has been defined to be of small deformations, if terms only up to order O(ε) are retained. It is a straightforward task to verify that the small deformation version of the micromorphic model proposed, reads as follows (we confine attention to static balance equations and omit the body and double body forces).

Equilibrium equations

∂T ij ∂X j = 0 in R R , (1)
∂T ijk ∂X k + T ij -Σ ij = 0 in R R . (2)
Boundary conditions

T ij n j = ti on ∂R t i R = ∂R R \ ∂R u i R , (3)
T ijk n k = t(d) ij on ∂R t (d) ij R = ∂R R \ ∂R h ij R , (4)
u i = ūi on ∂R u i R , (5)
A c c e p t e d m a n u s c r i p t

h ij = hij on ∂R h ij R . (6)
Kinematics

H ij = ∂u i ∂X j , (7)
β ij = 1 2 (h ij + h ji ) , ij = H ij -h ij , K ijk = ∂h ij ∂X k , (8)
β ij = (β e ) ij + (β p ) ij , ij = ( e ) ij + ( p ) ij , K ijk = (K e ) ijk + (K p ) ijk . (9)
Specific free energy

Ψ = Ψ e + Ψ is + Ψ k . ( 10
)
Elasticity laws

Ψ e =(1 -D) 1
T ijk = ∂Ψ e ∂(K e ) ijk = (1 -D)(C e ) ijkpqr (K e ) pqr . (13)
Yield function

f = 1 1 -D (T ij -T k ij ) D (A A A y ) ijpq (T pq -T k pq ) D + (Σ ij -Σ k ij ) D (B B B y ) ijpq (Σ pq -Σ k pq ) D +(T ijk -T k ijk ) D (C y ) ijkpqr (T pqr -T k pqr ) D -1 2 - R 1 -D -k0 , (15)
k := R 1 -D + k0 , k := R 0 + k0 . (16)

A c c e p t e d m a n u s c r i p t

Flow rule

( ˙ p ) ij = ṡ ζ ∂f ∂T ij , ( βp ) ij = ṡ ζ ∂f ∂Σ ij , ( Kp ) ijk = ṡ ζ ∂f ∂T ijk , (17)
ζ := ∂f ∂T ij ∂f ∂T ij + ∂f ∂Σ ij ∂f ∂Σ ij + ∂f ∂T ijk ∂f ∂T ijk . ( 18
)
Plasticity

L(t) := [ ḟ(t)] s=const. ( 19
) ṡ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ > 0 for f = 0 & L > 0 , = 0 otherwise , (20)
ṡ : to be determined from consistency condition ḟ = 0 .

Viscoplasticity

ṡ := f m η ≥ 0 , (22)
f : overstress . ( 23
)
Isotropic hardening

Ψ is = (1 -D) γ 2 (r 2 + 2r 0 r) , (24)
R = ∂Ψ is ∂r = (1 -D)γ(r + r 0 ) = (1 -D)(γr + R 0 ) , (25)
ṙ = (1 -βr) ṡ ζ . ( 26
)
Kinematic hardening

Ψ k =(1 -D) 1 2 (A A A k ) ijpq ( k ) ij ( k ) pq + 1 2 (B B B k ) ijpq (β k ) ij (β k ) pq +(D D D k ) ijpq ( k ) ij (β k ) pq + 1 2 (C k ) ijpqr (K k ) ijk (K k ) pqr , (27)
A c c e p t e d m a n u s c r i p t

(Σ k ) ij = ∂Ψ k ∂(β k ) ij = (1 -D){(B B B k ) ijpq (β k ) pq + (D D D k ) ijpq ( k ) pq } , ( 28
) (T k ) ij = ∂Ψ k ∂( k ) ij = (1 -D){(A A A k ) ijpq ( k ) pq + (D D D k ) ijpq (β k ) pq } , (29)
(T k ) ijk = ∂Ψ k ∂(K k ) ijk = (1 -D)(C k ) ijkpqr (K k ) pqr , (30)
˙ k = ˙ p - ṡ 1 -D {M k 1 (trT k )1 + M k 2 T k + M k 3 (T k ) T } , ( 31
) βk = βp - ṡ 1 -D {N k 1 (trΣ k )1 + 2N k 2 Σ k } , (32)
Kk = Kp - ṡ 1 -D Pk [T ] . ( 33
)
Evolution law for damage

Ḋ = -α 1 ṡ ∂Ψ ∂D . ( 34
)
The aim of the present paper is to demonstrate the capabilities of this model in describing size effects present in bending of beam specimens and in plates with a hole under tension loading. It should be remarked that Part I and II made it clear, that micromorphic constitutive theories are very complex and include a large number of material parameters. Therefore, we decided to make transparent capabilities of such theories only for small deformations, excluding from considerations geometrical nonlinearities. Also, several material parameters will be assumed to vanish, in order to reduce the effort of the analysis. Of course, this implies that important capabilities of the model may be not activated. However, the present investigation is not entitled to be complete and will be of qualitative character only. This also concerns the isotropic hardening rule. In fact, isotropic hardening effects due to strains and micromorphic curvature tensors are captured in a unified manner. There are, however, possibilities to account for isotropic hardening effects due to strain and micromorphic curvature effects separately. Such isotropic hardening rules

A c c e p t e d m a n u s c r i p t

have been elaborated by [START_REF] Grammenoudis | Isotropic hardening in micropolar plasticity[END_REF] in micropolar plasticity and are not pursuit here.

Examples

Examples illustrating the capabilities of the theory to capture size effects are given in this section and are taken from the doctoral thesis of [START_REF] Hofer | Simulation von größeneffekten mit mikromorphen theorien[END_REF],

where also more details about the implementation are given. Further examples and interesting results on this topic may be found in [START_REF] Dillard | Micromorphic continuum modelling of the deformation and fracture behaviour of nickel foams[END_REF]; [START_REF] Kirchner | A unifying treatise on variational principles for gradient and micromorphic continua[END_REF]; Neff and Forest (2007); Lazar and Maugin (2007) as well as [START_REF] Hirschberger | On deformational and configurational mechanics of micromorphic hyperelasticity -theory and computation[END_REF].

In the ensuing analysis, the chosen values of the material parameters do not reflect some responses of realistic material behavior, i.e., they are only of academic interest and serve to discuss basic features of the model. We set

A e 1 ≡ λ = 1, 21 • 10 5 N/mm 2 , A e 2 = μ + α , A e 3 = μ -α , (35)
μ = 8, 08 • 10 4 N/mm 2 , (36)
B e 1 ≡ λ , B e 2 ≡ μ + b 2 , b 2 = 10 • μ , (37)
D e 1 ≡ λ , D e 2 ≡ μ , (38)
C e i = 0 for i = 7 , C e 7 = c 7 ≥ 0 . ( 39
)
Although for the case (39), important aspects of the constitutive model may be retain inactive, we shall confine ourself on this special case in order to limit the discussion. For what follows, of particular interest is the internal length are imposed. The whole circular hole is assumed to be traction free, while the whole boundary is subjected to vanishing double traction t (d) .

l c := c 7 μ , (40)
For small circular hole, a nearly uniform stress component σ 0 in y direction, at y = b 2 , will be required to realize the given boundary conditions. In classical elasticity, attention is focussed on the so-called stress concentration factor

T yy σ 0 , T yy := T yy (x = r, y = 0) , (41)
which turns out to be equal to 3 (see e.g. Gould, 1983, p. 124) The effect of α, c 7 on the stress concentration factor is illustrated in Fig. 3.

A c c e p t e d m a n u s c r i p t

To obtain an insight into the size effects due to different, but otherwise similar boundary value problems, we ask for the stress concentration factor T yy σ 0 for the cases where α ≡ μ and geometry and boundary conditions of the specimens vary from each other according to a factor n = 1, 4, 20, 200. Corresponding results are displayed in Fig. 4, from which we deduce that all distributions are similar. In fact, if the T yy σ 0 values corresponding to the specimen according to factor n are plotted as a function of c 7 μ • n 2 , then all plots will coincide (see Fig. 5). In other words, for linear micromorphic elasticity, size effects may be visualized by varying the parameter c 7 , the other parameters being held fixed.

T * yy / σ 0 [ ] 3 l m /l c [ ] T * yy /σ 0 [ ]

A c c e p t e d m a n u s c r i p t

Further size effects may be elucidated by introducing a typical geometry length, as e.g. l m := 4r = 0, 4b. Again we concentrate ourself on specimen geometries and related boundary conditions, differing according to a factor n, with n being now n = 0, 0001, 0, 01, . . ., 400, 10000. On choosing c 7 μ = 0, 1mm 2 , the internal length l c becomes l c = 0, 31623mm. It can be seen in Figure 6, that the stress concentration factor T yy σ 0 is a function of the ratio l m l c (cf. also Mindlin (1963)).

Displacement controlled loading of cantilever rectangular beam

Further features of micromorphic elasticity may be illustrated with the aid of the cantilever rectangular beam shown in Fig. 7. We use Cartesian coordinates x, y and assume plane strain state to apply, with following boundary conditions,

x = 0 : u y = ūy , t x = 0 , t (d) = 0 , (42) x = l : u = 0 , h = 0 , (43) y = 0 : t = 0 , t (d) = 0 , (44) y = l : t = 0 , t (d) A c c e p t e d m a n u s c r i p t

c 7 /μ [mm 2 ] M c / T (class) xx • l [ ]
0 1 2 3 -2 -1 0 1 2 c 7 = 10 -8 • μ mm 2 x [mm] y [mm] 0 1 2 3 -2 -1 0 1 2 c 7 = 0.1 • μ mm 2 x [mm] y [mm] 0 1 2 3 -2 -1 0 1 2 c 7 = μ mm 2 x [mm] y [mm] 0 1 2 3 -2 -1 0 1 2 c 7 = 5 • μ mm 2 x [mm] y [mm] 0 1 2 3 -2 -1 0 1 2 c 7 = 10 • μ mm 2 x [mm] y [mm] 0 1 2 3 -2 -1 0 1 2 c 7 = 10 8 • μ mm 2 x [mm] y [mm]

A c c e p t e d m a n u s c r i p t 2.2 Micromorphic plasticity coupled with damage

In the following, we set

A y 1 ≡ 0 , A y 2 = 1, 5 , A y 3 = 0 , (46)
B y 1 ≡ 0 , B y 2 = 0 , (47)
C y i ≡ 0 for i = 7 , C y 7 = r 7 = 0 , (48)
k 0 = 350N/mm 2 (49)
in the yield function, and

β = 17 , γ = 4100N/mm 2 (50)
in the rule for isotropic hardening. Moreover, we fix the values of α and c 7 in the elasticity laws by

α = 0, 1 • μ , c 7 μ = 0, 1mm 2 . ( 51
)

Uniaxial loading

First, we present calculations for homogeneous uniaxial tension loading of a rectangular specimen (plane strain), according to Fig. 13. At the bottom of the specimen it is given u y = 0, t x = 0, t (d) = 0, while at the top it is u y = ūy ,

t x = 0, t (d) = 0.
The remaining boundaries are subject to the conditions t = 0 and t (d) = 0. The aim is to demonstrate the capabilities of the damage model.

To this end, it suffices to concentrate on isotropic hardening only. Further, as the deformations are homogeneous, no material parameters of terms related to micromorphic curvature tensors are involved. Fig. 14 shows the effect of the damage parameter α 1 (cf. Eq. ( 34)) on the responses of the uniaxial stress σ We consider again the boundary value problem of Sect. 2.1.1, but now with respect to the specimen geometry displayed in Fig. 15 (length l differs from width b). In order to elucidate the capabilities of the micromorphic theory in predicting size effects, four specimen geometries are considered, referred to as specimens 1,4,20 and 200 (see Tab. 1).

First, only isotropic hardening is addressed, with material parameters as given in Sect. 2.2, and r 7 = 10mm -2 . The discussion is referred to the stress com- Table 1 Specimen geometries. Comparison of Fig. 16 (α 1 = 0, 1) with Fig. 17 Next, we assume the micromorphic model material to exhibit kinematic hardening only, governed by the material parameters

r 7 = 10mm -2 , M k 2 = 50mm 2 /N , A k 2 = A k 3 = 200N/mm 2 , ( 52
) P k 7 = 500/N , C k 7 = 200N/mm , (53)
the remaining material parameters related to kinematical hardening being vanishing. Similar to the case of pure isotropic hardening, Fig. 19 suggests that, softening for large specimens begins earlier than for small ones. Fig. 20 confirms that this holds also for the case of combined isotropic and kinematic hardening. Part III. However, no comparison with experimental data is available, so that it is not possible to evaluate the appropriateness of the chosen constitutive functions. This concerns over all the yield function and the isotropic hardening, the latter being unifiedly postulated. Further studies, with reference to experimental results will help to clarify such issues, but this is beyond of the scope of the present paper. All the discussions in the three articles make clear, that phenomenological micromorphic theories (at least plasticity theories) are very complicated and involve a large number of material parameters. Therefore, it will be useful to clarify in future works, if it is possible to approximate the essential material responses predicted by micromorphic theories by some simpler gradient models, which deal with classical stresses only, and involve a smaller number of material parameters. Also it is of interest to answer the following question. Is the micromorphic model appropriate enough to describe all known size effects to the necessary degree?

  2 (A A A e ) ijpq ( e ) ij ( e ) pq + 1 2 (B B B e ) ijpq (β e ) ij (β e ) pq +(D D D e ) ijpq ( e ) ij (β e ) pq + 1 2 (C e ) ijpqr (K e ) ijk (K e ) pqr ,(11)Σ ij = ∂Ψ e ∂(β e ) ij = (1 -D){(B B B e ) ijpq (β e ) pq + (D D D e ) ijpq ( e ) pq } ,(12)T ij = ∂Ψ e ∂( e ) ij = (1 -D){(A A A e ) ijpq(e ) pq + (D D D e ) ijpq (β e ) pq } ,

  elasticity laws. Firstly, we shall discuss micromorphic elasticity without damage (pure micromorphic elasticity). 2.1 Pure micromorphic elasticity 2.1.1 Rectangular specimens with circular hole under tension loading Consider the plane strain problem in Fig. 1 where the quadratic section (length b) with a circular hole (radius r) located in the center of the section, is stretched in y direction. With respect to the Cartesian coordinate system x, y, the boundaries x = ± b 2 are assumed to be traction-free. At the boundary y = -b 2 the displacement u y and the traction t x are assumed to vanish, while at the boundary y = b 2 , given displacement u y and traction t x = 0

Fig. 1 .

 1 Fig. 1. Plane strain problem. The quadratic section with a circular hole is stretched in y direction. x ≥ r, or equivalently a := xr ≥ 0. For c 7 μ = 0, 1mm 2 this distribution, parameterized by α μ , is shown in Fig. 2. It can be recognized that increasing values of α μ cause decreasing values of T yy σ 0 in the neighborhood of a = 0, and consequently decreasing values of stress concentration factors T yy σ 0 for the micromorphic material. Note that all distributions intersect at a = 0, 13mm.

Fig. 3 .

 3 Fig. 2. Distribution of T yy σ 0 , T yy = T yy (y = 0, a ≥ 0), for c 7 μ = 0, 1mm 2 and varying values α μ .

  Fig. 4. Distributions of T yy σ 0 against c 7 μ , for α = μ and different specimens. Geometry and boundary conditions of the specimens differ by a factor n = 1, 4, 20, 200, the corresponding specimens being referred to as specimen 1, . . . , specimen 200, respectively.

Fig

  Fig. 6. Stress concentration factorT yy σ 0 as a function of the ratio l m l c at α = μ and c 7 μ = 0, 1mm 2 ; left: linear plot, right: semilogarithmic plot.

Fig. 7 .]Fig. 8 .Fig. 9 .

 789 Fig. 7. Displacement controlled loading of a cantilever rectangular beam, l = 3, 4375mm, b = 1, 25mm, ūy = 0, 01mm.

Finally,Fig. 10 .

 10 Fig. 10. Effect of material parameters α, c 7 on the response of stress component T xx for point B (indicating a larger distance than point A from the singular point x = l, y = b). The values of T xx are always smaller than T (class) xx .

Fig. 11 .

 11 Fig. 11. Responses of T xy and T yx (left), as well as M c (right), at point A (α = 1, 0•μ).

Fig. 12 .

 12 Fig. 12. Initial and deformed meshes of the rectangular beam for fixed α = μ and varying material parameter c 7 . Displacements u y are presented enlarged, by factor 100. The classical case is approached for c 7 → 0.

  Fig. 13. Displacement controlled uniaxial loading (1 element, plain strain).

Fig. 15 .

 15 Fig. 15. Displacement controlled tension loading of rectangular sections with circular hole. All stress responses in the following figures are referred to point (a) on the hole (x = b + d 2 , y = l 2 ).

(

  ponent T yy at point (a) (see Fig.15). It can be recognized from Fig.16that softening for large specimens begins earlier than for small ones. size factor n length l width b diameter d

]Fig. 16 .

 16 Fig. 16. Response of the stress component T yy on the hole at point (a) as a function of the global strain Δl l (α 1 = 0.1, r 7 = 10mm -2 ).

]Fig. 17 .

 17 Fig. 17. Response of the stress component T yy on the hole at point (a) as a function of the global strain Δl l (α 1 = 1.0, r 7 = 10mm -2 ).

(α 1

 1 = 1, 0) suggests that the form of the responses is strong dependent on the damage parameter α 1 . Moreover, Fig.18illustrates that maximal values of stresses and maximal global strains depend on the material parameter r 7 , present in the yield function.

]Fig. 18 .

 18 Fig. 18. Effect of the material parameter r 7 on the response of stress component T yy at (a) (specimen 4, α 1 = 0, 1).

]Fig. 19 .]Fig. 20 .

 1920 Fig. 19. Responses of T yy at (a) for pure kinematic hardening (α 1 = 1,k 0 = 500N/mm 2 ).

Weber, G., Anand, L., 1991. Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids.

Comp. Meth. Appl. Mech. Eng. 79, 173-202. Weber, G., Lush, A., Zavaliangos, A., Anand, L., 1990. An objective timeintegration procedure for isotropic rate-independent and rate-dependent elastic-plastic constitutive equations. Int. J. Plasticity 6, 701-744.