On the statistical determination of strength of random heterogeneous materials

Accepted Manuscript letter On the statistical determination of strength of random heterogeneous materials Geng Chen, Alexander Bezold, Christoph Broeckmann, Dieter Weichert PII: DOI: Reference:

S0263-8223(16)30269-0 http://dx.doi.org/10.1016/j.compstruct.2016.04.023 COST 7387

To appear in:

Composite Structures

Please cite this article as: Chen, G., Bezold, A., Broeckmann, C., Weichert, D., On the statistical determination of strength of random heterogeneous materials, Composite Structures (2016), doi: http://dx.doi.org/10.1016/ j.compstruct.2016.04.023

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

On the statistical determination of strength of random heterogeneous materials Geng Chena , Alexander Bezolda , Christoph Broeckmanna,1 , Dieter Weichertb a

Institute of Materials Application in Mechanical Engineering, RWTH Aachen University, Augustinerbach 4, 52062, Aachen, Germany b Institute of General Mechanics, RWTH Aachen University, Templergraben 62, 52062, Aachen, Germany

Abstract This paper introduces a numerically-based methodology for determining the load-bearing capacity of particulate reinforced metal matrix composites (PRMMC) under both monotonic and cyclic loading. A multi-scale approach combining shakedown analysis with homogenization was used to evaluate the influence exerted by the composite structure on the global behavior of PRMMCs. The results were interpreted using statistical methods in order to take account of the randomness associated with the composite structure of the material. The general work flow of the approach can be summarized as follows: first, a large number of representative volume elements (PRMMC) were constructed from real material images using an in-house code. Next, lower bound limit and shakedown problems were solved via the interior-point method. Finally, results were converted to their corresponding macro quantities and evaluated statistically. With this approach we investigated a representative PRMMC structure, WC/Co, with two different binder contents and scrutinized the contribution of the composite structure over composite strength. Keywords: A. Particle-reinforced composites, C. Statistics, B. Strength, C. Finite element analysis (FEA), Direct method (DM)

1

Email address: [email protected] (Christoph Broeckmann) Corresponding Author

Preprint submitted to Composite Structures

April 12, 2016

1. Introduction Mechanical structure failure is often caused by material deterioration accumulated over a large number of load cycles. Therefore, to guarantee structural integrity and to make economic use of material, it is important to know limit strength under both monotonic and long-term variable loading. If the part is made of composite material, it is important to assess the precise strength of the composite material under the prevailing loading conditions as well as the overall behavior of the mechanical structure. This is the subject of this paper. Since composite materials consisting of dissimilar phases can be regarded as mechanical systems on the mesoscale, their overall responses are greatly influenced by the alignment of the reinforcement phase and the underlying material morphology [1, 2, 3]. Understanding how the microstructure contributes to the global material behavior, particularly the macroscopic strengths, is one of the central tasks of currents studies on composites [4]. Some phenomenological failure criteria formulated in terms of stress [5, 6, 7] or strain [7, 8] have been introduced in order to evaluate the strengths of heterogeneous materials. These criteria predict the overall composite strength by evaluating several macroscopic descriptors, most of which require a large set of parameters to be collected from experiments. In addition to these phenomenological criteria, the rapid advances in the field of micromechanics allowed continuum mechanics to be applied directly to the mesoscale and macroscopic parameters to be extracted using homogenization techniques. In the light of such development, limit analysis was introduced to composite materials to study their macroscopic limits against plastic collapse [9, 10, 11]. In parallel to the progress of methods developed for monotonic loading, criteria for predicting fatigue strength have also been introduced in a large body of works [12, 13, 14, 15]. One common feature shared by these phenomenological and analytical methods is that the information about the actual composite structure was omitted. This shortcoming can be compensated by introducing numerical methods, predominantly the finite element (FE) method, to the mesoscale. This technique has been applied in numerous works to study the global elastic-plastic behaviour of the composite under both monotonic loading [16, 17, 18, 19, 20] and cyclic loading [21, 22]. For many engineering problems, only the load bearing capacity of a structure is of interest, rather than the entire load bearing history. In such circum2

stance, the effort required to perform cumbersome step-by-step simulation can be saved by using Direct Methods (DM). The theories behind DM are bounding theorems initially proved by Melan [23] and Koiter [24]. A number of numerical methods have been developed on the basis of these theorems in recent years, including linear programming with linearized yield surface [25], basis reduction [26, 27, 28], selective algorithm [29, 30], linear matching method (LMM) [31, 32], residual stress decomposition method (RSDM) [33], nonlinear programing (NP) [34, 35, 36, 37, 38]. In addition to being used to study conventional structural problems on a macroscale, many of these approaches were used also in conjunction with numerical homogenization techniques to predict the strengths of heterogeneous materials. The contributions made in this respect can refer to following studies: for instance, [39, 40, 41, 42, 43, 44, 45] in which where the static theorem was applied, and [46, 47, 48, 49, 50] where the kinematic theorem was applied. It is important to note that in these studies, the materials to be considered were confined to periodic heterogeneous materials with regular geometrical arrangement of constituents. Numerous pragmatic issues have to be addressed prior to expanding the scope of the study to encompass a wider range of materials such as particulate-reinforced metal matrix composites (PRMMCs) which often have a random and non-periodic microstructure. When numerical micromechanical approaches are applied to non-periodic random materials, one of the first steps is to determine a minimum size of the representative volume element (RVE) which explicitly takes account of all features of the microstructure [51, 18]. Although the problem of size effect and its solution has been extensively investigated in many studies [52, 53, 54, 55, 56, 57, 58], the question as to whether an RVE exists in a general sense [59] remains an open question. One possible solution to this problem is to replace an individual RVE by a number of statistical RVEs and view the effective properties from a statistical perspective [60, 61, 62, 63, 64]. This way, the goal of predicting the macroscopic material behavior can be achieved by evaluating statistical measures extracted from numerous model samples. One of the intrinsic difficulties in using SRVEs is that many wellestablished modeling techniques [65, 66, 67] are no longer suitable due either to the demand for user involvement or to the overall computational effort required. Alongside the difficulty raised by modeling RVE geometry, the computational power required for solving associated structural and NP problems present further challenges. To summarize, before DM can be used to study random PRMMCs with complex microstructures, it will be neces3

sary to find an optimal way to build RVE models and rely on an efficient algorithm to solve large scale NPs. These two tasks can be performed successively by applying the automatic model building technique [68] and the interior point method. The interior point method was first introduced by Karmarkar in the mid-1980s [69]. The method is insensitive to the problem scale and provides brilliant superlinear convergence [70, 71, 72]. This algorithm has been successfully introduced to DM by means of both problem-tailored code [73, 74] and customized problem formation handled by general-purpose solvers [75, 76]. In our previous study [77], based on the elaborated numerical scheme, we studied the overall plastic limit of WC reinforced with 30 Wt.% Co using 2D plane strain RVEs. The scope of the work has been expanded in the present study to include the following aspects: first, 2D plane strain models are replaced by thin-layered 3D models called 2.5D models to enhance objectivity. Next, the study predicts the shakedown limit of materials under both proportionally and disproportionately varied macro loadings. Finally, for each material composition, 50 RVE samples are evaluated and the statistical impact of composite structures are analysed in detail.

2. Statistical lower bound shakedown analysis for heterogeneous materials 2.1. Homogenization method In homogenization theory, the material is considered in two well-separated scales: the mesoscopic scale (local scale) prescribed by the coordinate system y is a scale small enough for the heterogeneities to be separately identified. The macroscopic scale (overall scale) prescribed by the coordinate system x is, in contrast, large enough for the heterogeneities to be smeared-out. The relationship between x and y writes y=

x 

(1)

here  is a small scale parameter determining the size of the RVE. It plays an important role in studying the heterogeneous material, especially for nonuniform structures. For a heterogeneous composite, once submitted to an external loading, its

4

microscopic stress field σ in y and its macroscopic counterpart Σ satisfy the relationship Z 1 Σ= σ(y)dV = hσ(y)i . (2) Ω Ω Here h·i stands for the averaging operator, and Ω indicates the RVE domain. A similar relationship also holds between two strain measures Z 1 E= ε(y)dV = hε(y)i . (3) Ω Ω The local strain ε can be further decomposed into two parts: average value E and a fluctuating part ε∗ ε(u) = E + ε = ε∗ (u)

(4)

with the average of ε∗ over the RVE equals zero, hε∗ i = 0 ,

(5)

u = E · y + u∗ ,

(6)

and where u∗ is the fluctuating displacement. When all phases of an RVE are elastic, the overall behavior of the RVE is elastic as well. In this case, Σ and E are correlated by the effective elastic tensor C Σ=C:E. (7) In contrast, if any constituents of the material exhibit plastic behavior, (7) no longer holds. In this case, the actual micro-stress can be divided into two parts: the first occurs if the material behaves perfectly elastically and the second belongs to a self-equilibrated residual stress field σ(y) = D(y) : Σ + ρ(y) .

(8)

Here D represents the stress concentration tensor. It relates the global stress Σ to the local stress in a purely elastic reference body. In the present study, all constituents are von Mises type materials, thus the effective stress is expressed as r 3 0 σeq = σ : σ0 . (9) 2 The apostrophe in (9) indicates the deviatoric part of a stress tensor. 5

2.2. Lower-bound limit and shakedown analysis The lower-bound shakedown analysis used in this paper is based on Melan’s static theorem; limit analysis is considered as a particular case of shakedown where load is enforced to increase monotonically. The original Melan’s theorem is formulated as follows: if there exist a safety factor α > 1 ¯ which, superposed with and a time-independent field of residual stresses ρ, e the pure elastic stresses σ does not exceed the yield condition F at any time t > 0, ¯ σY ) ≤ 0, ∀t , F (ασ e (t) + ρ, (10) with ∇ · ρ¯ = 0 n · ρ¯ = 0

in V , on ∂V .

(11a) (11b)

Here, σ e is the field of the pure elastic stresses satisfying the equilibrium and boundary conditions for given external loading, σY is the yield stress, n is the outer normal on the boundary ∂V . By neglecting the influence of the mass forces, when shakedown problem is considered in the micro scale, the micro stress fields σ e and ρ¯ are required to satisfy following relationships  ∇ · σ e = 0 in V ,      σ e = C : (E + ε∗ ) in V , elas σ e · n anti-periodic on ∂V , P (12)  ∗  u periodic on ∂V ,    hεi = E .  ∇ · ρ¯ = 0 in V , P res (13) ρ¯ · n anti-periodic on ∂V . Here ρ¯ is a self-equilibrium residual stress field. In accordance with the principle of virtual work, it does not contribute to the virtual work of the system Z ¯ δε : ρdV =0.

(14)

V

In (14), δε indicates the virtual strain and is compatible with virtual displacement δu. Using finite element method, (14) can be discretized and converted to: NE N GE X X

¯ = [C]{ρ} ¯ =0 wm |Jm | [Bnm ]T {ρ}

n=1 m=1

6

(15)

Figure 1: Combination of elastic stresses

Here, C is the equilibrium matrix which relates the stress and the nodal forces; N E indicates the number of elements; N GE the number of Gauss points per element; N G = N GE · N E the number of Gauss points in the whole system; w the weight factor of integration points; J the determinant of the Jacobian matrix. For FE models having N K nodes and N DoF degree of freedoms for each node, one obtains C ∈ RN DoF ·N K×6N G with N DoF equals 3 for 3D case, and 2 for 2D case. Equation (15) represents the equilibrium condition (13) in a weak, discretized form. In most cases, C is a matrix with full row rank, this means the solution vector ρ¯ to this linear system is not unique. In shakedown analysis, the aim is to find a ρ¯ which leads to the maximal α. Here, due to the convexity of (10), it can be proven that such ρ¯ is unique. 2.3. Formulation of the mathematical programming For heterogeneous materials, the external loads can be either macroscopic stresses Σ or macroscopic strains E. When the RVE is submitted to N L linear, independently varying macroscopic loads Pˆn (x, t), Pˆn (x, t) can be separated for convenience into a time independent base Pˆ0n (x) characterizing the load pattern and an associated coefficient µn describing the load magnitude. Thus, as shown in Figure 1, any specific load history H can be defined as H(x, t) =

NL X

Pˆn (x, t) =

n=1

NL X

µn (t)Pˆ0n (x) .

(16)

n=1

K¨onig [78] proves that for materials with a convex yield surfaces shakedown will happen over any load path H within a given load domain L if it happens 7

Figure 2: Shakedown under two proportionally varied loads (2P) and two independently varied loads (4P)

over a cyclic load path containing all vertices of L. This means one can formulate shakedown condition (10) only in regards to loading vertices ¯ σY ) ≤ 0, F (ασ e (Pk ) + ρ,

k ∈ [1..N V ]

(17)

This way, shakedown problem writes (PORI ) minimize ¯ ρ,α

subject to

−α NG X

Ci ρ¯i = 0 ,

(18)

i=1 e F (ασik + ρ¯i ) − σY2 i ≤ 0 ∀i ∈ [1, N G] ; k ∈ [1, N V ] .

In the present study, in order to expose the anisotropy of the RVE, the strength is evaluated on multiple directions by following the subsequent procedure (c.f. Figure 1): first orthogonal stresses Σ11 and Σ22 are prescribed alternately on a purely elastic reference RVE, and each of the microscopic e e stress fields σ11 (y) and σ22 (y) is calculated. A combined loading can then be formed as a joint effect of Σ11 and Σ22 by introducing an angle θ. As demonstrated in Figure 2, the strength in three different loading scenarios was considered in the present study and α was calculated with inequality constraints in (18) satisfied at different set of vertices: • Plastic limit (ΣU = α1P · P1 ): satisfied at P1 , • Shakedown limit under proportionally varied stresses (Σ2P ∞ = α2P · P1 ): satisfied at P1 , P2 , 8

• Shakedown limit under disproportionally 4P (Σ∞ = α4P · P1 ): satisfied at P1 , P2 , P3 , P4

varied

stresses

To improve the efficiency of the calculation, Akoa [75] suggested converting the convex quadratic constraints in (18) into Euclidean ball constraints. When this transformation is adopted, the shakedown problem (18) with N V = 2 can be reformulated into an equivalent form −α

minimize

ui1 ,ui2 ,y1 ,α

subject to

NG X

Mi ui1 + N y1 − αw1 = 0 ,

i=1

ui1 − ui2 = γi ,

(19)

kui1 k2 ≤ 1 , kui2 k2 ≤ 1 , i = 1, . . . , N G . Here, α ∈ R, uik ∈ R5 , y1 ∈ RN G , Mi ∈ R3N K×5 , N ∈ R3N K×N GS , w1 ∈ R3N K . It is important note that the linear constraints in (19) vanish when their corresponding degrees of freedom are restrained. When more than 2 load vertices are considered, (19) can be extended by adding uik (k = 3, . . . , N V ) as variables, and introducing the differences between them ui1 as linear equality constraints whilst their Euclidean norms as inequality constraints. (For more details c.f. [75].) 2.4. Solving the optimization problem by the interior-point method In order to obtain the shakedown factor α, the nonlinear programming (19) has to be solved. In the present study, this problem is solved by a general purposed nonlinear solver IPOPT [79, 80] using a parallel linear algebra package PARDISO [81, 82]. The algorithm behind IPOPT is the interiorpoint line-search filter method. IPOPT accepts problem formulations in a general form minimize n ˆ x∈R

subject to

9

ˆ f (x)

(20a)

ˆ =0, c(x) ˆ>0. x

(20b) (20c)

Here, f : Rn → R, and c : Rn → Rm are twice continuously differentiable functions. To avoid the complication of directly dealing with the inequality constraints, IPOPT solves the problem (20) by introducing a gradually decreasing barrier parameter µ > 0 and computing solutions to a sequence of barrier problems ˆ = f (x) ˆ −µ ϕµ (x)

minimize n ˆ x∈R

n X

log xˆi

(21a)

i=1

ˆ =0. c(x)

subject to

(21b)

For a sub-barrier problem with fixed µ, the first-order Karush-Kuhn-Tucker (KKT) condition (21) writes ˆ + ∇f (x)

m X

ˆ −1 e = 0 , ˆ − µX λi ∇ci (x)

(22a)

i=1

ˆ =0, c(x) ˆ ≥ 0) . (x

(22b) (22c)

Here, λ represents a set of lagrangian multipliers associated with the equality ˆ −1 e, then the first-order KKT system can constraint. By denoting z = −µX be expressed as ˆ + ∇f (x)

m X

ˆ −z =0, λi ∇ci (x)

(23a)

ˆ =0, c(x)

(23b)

i=1

ˆ XZe − µe = 0

ˆ z ≥ 0) . (x,

(23c)

ˆ = diag(x), ˆ Z = diag(z), and e = (1, ...1)T . The nonlinear primalwith X dual KKT system (23) can be solved by a Newton-typed iterative approach. ˆ 0 , z0 ≥ 0 obtained from the previous barrier probStarting from a feasible x lem, a search direction (dxkˆ , dλk , dzk ) is calculated from the linearizion of (21) ˆ k , λk , zk ) at point (x    xˆ    ˆ k ) + A(x ˆ k )λk − zk ∇f (x Wk Ak −I dk  ATk 0 . ˆk) c(x 0  dλk  = −  (24) z ˆ k Zk e − µj e Zk 0 0 dk X 10

ˆ k ) is the Jacobian matrix of the constraints, and Wk Here, Ak = ∇c(x ˆ k , λk , zk ) of the Lagrangian function, denotes the Hessian ∇Lxˆxˆ (x ˆ + cT (x)λ ˆ −z Lxˆxˆ = f (x)

(25)

Instead of directly solving the nonsymmetric linear system (24), IPOPT tries to solve a smaller, simplified, symmetric linear system. After a search direction has been obtained by solving this linear system, the step size will be determined by using a line-search variant of Fletcher and Leyffer’s filter method [83]. The sub-barrier problem-solving process (22) is terminated when the user-defined tolerances are met. After that, the barrier parameter µ is updated and the solution to the evolved barrier problem is calculated by the same iterative scheme. It is apparent that with µ ↓ 0, the solution to the barrier problem is exactly the same as the original problem. When implementing IPOPT, the user is required to provide the parameters determining the acceptable tolerances in addition to the gradient of the objective function ∇f , Jacobian matrix of the constraints Ak = ∇c(xk ), and the Hessians of both f and ci . According to the reformulated static problem (19), primal variables can be noted as  T x = (u11 )T . . . (uN G 1 )T (u12 )T . . . (uN G 2 )T (y1 )T α . (26) Gradient of the objective function ∇f yields  T ∇f = 0T . . . 0T 0T . . . 0T 0T −1 .

(27)

Jacobian matrix of the constraints Ak is expressed as   M1 ... MN G 0 ... 0 N w1  I 0 0 −I 0 0 0 0   ..     .    0 0 I 0 0 −I 0 0   −2(u11 )T 0 0 0 0 0 0 0 T   . Ak =  ..  .    0 0 −2(uN G 1 )T 0 0 0 0 0    T 0 0 0 −2(u12 ) 0 0 0 0     . .   . T 0 0 0 0 0 −2(uN G 2 ) 0 0 (28) 11

Here, I are 5×5 unit matrices. It is worthy to note, that the IPOPT enforces the original optimization problems to have the form (20), and this requires the inequality constraints in (19) to be written as equality constraints with the additional slack variables. When the slack variables associated with a certain constraint are considered, the corresponding row in the Jacobian matrix has to be extended by adding columns with component 1. In practice, IPOPT does not require the user to include derivatives relating to the slack variables in the matrix. Instead, it requires the Jacobian matrices for both equality and inequality constraints to be provided separately, and IPOPT introduces the slack variables and extends the matrix (28) in the background. In order to assemble the Hessian of the Lagrangian function, ∇2xx L(xk , λk , zk ), IPOPT asks the Hessian matrix for each function separately, and it multiplies these matrices with corresponding dual variables and sum them into an integrity. It is noticeable from the optimization problem (19), that with the exception of the inequality constraints, all the functions in this problem are linear and their Hessians are zero matrices. In the case of inequality constraints, taking advantage of the reformulation, their Hessians are   2  2    2   . 2 ∇ c(uir ) =  (29)    2 2 It is immediately apparent that by multiplying the element matrices in (29) with corresponding lagrange multiplier and then assembling the products into a global Hessian, the resulting matrix will become a diagonal matrix. 2.5. Statistical models for result interpretation Weibull distribution introduced by Frecht and Weibull [84] has frequently been applied to describe the fatigue strength of material samples. The probability density function of a Weibull random variable follows: (
k x k−1 −(x/λ)k e x ≥ 0, f (x; λ, k) = λ λ (30) 0 x < 0, Here, k > 0 is the shape parameter, and λ is the scale parameter of the distribution. The cumulative distribution function for the Weibull distribution is expressed as: k F (x; k, λ) = 1 − e−(x/λ) (31) 12

for x ≥ 0, and F (x; k, λ) = 0 for x < 0. In the present work, after shakedown limits are calculated for each RVE sample, the results are fitted to (30) using the least square method to obtain k and λ. The consistency between the observations and the fitted Weibull distribution is checked using the χ2 goodness of fit test with following null and alternative hypotheses H0 : The data are consistent with a specified distribution Ha : The data are not consistent with a specified distribution Here, the specified distribution indicates the Weibull distribution with regressed k and λ values. Besides the χ2 goodness of fit test, there are two additional hypothesis tests, namely the two-sample Kolmogorov-Smirnov test (K-S test) and Wilcoxon rank sum test (rank sum test), applied in the present study to interpret the results obtained by shakedown analysis. The purpose of using these hypothesis tests is to examine if the data retrieved from two orthogonal directions, where θ = 0 and θ = 1/2π, are consistent. Because the material to be studied is globally isotropic given that the RVE size is sufficient, therefore the inconsistency between the data obtained from two directions can be used as an indicator of the insufficiency of the RVE size. K-S test examines if two samples X and Y are from the same continuous distribution. Null and alternative hypotheses of this test are H0 : Two samples are from the same continuous distribution Ha : Two samples are from the distinctive continuous distribution The decision of a two samples K-S test is made based on the distance between the empirical distribution functions of two samples, where the empirical distribution function indicates the culmulative distribution function of a sample that jumps up by 1/n at each of the n data points. The rank sum test, on the other hand, can be seen as a nonparametric equivalent to t-test which does not require the data to be subjected to the normal distribution. Null and alternative hypotheses of rank sum test are H0 : Two samples are from continuous distributions with equal medians Ha : Two samples are not from continuous distributions with equal medians As can be seen from the hypotheses statements, when K-S test and rank sum test are applied to the data in two directions, their results can be used to assess if the RVE size is sufficient and if it is justifiable to merge the data. 13

Table 1: Material properties of both phases

WC Co

E [GPa] ν [−] σY [MPa] 700 0.24 2000 210 0.30 683

3. Results and Discussion 3.1. RVE models for the statistical analysis Using the automatic mesh generator, we created 100 40µm-by-40µm models from different SEM images. These models are built by the generator using a uniform mesh configuration: the RVE models are meshed by linear wedge elements (C3D6). The elements covering non-critical regions were assigned with a global size of 0.8µm; near the phase boundaries a finer mesh with an edge size of 0.2µm (Fig. 3) was generated. Under this configuration, the number of elements for each RVE sample varies between 15,000-20,000. The material samples fall into two groups: the first group consists of 50 RVEs of WC with 30 Wt.% Co, and the second group of 50 models from the same material with 20 Wt.% Co (Fig. 4). Theoretically, when weight percents are converted to the volume percent, 20 Wt.% Co results in 30.5 Vol.% Co, while 30 Wt.% Co results in 43.9 Vol.% Co. However, according to our analyses of 100 SEM images, in average, samples from 20 Wt.% Co results in 37.5 Vol.% Co, while samples from 30 Wt.% Co in 45.8 Vol.% Co. The discrepancy between the theoretical value and our own data is assumed to be caused by the impurities of the binder phase: due to the production process of the material, the binder is a mixture of Co, W and C, thus the density is different compared to the pure Co. In all models, both constituents are considered as elastic-perfect plastic materials with the parameters given in Table 1. The study is restricted to small deformation and plastic failure is assumed to be the only failure mechanism. RVE models are built using ABAQUS [85] finite element software, which is commercially available. Using a Python script developed in-house, these models are assigned global loadings for calculating the retained elastic stress field. Once elastic stresses have been calculated, the geometrical setup of models and their associated results are transferred to MATLAB [86] to construct optimization problem (19). This problem is formulated as the transformation introduced in section 2.3 and then submitted to IPOPT for calculating the loading factor α. By altering the angle illustrated in Fig. 1, 14

Figure 3: Conversion of SEM images into 2.5D models

Figure 4: Materials with different binder contents

15

Limit Analysis

Σ22 [Mpa ]

Σ22 [Mpa ]

Shakedown Analysis

0 -1600

-800

0 Σ11 [MPa ]

800

1600

2.5D Coarse 2.5D Fine

-2000

-1500

P.Strain Coarse P.Strain Fine

-1000

-500 0 500 Σ11 [MPa ]

1000

1500

P.Stress Coarse P.Stress Fine

Figure 5: Mesh Insensitivity of the 2.5D RVE

the strength is evaluated in many different load combinations. Following the application of DM to each sample RVE, the results are collected, converted to macroscopic strength and finally regressed by Weibull distribution. Furthermore, the quality of the regression is assessed by the χ2 goodness of fit test. In addition to the distribution parameters λ and k, several statistical measurements including the mean value x¯, standard deviation s, as well as first, second and third quantile Q1 , Q2 , Q3 were likewise evaluated. Since both plane stress and plane strain models are extreme about the stiffness in the 3rd direction and, on the other hand, full sized 3D models are demanding in terms of both computational power and user expertise, a special type of thin-layered 3D model is deployed in this study. These models are generated by extracting elements of a 2D mesh for 1 µm along z axis and are intuitively referred to as 2.5D models. Compared to 2D models, other than being more realistic about the out-of-plane stiffness, 2.5D models are notable for their advantages in mesh insensitivity (c.f. Fig.5). It is important to be aware that the benefits accrued by adopting a 2.5D model are at the expense of a significant increase in the scale of the optimization problem. An arbitrary model is picked to illustrate this effect: the geometry consists of 17,739 integration points and 17,915 nodes in 2D. For the same model in 2.5D the number of integration points increases to 35,478 and the number of nodes to 35,830. When this model is used to study the shakedown limit under both uni-axial and bi-axial stresses, the increase in the size of the problem can be seen from Table 2. 16

Table 2: Scale of optimization problem for a FE mesh with different element types

Shakedown (Proportional) P.Stress P.Strain Num.Var 106,435 124,175 390,259 Num.EQ 89,047 89,047 284,880 Num.InEQ 35,478 35,478 70,956

Shakedown (Disproportional) 2.5D 745,039 639,660 141,912

3.2. Statistical Study of the Composite Strength The feasible load domains are shown in Fig. 6 and 7. In these figures as well as in the rest of the paper, the results are presented after normalization with respect to the binder strength in order to accentuate the contribution of the reinforcement particles. In these figures, the rigid lines represent the mean value and the shadowed regions the strength intervals [−2s, 2s]. The K-S test and rank sum test were applied to ΣU and Σ∞ calculated from two normal directions, and as shown in Table 3, for α fixed to 0.05 there is insufficient evidence to reject the H0 that the data belong to the same continuous distribution and the medians are the same. The results in Table 3 suggest that the RVEs meet the basic size requirement and it is justifiable to merge data in two directions to form one uniform strength statistic. The histogram and regressed Weibull distribution of ΣU and Σ2P ∞ are given in Figure 8 and 9, respectively. The endurance limit evaluated under two independently varying normal stresses Σ4P ∞ and the result is presented in 10. In addition to the graphical representations, the statistical parameters were calculated from the data and are summarized in Table 4. In this table, k and λ are Weibull parameters calculated from the data. In order to check the viability of describing the data via the Weibull distribution, χ2 goodness of fit test were performed. The results of the χ2 tests are collated in Table

Table 3: p values of hypothesis tests on RVEs in Groups 1 and 2

Wt.% Co 30 20

22 Σ11 U vs ΣU K-S Test Rank Sum 0.612 0.811 0.437 0.392

17

22 Σ11 ∞ vs Σ∞ K-S Test Rank Sum 0.789 0.784 0.553 0.437

Figure 6: Load domains pertained to WC - 30 Wt.% Co RVEs

Figure 7: Load domains pertained to WC - 20 Wt.% Co RVEs

18

100% 90%

WC-20 Wt.% Co WC-20 Wt.% Co WC-30 Wt.% Co WC-30 Wt.% Co

80%

FX (x) [-]

70% 60% 50% 40% 30% 20% 10% 0% 1.3

1.4

1.5

1.6

1.7

'U =
1.8

1.9

2

2.1

Figure 8: Culmulative distribution of ΣU and regressed Weibull distribution

100% 90%

WC-20 Wt.% Co WC-20 Wt.% Co WC-30 Wt.% Co WC-30 Wt.% Co

80%

FX (x) [-]

70% 60% 50% 40% 30% 20% 10% 0% 1.2

1.25

1.3

1.35

1.4

Co '2P 1 =
1.45

1.5

Figure 9: Culmulative distribution of Σ2P ∞ and regressed Weibull distribution

19

1.55

100% 90%

WC-20 Wt.% Co WC-20 Wt.% Co WC-30 Wt.% Co WC-30 Wt.% Co

80%

FX (x) [-]

70% 60% 50% 40% 30% 20% 10% 0% 0.7

0.75

0.8

0.85

0.9

Co '4P 1 =
0.95

1

1.05

Figure 10: Culmulative distribution of Σ4P ∞ and regressed Weibull distribution

5. According to this table, the p value of the tests when only the 5th to 95th percentile of the data, p0χ2 (0.05 − 0.95), are taken into account , is constantly greater than the predefined α = 0.05. Thus there was not sufficient evidence to reject H0 which states that the data is consistent with the Weibull distribution with k and λ given in Table 4. 4P When the statistics of ΣU , Σ2P ∞ and Σ∞ are compared to each other, it is immediately noticeable from Table 4 that in the case of RVEs with fixed geometries, s decrease as N V increases. However, when two RVE groups are compared, there are no noticeable dissimilarities between their λ. Concerning the statistics of ΣU , one noteworthy phenomenon is, that for both RVE groups there is some overlap between the limit and shakedown domains. This phenomenon is not unique and is frequently described in literature relating to the study of periodic composites, c.f. [87, 88]. When incremental FE simulation was used, it became clear that this phenomenon is attributable to the local effect: a band consisting of the weak binder phase perpendicular to the loading direction can be identified in the case of a few models. When those models are subjected to external load, severe localized plastic strain is exerted on this band. When this is the case, not the entire material but just one weak band on its own, determines the strength of the entire model.

20

Table 4: Statistics of strengths in RVE Groups 1 and 2

[-] ΣU /σYCo

Co Σ2P ∞ /σY

Co Σ4P ∞ /σY

Wt.% Co 30 30 20 20 30 30 20 20 30 20

θ 0 1/2π 0 1/2π 0 1/2π 0 1/2π 1/4π 1/4π

k 1.682 1.677 1.856 1.860 1.387 1.382 1.487 1.458 1.197 1.290

λ 15.1 14.9 15.8 15.4 15.7 16.3 14.3 16.4 18.4 18.5

x¯ 1.624 1.616 1.813 1.788 1.343 1.342 1.435 1.422 1.164 1.251

s 0.133 0.146 0.134 0.144 0.097 0.086 0.110 0.103 0.072 0.088

Q1 1.580 1.544 1.755 1.676 1.277 1.283 1.375 1.358 1.113 1.172

Q2 1.636 1.662 1.836 1.824 1.325 1.338 1.455 1.427 1.169 1.275

Q3 1.714 1.706 1.895 1.880 1.418 1.396 1.520 1.493 1.210 1.320

Table 5: Results of χ2 tests on strength data in RVE Groups 1 and 2

Wt.% Co 30 20

p χ2 0.020 0.038

ΣU 0 pχ2 (0.05 − 0.95) 0.168 0.292

p χ2 0.079 0.033

Σ2P Σ4P ∞ ∞ p0χ2 (0.05 − 0.95) pχ2 0.137 0.107 0.425 0.236

Presumably, for current materials when the model size increases, the chance diminishes that a band of this nature exists. Nevertheless, this phenomenon cannot be completely ruled out due to the randomness of the material. 4P For endurance limits Σ2P ∞ and Σ∞ , the scatter of the strengths, as schematically depicted as the transparent shadows in Fig. 6 and 7, has been found to relate to the failure modes. It is well known that failure can be caused by either incremental collapse or alternating plasticity, and the shakedown limit is equal to the minimum value of the two. In the case of RVEs considered in the present study, alternating plasticity is the major cause of the failure and in overall more than 90% RVEs fail due to this mechanism under varied uniaxial stress. Despite the importance of alternating plasticity, incremental collapse occasionally leads to failure of RVEs as well. It can be inferred on the basis of numerous case studies that in principle it is very hard, if not impossible, to distinguish between two failure modes from the elastic stress field without performing the direct method calculation. When RVEs fail due to incremental collapse, their Σ∞ will be much smaller than twice the elastic limit; and statistically it also results a greater s(Σ∞ ). The failure modes from 21

Fig. 6 and 7 can, therefore, be estimated roughly: more RVEs are expected to fail due to incremental collapse when a wider domain is observed. Finally, the influence of the binder contents can be summarized as follows: RVEs with 20 Wt.% Co have greater ΣU and Σ∞ compared to RVEs with 30 Wt.% Co, but their statistical characteristics show no clear differences. For both ΣU and Σ∞ , a large percentile of the data can be well described by Weibull distribution featured by (k, λ) calculated from the sample. The small ΣU is due to the local effect and strain concentration, while the reason for small Σ∞ is the incremental collapse. 4. Conclusion This paper shows how the lower bound direct method can be combined with statistical homogenization techniques to predict the strength of random heterogeneous materials. In the present study, a large amount of RVE samples were built and their global plastic behavior was studied via limit and shakedown analyses. Using an advanced nonlinear programing solver and linear algebra packages, large-scale nonlinear programmings originating from complicated real composite structures are formulated and solved. The success of solving DM problems of such scale and complexity proves that the method is a powerful tool for studying composite materials. By analyzing statistically the global behavior of various RVE samples under different types of loadings, the study reveals how the behaviors of PRMMCs are influenced by their composite structures. In addition, comparative study of samples with different binder contents revealed that most of the statistical indicators vary only slightly, which confirms the objectivity of the result. Finally, it is important to emphasize that the cause of the failure assumed in this paper is restricted to plastic dissipation. In our future studies additional material behaviors such as plastic hardening and damage caused by locally accumulated plasticity will be taken into account. To solve the mathematical programming associated with these more complicated problems, additional effort will be spent on improving the capability of the numerical method. References [1] N. Chawla, J.W. Jones, C. Andres, and J.E. Allison. Effect of SiC volume fraction and particle size on the fatigue resistance of a 2080 22

Al/SiCp composite. Metallurgical and Materials Transactions A-physical Metallurgy and Materials Science, 29(11):2843–2854, 1998. [2] N. Chawla, U. Habel, Y.-L. Shen, C. Andres, J.W. Jones, and J.E. Allison. The effect of matrix microstructure on the tensile and fatigue behavior of SiC particle-reinforced 2080 Al matrix composites. Metallurgical and Materials Transactions A-physical Metallurgy and Materials Science, 31(2):531–540, 2000. [3] S.-Y. Fu, X.-Q. Feng, B. Lauke, and Y.-W Mai. Effects of particle size, particle/matrix interface adhesion and particle loading on mechanical properties of particulate-polymer composites. Composites Part BEngineering, 39(6):933–961, 2008. [4] J. F¨ ussl and R. Lackner. Homogenization of strength: A numerical limit analysis approach. In J. Eberhardsteiner, C. Hellmich, H.A. Mang, and J. P´eriaux, editors, ECCOMAS Multidisciplinary Jubilee Symposium, volume 14 of Computational Methods in Applied Sciences, pages 183– 201. Springer Netherlands, 2009. [5] S.W. Tsai and E.M. Wu. A general theory of strength for anisotropic materials. Journal of Composite Materials, 5(1):58–80, 1971. [6] Z. Hashin. Failure criteria for unidirectional fiber composites. Journal of Applied Mechanics, 47(2):329–334, 1980. [7] R.M. Christensen. Tensor transformations and failure criteria for the analysis of fiber composite materials. Journal of Composite Materials, 22(9):874–897, 1988. [8] W.W. Feng. A failure criterion for composite materials. Journal of Composite Materials, 25(1):88–100, 1991. [9] P. De Buhan and A. Taliercio. A homogenization approach to the yield strength of composite materials. European Journal of Mechanics Asolids, 10(2):129–154, 1991. [10] A. Taliercio. Lower and upper bounds to the macroscopic strength domain of a fiber-reinforced composite material. International Journal of Plasticity, 8(6):741–762, 1992.

23

[11] A. Taliercio and P. Sagramoso. Uniaxial strength of polymeric-matrix fibrous composites predicted through a homogenization approach. International Journal of Solids and Structures, 32(14):2095–2123, 1995. [12] Z. Hashin and A. Rotem. A fatigue failure criterion for fiber reinforced materials. Journal of Composite Materials, 7(4):448–464, 1973. [13] K.L. Reifsnider and Z.J. Gao. A micromechanics model for composites under fatigue loading. International Journal of Fatigue, 13(2):149–156, 1991. [14] J. Degrieck and W. Van Paepegem. Fatigue damage modeling of fibrereinforced composite materials: Review. Applied Mechanics Reviews, 54(4):279–300, 2001. [15] R. Talreja. Damage and fatigue in composites-a personal account. Composites Science and Technology, 68(13):2585–2591, 2008. Directions in Damage and Durability of Composite Materials, with regular papers. [16] C.M. Friend. The effect of matrix properties on reinforcement in short alumina fibre-aluminium metal matrix composites. Journal of Materials Science, 22(8):3005–3010, 1987. [17] A. Borbe, H. Biermann, and O. Hartmann. FE investigation of the effect of particle distribution on the uniaxial stress strain behaviour of particulate reinforced metal-matrix composites. Materials Science and Engineering: A, 313:34–45, 2001. [18] S. Torquato. Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer, 2002. [19] B.A. Bednarcyk, S.M. Arnold, J. Aboudi, and M.-J. Pindera. Local field effects in titanium matrix composites subject to fiber-matrix debonding. International Journal of Plasticity, 20(8-9):1707–1737, 2004. [20] Y. Bansal and M.-J. Pindera. Finite-volume direct averaging micromechanics of heterogeneous materials with elastic-plastic phases. International Journal of Plasticity, 22(5):775–825, 2006. [21] J. Llorca, S. Suresh, and A. Needleman. An experimental and numerical study of cyclic deformation in metal-matrix composites. Metallurgical Transactions A, 23(3):919–934, 1992. 24

[22] O. Pierard, J. LLorca, J. Segurado, and I. Doghri. Micromechanics of particle-reinforced elasto-viscoplastic composites: Finite element simulations versus affine homogenization. International Journal of Plasticity, 23(6):1041–1060, 2007. [23] E. Melan. Zur Plastizit¨at des r¨aumlichen Kontinuums. Ingenieur Archiv, 9(2):116–126, 1938. [24] W.T. Koiter. A new general theorem on shake-down of elastic-plastic structures. In Proceeding of Eighth International Congress of Appllied Mechanics, pages 220–230, 1956. [25] S. W. Sloan. Lower bound limit analysis using finite elements and linear programming. International Journal for Numerical and Analytical Methods in Geomechanics, 12(1):61–77, 1988. [26] W.P. Shen. A general method for shakedown analysis. Computers & Structures, 30(4):901–903, 1988. Special Issue: Computational Engineering Mechanics. [27] E. Stein, G.B. Zhang, and J.A. K¨onig. Shakedown with nonlinear strainhardening including structural computation using finite element method. International Journal of Plasticity, 8(1):1–31, 1992. [28] M. Heitzer, G. Pop, and M. Staat. Basis reduction for the shakedown problem for bounded kinematic hardening material. Journal of Global Optimization, 17(1-4):185–200, 2000. [29] A. Hachemi, S. Mouhtamid, A.D. Nguyen, and D. Weichert. Application of shakedown analysis to large-scale problems with selective algorithm. In D. Weichert and A. Ponter, editors, Limit States of Materials and Structures, pages 289–305. Springer Netherlands, 2009. [30] J.-W. Simon, M. Kreimeier, and D. Weichert. A selective strategy for shakedown analysis of engineering structures. International Journal for Numerical Methods in Engineering, 94(11):985–1014, 2013. [31] Alan R.S. Ponter, Paolo Fuschi, and Markus Engelhardt. Limit analysis for a general class of yield conditions. European Journal of Mechanics A-solids, 19(3):401 – 421, 2000.

25

[32] H.F. Chen and A.R.S. Ponter. Shakedown and limit analyses for 3-D structures using the linear matching method. International Journal of Pressure Vessels and Piping, 78(6):443–451, 2001. [33] K.V. Spiliopoulos and K.D. Panagiotou. The residual stress decomposition method (RSDM): A novel direct method to predict cyclic elastoplastic states. In K.V. Spiliopoulos and D. Weichert, editors, Direct Methods for Limit States in Structures and Materials, pages 139–155. Springer Netherlands, 2014. [34] F. Pastor and E. Loute. Solving limit analysis problems: an interiorpoint method. Communications in Numerical Methods in Engineering, 21(11):631–642, 2005. [35] C.D. Bisbos, A. Makrodimopoulos, and P.M. Pardalos. Second-order cone programming approaches to static shakedown analysis in steel plasticity. Optimization Methods & Software, 20(1):25–52, 2005. [36] K. Krabbenhøft, A.V. Lyamin, S.W. Sloan, and P. Wriggers. An interiorpoint algorithm for elastoplasticity. International Journal for Numerical Methods in Engineering, 69(3):592–626, 2007. [37] Giovanni Garcea and Leonardo Leonetti. Numerical methods for the evaluation of the shakedown and limit loads. In 7th EUROMECH Solid Mechanics Conference, Lisbon, Portugal, 2009. [38] J.-W. Simon. Direct evaluation of the limit states of engineering structures exhibiting limited, nonlinear kinematical hardening. International Journal of Plasticity, 42(0):141–167, 2013. [39] D. Weichert, A. Hachemi, and F. Schwabe. Application of shakedown analysis to the plastic design of composites. Archive of Applied Mechanics, 69(9-10):623–633, 1999. [40] A. Hachemi, F. Schwabe, and D. Weichert. Failure investigation of fiberreinforced composite materials by shakedown analysis. In D. Weichert and G. Maier, editors, Inelastic Analysis of Structures under Variable Loads, volume 83 of Solid Mechanics and Its Applications, pages 107– 119. Springer Netherlands, 2000.

26

[41] H. Magoariec, S. Bourgeois, and O. D´ebordes. Elastic plastic shakedown of 3D periodic heterogeneous media: a direct numerical approach. International Journal of Plasticity, 20(8):1655–1675, 2004. [42] S. Bourgeois, H. Magoariec, and O. D´ebordes. A direct method for the determination of effective strength domains for periodic elastic-plastic media. In Weichert D. and Ponter A., editors, Limit States of Materials and Structures, pages 67–86. Springer Netherlands, 2009. [43] H. Zhang, Y.H. Liu, and B.Y. Xu. Plastic limit analysis of ductile composite structures from micro-to macro-mechanical analysis. Acta Mechanica Solida Sinica, 22(1):73–84, 2009. [44] J.-H. You, B.Y. Kim, and M. Miskiewicz. Shakedown analysis of fibre-reinforced copper matrix composites by direct and incremental approaches. Mechanics of Materials, 41(7):857–867, 2009. [45] M. Chen and A. Hachemi. Progress in plastic design of composites. In K. Spiliopoulos and D. Weichert, editors, Direct Methods for Limit States in Structures and Materials, pages 119–138. Springer Netherlands, 2014. [46] A.R.S. Ponter and F.A. Leckie. Bounding properties of metal-matrix composites subjected to cyclic thermal loading. Journal of The Mechanics and Physics of Solids, 46(4):697–717, 1998. [47] V. Carvelli. Shakedown analysis of unidirectional fiber reinforced metal matrix composites. Computational Materials Science, 31(1-2):24–32, 2004. [48] H.F. Chen and A.R.S. Ponter. On the behaviour of a particulate metal matrix composite subjected to cyclic temperature and constant stress. Computational Materials Science, 34(4):425–441, 2005. [49] H.X. Li and H.S. Yu. Limit analysis of composite materials based on an ellipsoid yield criterion. International Journal of Plasticity, 22(10):1962– 1987, 2006. [50] O. Barrera, A.C.F. Cocks, and A.R.S. Ponter. The linear matching method applied to composite materials: a micromechanical approach. Composites Science and Technology, 71(6):797–804, 2011. 27

[51] S. Torquato. Random heterogeneous media: Microstructure and improved bounds on effective properties. Applied Mechanics Reviews, 44(2):37–76, 1991. [52] W.J. Drugan and J.R. Willis. A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites. Journal of The Mechanics and Physics of Solids, 44(4):497–524, 1996. [53] A.A. Gusev. Representative volume element size for elastic composites: A numerical study. Journal of The Mechanics and Physics of Solids, 45(9):1449–1459, 1997. [54] T. Kanit, S. Forest, I. Galliet, V. Mounoury, and D. Jeulin. Determination of the size of the representative volume element for random composites: statistical and numerical approach. International Journal of Solids and Structures, 40(13):3647–3679, 2003. [55] M. Stroeven, H. Askes, and L.J. Sluys. Numerical determination of representative volumes for granular materials. Computer Methods in Applied Mechanics and Engineering, 193(30):3221–3238, 2004. [56] C. Pelissou, J. Baccou, Y. Monerie, and F. Perales. Determination of the size of the representative volume element for random quasi-brittle composites. International Journal of Solids and Structures, 46(14-15):2842– 2855, 2009. [57] A. Salahouelhadj and H. Haddadi. Estimation of the size of the rve for isotropic copper polycrystals by using elastic–plastic finite element homogenisation. Computational Materials Science, 48(3):447–455, 2010. [58] C. Heinrich, M. Aldridge, A.S. Wineman, J. Kieffer, A.M. Waas, and K. Shahwan. The influence of the representative volume element (RVE) size on the homogenized response of cured fiber composites. Modelling and Simulation in Materials Science and Engineering, 20(7):075007, 2012. [59] I.M. Gitman, H. Askes, and L.J. Sluys. Representative volume: Existence and size determination. Engineering Fracture Mechanics, 74(16):2518–2534, 2007. 28

[60] M.M. Kaminski. Computational mechanics of composite materials: sensitivity, randomness and multiscale behaviour. Springer Science & Business Media, 2006. [61] S. Swaminathan, S. Ghosh, and N.J. Pagano. Statistically equivalent representative volume elements for unidirectional composite microstructures: Part I-Without damage. Journal of Composite Materials, 40(7):583–604, 2006. [62] M. Ostoja-Starzewski. Material spatial randomness: From statistical to representative volume element. Probabilistic Engineering Mechanics, 21(2):112–132, 2006. [63] X.L. Yin, W. Chen, A. To, C. McVeigh, and W.K. Liu. Statistical volume element method for predicting microstructure - constitutive property relations. Computer Methods in Applied Mechanics and Engineering, 197(43-44):3516–3529, 2008. Stochastic Modeling of Multiscale and Multiphysics Problems. [64] T.J. Vaughan and C.T. McCarthy. A combined experimental-numerical approach for generating statistically equivalent fibre distributions for high strength laminated composite materials. Composites Science and Technology, 70(2):291–297, 2010. [65] M. Nyg˚ ards and P. Gudmundson. Three-dimensional periodic Voronoi grain models and micromechanical fe-simulations of a two-phase steel. Computational Materials Science, 24(4):513–519, 2002. [66] J. Mart´ın-Herrero and C. Germain. Microstructure reconstruction of fibrous C/C composites from X-ray microtomography. Carbon, 45(6):1242–1253, 2007. [67] R. Quey, P.R. Dawson, and F. Barbe. Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing. Computer Methods in Applied Mechanics and Engineering, 200(1720):1729–1745, 2011. [68] G. Chen, U.A. Ozden, A. Bezold, and C. Broeckmann. A statistics based numerical investigation on the prediction of elasto-plastic behavior of WC-Co hard metal. Computational Materials Science, 80(0):96–103, 2013. 29

[69] N. Karmarkar. A new polynomial-time algorithm for linear programming. Combinatorica, 4(4):373–395, 1984. [70] A.S. El-Bakry, R.A. Tapia, T. Tsuchiya, and Y. Zhang. On the formulation and theory of the Newton interior-point method for nonlinear programming. Journal of Optimization Theory and Applications, 89(3):507–541, 1996. [71] E.D. Andersen, C. Roos, and T. Terlaky. On implementing a primal-dual interior-point method for conic quadratic optimization. Mathematical Programming, 95(2):249–277, 2003. [72] M. Wright. The interior-point revolution in optimization: history, recent developments, and lasting consequences. Bulletin of The American Mathematical Society, 42(1):39–56, 2005. [73] J.-W. Simon and D. Weichert. Numerical lower bound shakedown analysis of engineering structures. Computer Methods in Applied Mechanics and Engineering, 200(41):2828–2839, 2011. [74] J.-W. Simon and D. Weichert. Interior-point method for lower bound shakedown analysis of von Mises-type materials. In G. de Saxc´e, A. Oueslati, E. Charkaluk, and J.-B. Tritsch, editors, Limit State of Materials and Structures, pages 103–128. Springer Netherlands, 2013. [75] F. Akoa, A. Hachemi, M. An, L.T.H. Said, and P.D. Tao. Application of lower bound direct method to engineering structures. Journal of Global Optimization, 37(4):609–630, 2007. [76] A.D. Nguyen, A. Hachemi, and D. Weichert. Application of the interiorpoint method to shakedown analysis of pavements. International Journal for Numerical Methods in Engineering, 75(4):414–439, 2008. [77] A. Hachemi, M. Chen, G. Chen, and D. Weichert. Limit state of structures made of heterogeneous materials. International Journal of Plasticity, 63(0):124–137, 2014. [78] J.A. K¨onig. Shakedown of elastic-plastic structures. Fundamental studies in engineering. Elsevier, 1987.

30

[79] A. W¨achter. An Interior Point Algorithm for Large-Scale Nonlinear Optimization with Applications in Process Engineering. PhD thesis, Carnegie Mellon University, 2002. [80] A. W¨achter and L.T. Biegler. On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming, 106(1):25–57, 2006. [81] O. Schenk. On fast factorization pivoting methods for sparse symmetric indefinite systems. Electronic Transactions On Numerical Analysis, 23:158–179, 2006. [82] O. Schenk and K. G¨artner. Solving unsymmetric sparse systems of linear equations with PARDISO. Future Generation Computer Systems-the International Journal of Grid Computing-theory Methods and Applications, 20(3):475–487, 2004. [83] R. Fletcher and S. Leyffer. Nonlinear programming without a penalty function. Mathematical Programming, 91(2):239–269, 2002. [84] W. Weibull. A statistical distribution function of wide applicability. Journal Of Applied Mechanics-Transactions of the ASME, 18(3):293– 297, 1951. [85] ABAQUS. ABAQUS/CAE user’s manual : version 6.13. Simulia, Dassault Syst´emes, 2013. [86] MATLAB. version 8.4.0 (R2014b). The MathWorks Inc., Natick, Massachusetts, 2014. [87] F. Schwabe. Einspieluntersuchungen von Verbundwerkstoffen mit periodischer Mikrostruktur. PhD thesis, RWTH Aachen, 2000. [88] M. Chen, A. Hachemi, and D. Weichert. Shakedown and optimization analysis of periodic composites. In G. de Saxc´e, A. Oueslati, E. Charkaluk, and J.-B. Tritsch, editors, Limit State of Materials and Structures, pages 45–69. Springer Netherlands, 2013.

31