Numerical limit and shakedown analysis method for kinematic hardening structure made of arbitrary inhomogeneous material

Composite Structures xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Numerical limit and shakedown analysis method for kinematic hardening structure made of arbitrary inhomogeneous material ⁎

Song Huanga, Hu Huia, , Zhiping Chenb a b

School of Mechanical and Power Engineering, East China University of Science and Technology, Shanghai 200237, PR China Institute of Process Equipment, Zhejiang University, 38 Zheda Road, Hangzhou, Zhejiang 310027, PR China

A R T I C LE I N FO

A B S T R A C T

Keywords: Inhomogeneity Kinematic hardening Limit load Shakedown Basis reduction method

The prediction of structure’s ultimate load has a lot of practical applications in composite structure. In the present work, limit and shakedown load of structure made from arbitrary inhomogeneous material was predicted by a novel solution algorithm named Back Stress Iteration Method, where material’s kinematic hardening behavior is involved as well. Firstly, details about the proposed algorithm were introduced. Basis reduction method was employed as framework and its advantage in computational efficiency is inherited. Secondly, the influence of inhomogeneous kinematic hardening was solved by a new invented Back Stress Iteration process. Then, three examples with increasing complexity were investigated numerically for validation purpose. The results were discussed intensively with respect to the accuracy, the effect of inhomogeneity and the inhomogeneous kinematic hardening behavior to verify the rationality of the proposed method. At last, the paper was summarized and conclusions were drawn. The proposed method is promising to provide a powerful tool for the accurate prediction of arbitrary inhomogeneous structure’s load-bearing capacity, such is applicable to many kinds of composite structures.

1. Introduction Load-bearing structures with inhomogeneous mechanical properties or made from composites are common scenarios in modern industries. Material strength in structure changes according to several mechanisms. For instances, material properties in a hydrogenation reactor shell changes through the thickness due to manufacturing residual influence [1]. For underwater structures such as submarine pipeline and ship structure, corrosion is the main reason of non-uniform material degeneration [2]. Servicing tailored structures made of function gradient material (FGM) are widely used in automobile industry [3] and aircraft industry [4]. For composites such as fiber-reinforced [5] or particulateenhanced material [6], the heterogeneity nature of the composite structures must be considered in strength analysis. For these structures, the prediction of ultimate load (limit load under monotonic load or elastic shakedown load under cyclic loads) requires a reasonable consideration about the impact of material inhomogeneity. Meanwhile, failure mode of steel structure subjected to complex loads depends strongly on material’s Bauschinger effect, i.e. kinematic hardening behavior. Therefore, a general ultimate load prediction method which involves in inhomogeneous material properties and kinematic hardening behavior simultaneously benefits the accuracy of industry

⁎

product design, serving performance prediction and life assessment, such can significantly promote the development of modern industry. In most engineering cases, detailed load history applied on structures components is usually difficult to be determined. Consequentially, conventional step-by-step analysis which is based on incremental finite element (FE) analysis is not applicable. In these cases, direct method based on Melan theorem [7] and Koiter theorem [8] provides alternative solution which predicts structure’s ultimate load with only the knowledge of load domain. For limit and shakedown analysis made of homogeneous kinematic hardening material, Weichert & Groβ-weege [9] firstly presented the famous Two Surface Model involving limited kinematic hardening into Melan shakedown theorem. Since then, numerous techniques were developed to solve Two Surface Model. Stein et al. [10] proposed an overlay model, which decomposes the kinematic hardening shakedown problem into several overlapped elasto-perfectly plastic shakedown problems. Heitzer et al. [11] developed a technique which constructs back stress of the kinematic hardening shakedown problem with the solution of a pre-solved elasto-perfectly plastic problem. Huang et al. [12] further develop Hetizer et al.’s method to complex multi-axial load condition. Simon et al. [13] invented a high efficient solving algorithm based on problem-tailored internal point algorithm, which can solve the nonlinear mathematic programming

Corresponding author. E-mail address: [email protected] (H. Hui).

https://doi.org/10.1016/j.compstruct.2019.111641 Received 5 May 2019; Received in revised form 10 September 2019; Accepted 30 October 2019 0263-8223/ © 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: Song Huang, Hu Hui and Zhiping Chen, Composite Structures, https://doi.org/10.1016/j.compstruct.2019.111641

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

(NLP) of Two Surface Model directly. These contributions provide solving techniques with acceptable computational consume, but they cannot be applied to the structure made of inhomogeneous material. Homogenization idea is usually employed to connect the inhomogeneous micromechanical with macro mechanical properties. After homogenization, the mechanical behavior represented by Representative Volume Element (RVE) can be applied on macroscopical analysis methods. The contributions with respect to homogenization method and shakedown analysis can be found in literatures [14–16]. However, the feasibility of homogenization method is confined to the material whose inhomogeneity is periodical. In the case of randomly distributed or non-periodical inhomogeneity, an alternative solution is to set up numerous RVEs from the material and use statistical based method [17]. If the inhomogeneity changes along one dimension of the structure, overlay model (OLM) along the concerned dimension can be employed, representative case is carbon fiber composite structures [18,19]. If the inhomogeneity is 2 or 3-dimensional, the problem may be solved by sophisticated software reengineering techniques. Nevertheless, investigations concerning shakedown problem with arbitrary inhomogeneity remains scarce and relevant numerical technique is still in sorely needed. The aim of this paper is to propose a novel solving algorithm for limit and shakedown analysis involving arbitrary material inhomogeneity and kinematic hardening behavior simultaneously. The present work is organized as following: Section 2 introduces theoretical fundamentals of inhomogeneous kinematic hardening shakedown analysis. In Section 3, details of the proposed method are illustrated elaborately. In Section 4, three examples with 1-dimensional, 2-dimensional and randomly inhomogeneity are calculated by the proposed method. In Section 5, the results are compared and discussed for validation purpose. At last, conclusions are drawn in Section 6. The proposed method is promising to provide a powerful tool for the accurate prediction of structure’s load-bearing capacity, whose potential application scenario includes most of composite structures.

Fig. 1. Illustration of load corners in load domain.

max α s.t. F [ασiE (j ) + ρi − πi] ⩽ σs Gs (xi ) F [πi] ⩽ σu Gu (xi ) − σs Gs (xi ) N

∑i =G1 Ci·ρi = 0 i = 1, 2, …NG, j = 1, 2, 3, …NL

where xi represents the coordinate in structure, Gs(x) and Gu(x) are the distribution function of yield and ultimate strength, which depends on a pre-determined mechanism of strength inhomogeneous. Problem (2) represents a group of large-scaled NLP, who’s solution technique is the main challenge for engineering application and the interest of present work. 3. Back stress iteration method In a three-dimensional discrete system, problem (2) is an NLP with 12NG + 1 independent variables and (NL + 1) × NG inequality constraints. Both the number of independent variables and constraints will be reduced by the proposed method.

2. Theoretical fundamentals Two Surfaces Model proposed by Weichert & Groβ-weege [9] provides the necessary condition of shakedown for kinematic hardening structure. Later on, modified two surfaces model was further developed by Simon [20] and Kanno [21], which is the necessary and sufficient condition of kinematic hardening shakedown. It is stated as

3.1. Framework based on basis reduction method Basis reduction method employed by Stein et al. [10], Liu et al. [22] and Chen et al. [23] is adopt as the framework of the proposed solution technique. The idea is to solve a series of sub-problem with small scale iteratively to approximate the solution of its large-scaled original problem. In basis reduction method, residual stress field ρ in Melan theorem is stated as

max α s.t. F [ασiE (j ) + ρi − πi] ⩽ σs, i F [πi] ⩽ σu, i − σs, i

ρ = Bm ·Xm

N

∑i =G1 Ci·ρi = 0 i = 1, 2, …NG, j = 1, 2, 3, …NL

(2)

(3)

where Bm = [b1, b2,…bm] is the matrix of self-equilibrated stress basis vector, bm is a basis vector of the stress space in m dimensions, Xm = [x1, x2,…xm]T denotes the coefficient matrix of Bm. Bm describes a self-equilibrated stress space and therefore ρ = Bm∙Xm satisfies ∇∙ρ = 0 automatically. Then, a series of sub-problems are considered to replace problem (2). The sub-problem constructs residual stress ρ in a r dimensional sub-space of Bm(r ≪ m), i.e. ρ = Br∙Xr, which leads to a subproblem

(1)

where α represents load multiplier, ρ is the residual stress field which is self-equilibrate, F(∙) is yield function in Mises sense, i ∈ NG denotes i-th Gaussian point in the discrete system, j ∈ NL denotes the j-th load corner of the applied load domain, as is shown in Fig. 1, Ci is the convert matrix depending on the discrete system, σiE (j ) represents the virtual elastic stress at Gaussian point i under the function of load corner j, πi is the back stress, σs,i and σu,i denote the yield and ultimate strength respectively. For engineering structures made of steel with sufficient ductility, the failure is controlled by material’s strength. Accordingly, the influence of material inhomogeneity is represented by non-uniform yield and ultimate strengths in the structure. Taken inhomogeneous material strengths into problem (1), one can derive

max α k s.t. F [α kσiE (j ) + Br , i ·Xr − πi] ⩽ σs Gs (xi ) F [πi] ⩽ σu Gu (xi ) − σs Gs (xi ) i = 1, 2, …NG, j = 1, 2, 3, …NL

(4)

where k represents the k-th sub-problem. Sub-problem (4) has r + 6NG + 1 independent variables and NL × NG inequality constraints. By solving sub-problem (4) iteratively, the object variable αk searches for the optimal ascending direction in a series of sub-spaces and 2

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

approaching the global optimum of original problem (2) gradually. The convergence of iteration is achieved when

|α k − α k − 1| ⩽ εerr α k−1

(5) −6

where εerr is the tolerance of convergence, which is set as εerr = 10 in the present work. The key of basis reduction method is the construction of basis vector matrix Bm. In the present work, technique employed by Chen et al. [23] is adopted. For a specific load vector applied on the structure, the stress state after k-1 th sub-problem is represented as σk − 1 = αk − 1σE + Brk − 1Xrk − 1, for k th sub-problem, the basis vector matrix Brk can be constructed as following: Applying a load increment ΔαP, the structure will yield further and the resulting elasto-plastic stress updating problem can be solved by an mN-R iteration method. The equilibrium equation in the first mN-R iteration is NG

(α k − 1 + Δα ) P −

NG

∑ Ci·σik −1 = ∑ Ci·D·Δεi1 i=1

i = 1, 2, 3, …NG

i=1

(6)

where P denotes the reference load vector, D is the elastic matrix of the material, Δε represents the increment of strain. In the q th mN-R iteration, Eq. (6) yields NG

(α k − 1 + Δα ) P −

NG

∑ Ci·σik −1,q−1 = ∑ Ci·D·Δεiq i=1

i = 1, 2, 3, …NG

i=1

Fig. 2. Relationship between total stress and back stress on failure point.

(7)

Subtracting Eq. (6) from Eq. (7), one have

components simultaneously, back stress have to be collinear with total stress in the limit state, as is shown in Fig. 2(b) & (d). Considering the iterative nature inheriting from basis reduction method and relationship Eq. (10), BSIM suggests a back stress field for k th sub-problem (4) as

NG

∑ Ci [σik −1,q−1 − σik −1 + D (Δεiq − Δεi1)] = 0

i = 1, 2, 3, …NG

i=1

σik − 1, q − 1

σik − 1

D (Δεiq

(8)

Δεi1)

− + − is a basis Eq. (8) indicates bq − 1 = vector. During this procedure, q times N-R iteration is required for constructing a matrix Brk with q−1 basis vectors. Then, for all active load corners j∗ ∈ NLA (active load leads to elasto-plastic response of the structure, NLA represents the collection of all active load corners), basis vectors will be generated, Brk will consist of the basis vectors from all k k , b2,1 , …brk− 1, j∗, brk− 1, j∗ brk, j∗]. j∗ ∈ NLA , i.e. Brk = [b1,1

πik = kKH , i Πik

i = 1, 2, …NG

(11)

where superscript k denotes the iteration time of basis reduction method, kKH,i is the compensation factor, Πk is effective stress field depending on the current stress state. Based on the solution of k−1 th sub-problem σk−1, i.e. αk−1σE + ρk−1, Πk is constructed according to NL

3.2. Back stress iteration technique

Πk =

∑ δ (j)[α k −1σ E (j) + Brk −1·Xrk −1]

j = 1, 2, …NL (12)

j=1

In sub-problem (4), back stress π leads to 6NG independent variables and NG inequality constraints, which increases additional computation burden. Further simplification of sub-problem (4) can be achieved by a novel Back Stress Iteration Method (BSIM) proposed in the present work. The idea of BSIM came from following observation: in order to maximize the external loads applied on the structure, back stress at failure point of the structure have to be collinear with its total stress in limit state. This phenomenon is illustrated in Fig. 2. and proved as following. According to Two surfaces model, the constraints of stress components at a given Gaussian point are

F [ασiE + ρi ] ⩽ σu

where δ(j) is stated as

1 j ∈ NLA δ (j ) = ⎧ ⎨ 0 j ∈ NL − NLA ⎩

In Eq. (11), compensation factor kKH,i at Gaussian i is defined as

kKH , i = min[kex1 (1), kex1 (2), …kex1 (j ), kex 2] j ∈ NLA

rmax = max( ∥Πik ∥e ) kex 2 =

σu Gu (xi) − σs Gs (xi ) rmax

(15)

where ||σ||e = F[σ] represents the norm of stress vector σ in Von-Mises stress space. Eq. (11) implies the action that generate back stress field which is collinear with current effective stress field. The scale of transformation kex2 guarantees the maximal stress in the structure stays in ultimate strength boundary. Then, following inequality holds

(9)

These relationships are shown in Fig. 2(a) & (c), where the direction of virtual elastic stress σiE is fixed. Load multiplier α, residual stress ρi and back stress πi are adjustable. In limit state, Mises stress equal to ultimate stress σu while α reaches its maximum, i.e. F [ασiE + ρi ] = σu . Meanwhile, Eq. (10) holds according to the triangle law between stress vectors.

F [ασiE + ρi ] = F [ασiE + ρi − πi] + F [πi]

(14)

kex2 is

for total stress

F [ασiE + ρi − πi] ⩽ σs for frictional stress F [πi] ⩽ σu − σs for back stress

(13)

F [πik ] = kKH , i F [Πik ] ⩽

σu Gu (x i ) − σs Gs (x i ) ∥Πik ∥e ⩽ σu Gu (x i ) − σs Gs (x i ) max( ∥Πik ∥e ) (16)

As a result, the proposed back stress satisfies the constraints of back stress in Eq. (9). In Eq. (11), kex1(j) is defined as the positive root of following

(10)

Eq. (10) implies that in order to satisfy the constraints of stress 3

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 3. Restriction of back stress movement.

equation

∥Πik ∥2e kex2 1 − 2σik − 1 (j )·Πik kex1 + ∥σik − 1 (j ) ∥2e − [σs Gs (x i )]2 = 0

Fig. 4. Stress evolution during BSIM iteration (NL = 1).

(17) It’s nature to notice that process shown in Fig. 4 suffers numerical unstable when the structure is subjected to a set of complex load combinations (NL > 1) because the stress field generated by different load vector may lead to a unfeasible back stress. As is shown in Fig. 5(a), both stress field σk−1(1) and σk−1(2) generated by load vector 1 and load vector 2 affect the movement of SYS, in which case the ignorance of either load stress field will lead to fatal infeasibility of the algorithm (stress yields beyond SYS). In BSIM, this obstacle is overcome by the action of effective stress field Π. The effective stress Πk is the vector summation of stress fields which contribute to plastic flow, therefore, the back stress constructed based on Πk guarantees that all mass points stay in the SYS under any applied load vector so that yield condition will not be violated in any case, as is shown in Fig. 5(b), such guarantee the robustness of BSIM.

Derivation of Eq. (17) is illustrated in Fig. 3. The shifted yield surface (SYS) transforms in stress space along the direction of Πk and kKH,i determines the movement of SYS. Fig. 3 demonstrates three possible SYSs named as SYS1, SYS2 and SYS3 as results of three possible back stress fields kex1Πk. It’s noticed that SYS3 is infeasible as current stress field σk−1(j) moves out of SYS3 in k th sub-problem. Therefore, maximal back stress should not violate the critical condition represented by SYS2, i.e. σk−1(j) lies on SYS. Accordingly, one can derive that in the critical state (SYS2), back stress and current stress satisfies

[kex1 ∥Πik ∥e ]2 + ∥σik − 1 (j ) ∥2e − 2 cos θkex1 ∥Πik ∥e ∥σik − 1 (j ) ∥e = [σs Gs (x i )]2 (18) With the definition of Π , θ can be calculated by k

cos θ =

Πik ·σik − 1 (j ) ∥Πik ∥e ∥σik − 1 (j ) ∥e

(19) 3.3. Numerical implementation

Substituting Eq. (19) into Eq. (18), one can obtain Eq. (17). By using Eq. (11), sub-problem (4) is transformed into

max α k s.t. F [α kσiE (j ) + Bri, i ·Xrk − πik ] ⩽ σs Gs (x i ) i = 1, 2, 3, …NG j = 1, 2, …NL

The proposed method was implemented on CODE MATLAB, while CODE ANSYS is employed for the pre-process and post-process. The flow chart of the proposed method is shown in Fig. 6. (20) 4. Numerical examples

Sub-problem (20) consists of only r + 1 independent variables and NL × NG constraints, leading to an NLP whose scale is much reduced compared with sub-problem (4) as well as original problem (2). Noteworthy, in order to construct basis vector matrix with Eqs .6–8 which is originally designed for elasto-perfectly plastic problem, one must generate basis vector based on the frictional stress component, i.e. σ-π, as

σ k − 1 (j ) = α k − 1σ E (j ) + Brk − 1·Xrk − 1 − π k − 1

j = 1, 2, …NL

Three typical numerical examples were investigated for validation purpose. These examples are with increasing dimensional inhomogeneous. For each instance, the discussion is organized as following: Firstly, conventional homogeneous kinematic hardening shakedown analysis is performed. The results are compared with available results in literatures. This would rationalize the numerical accuracy of the established program and provide precondition for further discussions. Then, same problem is evaluated under the influence of pre-defined inhomogeneity and compared with the results from step-by-step analysis. Finally, the impact of inhomogeneity and the rationality of the proposed method are discussed. The step-by-step analysis mentioned here represents the conventional elasto-plastic incremental FE analysis conducted on Commercial software ANSYS. The material behavior is considered to be multi-linear isotropic/kinematic hardening combined with J2 flow theory. In all examples, material properties in Table.1 will be used.

(21)

In order to explain the mechanism of BSIM, limit load case (NL = 1) is considered and illustrated in Fig. 4, where UYS, IYS and SYS denote Ultimate Yield Surface (UYS), Initial Yield Surface (IYS) and Shifted Yield Surface (SYS) respectively. In the iteration process of BSIM, SYS shifts according to the current effective direction defined by Πk. Consequently, back stress always tries to catch up total stress, as is shown in k = 1, 2, 3. When the structure reaches limit state (k = 4), another iteration cannot provide a more acceptable stress field compared with its current state. As a result, both total stress and back stress will cease to evolve during iterating, which triggers convergence criterion Eq. (5) and stops BSIM algorithm. 4

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 5. Back stress movement during BSIM iteration (NL > 1). Table 1 Material properties in numerical examples. Property

Value

elastic module E Poisson’s ratio ν Thermal expansion γ

200 GPa 0.3 1.5 × 10−4 K−1

instance is considered to be in radius direction (1-dimensional), which can be taken as a pipe made from FGM. For the contributions about shakedown load of this structure, one can refer literatures [11,20]. 4.1.1. Model setup Model description is shown in Fig. 7. The inner pressure p and inner temperature T1 can change independently while outside temperature T0 is equal to consistent environment temperature. The established 1/4 Finite Element (FE) model is shown in Fig. 8, which is modeled with ANSYS’s 3D brick element SOLID185. Normal displacement is restrict to 0 on the transection in x, y and z direction. The model consists of 384 elements and 612 nodes, resulting in an NLP with 36,865 independent variables and 9216 constraints if raw problem (2) is considered. With the application of the proposed method (r = 4), the scale of this problem is reduced significantly as 5 independent variables and 6144 constraints. 4.1.2. Homogeneous case Homogeneous material with elasto-perfectly plastic(EP) and kinematic hardening(KH) σu = 1.35σs, σu = 1.5σs, σu = 2σs are considered. The resulting shakedown load domain is shown in Fig. 9. The results are normalized by the shakedown load of EP case under pure inner pressure p0 and pure thermal load ΔT0. 4.1.3. Inhomogeneous case In inhomogeneous case, strength described by Eqn.(22) is applied to the pipe. Strength contour of the pipe is shown in Fig. 10. The nonuniform strength changes through the radius direction. x 2 100

3 y 2 100

( ) + ( ) ⎤⎦ σ G (x , y, z ) = 300 + 150 ⎡ ( ) + ( ) ⎤ ⎣ ⎦ σs Gs (x , y, z ) = 200 + 100 ⎡ ⎣ u

u

x 2 100

y 100

2 3

MPa MPa

(22)

Shakedown load of the inhomogeneous pipe was predicted by the proposed method. Overlay Model (OLM) based on step-by-step analysis was also employed. OLM pipe is a cylinder consisting of 10 layers stacking on each other. For each layer, homogeneous behavior is considered, as is shown in Fig. 11. The material strengths in layers are listed in Table 2. Specific load history is required for OLM analysis. The applied load spectrum is shown in Fig. 12. A total of 120 load steps were calculated.

Fig. 6. Flow chart of BSIM algorithm.

4.1. Example I: thin-walled pipe under thermal-mechanical loads The first example is a pipe subjected to variable inner pressure and thermal load, i.e. the famous Bree problem. Inhomogeneity in this

5

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 7. Dimensions and loads of the pipe under thermal-mechanical load.

4.2. Example II: holed plate under biaxial tensile loads Classical holed plate problem is chosen as the second example. The material inhomogeneity is set as 2-dimensional this model (x and y direction). Limit and shakedown load of this structure can be found in Refs. [10,11,20,23]. 4.2.1. Model setup As is shown in Fig. 15, the plate is subjected to variable tensile pressure P1 and P2 in x and y direction. Both P1 and P2 can change independently. A 1/4 sub-model of the plate was considered. Normal displacement is restricted to 0 on x & y transections. Besides, uz = 0 is assigned to a random chosen node to restrict rigid body deformation. Fig. 16 shows the FE model of the plate, which has 864 elements and 1425 nodes, resulting in an NLP with 82,945 independent variables and 34,560 constraints for problem (2) or an NLP with 5 independent variables and 27,648 constraints (r = 4) for the proposed method.

Fig. 8. Finite element model of the pipe under thermal-mechanical load.

4.2.2. Homogeneous case In homogeneous case, both EP behavior σu = σs and KH behavior σu = 1.5σs are considered. Limit and shakedown loads predicted by the proposed method are shown in Fig. 17, the Mises stress contours under different load cases are shown in Fig. 18. 4.2.3. Inhomogeneous case Inhomogeneity is introduced into the structure through a pre-defined distribution function, as is illustrated in Eq. (23). Contour of material strengths are shown in Fig. 19.

σs Gs (x , y ) = 200 + σu Gu (x , y ) = 300 +

x2 500

+

x2 250

y2 125

+

y2 250

MPa MPa

(23)

Limit load of the non-uniform plate was calculated by the proposed method. For comparison purpose, Field Variable Method (FVM) was employed. FVM is a step-by-step elasto-plastic analysis method where material properties depends on field variables. In the present work, FVM was implemented on CODE ABAQUS with user-subroutine USDFLD and field variables were chosen as the coordinates of mass points. It should be pointed out that FVM cannot simulate kinematic hardening behavior. Alternatively, isotropic hardening (IH) behavior was chosen to calculate the limit load because the difference between KH and IH disappears under monotonous load. Results of the proposed method and FVM are listed in Table 4.

Fig. 9. Shakedown domain of the homogenous pipe.

Shakedown domain resulting from the proposed method is shown in Fig. 13. The mechanical response of the pipe under some specific load combinations (validation points) were calculated by OLM. The validation points are also shown in Fig. 13 and the material behaviors considered are listed in Table.3. The resulting maximal equivalent plastic strain at validation points are shown in Fig. 14. These validation points are expected to be located at both sides of the shakedown boundary, therefore, the validity of shakedown boundary predicted by the proposed method can be evaluated based on the strain response at these points.

4.3. Example III: smooth tensile sample under uniaxial load In the last example, the potential of the proposed method is demonstrated by a smooth tensile sample subjected to uniaxial tensile load. In this example, the material strength changed randomly at each Gaussian point, corresponding to a particulate reinforced metal matrix composite [17]. Limit load of the smooth sample was investigated. The inhomogeneous case was compared with the results from 6

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 10. Strength distribution in the inhomogeneous pipe.

homogenization method.

Table 2 Material strength in overlay model.

4.3.1. Model setup The sketch of smooth tensile sample is illustrated in Fig. 22. 1/4 FE model and the corresponding boundary conditions of the model are shown in Fig. 23. The model consists of 2623 brick elements and 3432 nodes. 4.3.2. Material properties description In homogeneous case, yield and ultimate strength of the material are set as 200 MPa and 300 MPa respectively. When material inhomogeneity is considered, yield and ultimate strength of the sample are described as

σs Gs (x , y, z ) = 150 + 50 × gs (x , y, z ) MPa σu Gu (x , y, z ) = σs Gs (x , y, z ) + 100 × gu (x , y, z ) MPa

Number of layer

Yield strength (MPa)

Ultimate strength (MPa)

1 2 3 4 5 6 7 8 9 10

255 259 263 267 271 276 281 286 291 297

382 388 394 400 407 414 421 429 437 446

homogenized material strengths at a given transection are determined according to

(24)

where gs(x, y, z) and gu(x, y, z) are stated as

random number in (0, 1) for x ⩽ 100 gs (u) (x , y, z ) = ⎧ 1 for other points ⎨ ⎩

20

20

1

20

20

⎧ Ngps ∑y = 0 ∑z = 0 σs Gs (x , y, z ) x < 100 ⎨ 200 x > 100 ⎩

σ¯u (x ) =

⎧ Ngps ∑y = 0 ∑z = 0 σu Gu (x , y, z ) x < 100 ⎨300 x > 100 ⎩

(25)

Strength distribution in the sample are shown in Fig. 24. It is demonstrated that the lower boundary of yield and ultimate strength are 150.017 MPa and 150.076 Mpa respectively, whereas the upper boundary of yield and ultimate strength are 200 MPa and 300 MPa. It’s also noticed that yield strength and ultimate strength changes according to different patterns. For comparison purpose, homogenization concept is applied. The

1

σ¯s (x ) =

(26)

where Ngps denotes the number of Gaussian points in transection. Based on Eqn.(25), the overall material strengths of the homogenized sample are

Fig. 11. Overlay model of the pipe. 7

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 12. Load spectrum for overlay model.

5. Discussions 5.1. Numerical accuracy in homogeneous case It’s demonstrated in Fig. 9 that for example I, the shakedown load domain predicted by the proposed method coincides well with those proposed by Heitzer et al. [11] and Simon [20]. When thermal load controls the failure process of the pipe, the predicted result shows insignificant difference between EP case and KH cases. On the other hand, KH cases exhibit an enlargement of the shakedown load as a result of kinematic hardening when mechanical load (inner pressure) is the primary load. The predicted result indicates kinematic hardening behavior enhances the resistance against ratchet failure of the pipe while thermal fatigue resistance cannot be improved. Same conclusion can be found in literatures [11,20,24]. In example II shown in Fig. 17, it is obvious that the limit and shakedown load domain resulting from the proposed method agree well with the results from Chen et al. [23] and Akoa et al. [24] under both EP and KH cases. The predicted result exhibits an enlarged limit load boundary and an unchanged shakedown load domain after involving kinematic hardening behavior, which is also reported in Refs. [10,11,20,23]. In example III, results of the proposed method match well with the theoretical solution in homogeneous case, as is illustrated in Table.5. As a result, it’s concluded that proposed method is of sufficient accuracy when inhomogeneous problem degenerates into homogeneous case. Given that the homogeneous case is a special case of the proposed method, therefore, the numerical accuracy in homogeneous cases provides sound fundamental for further discussion about the effect of inhomogeneity.

Fig. 13. Shakedown domain of the inhomogeneous pipe. Table 3 Material behaviors of the validation points. Elasto-perfectly plastic

Kinematic hardening

A/B/E/F/I/K/M/O

C/D/G/H/J/L/N/P

σ¯s = min[σ¯s (x )] σ¯u = min[σ¯u (x )]

(27)

Eqs. (26) and (27) indicate that the homogenized material properties depending on the weakest transection in the sample which failures firstly.

5.2. Effect of inhomogeneity Plastic collapse and elastic shakedown are the main failure modes for engineering structure made from steel with high ductility, where failure is controlled by material strength. In these failure modes, a mass point in structure sustains external loads until equivalent stress reaching yield or ultimate strength. After that, the mass point suffering unconfined plastic flow or progressive ductile damage is considered failure. In example II, by comparing following contour pairs: Fig. 18(a)/ Fig. 20 (a), Fig. 18(b)/Fig. 20(c), Fig. 18(c)/Fig. 21 (a), Fig. 18(d)/ Fig. 21(c), one can find stress redistribution in inhomogeneous plate resulting from strength inhomogeneity. According to Figs. 20 and 21

4.3.3. Limit load According to the model shown in Fig. 22, the theoretical solution of the sample is given as

PL =

σs (u) MPa 4

(28)

Limit load of the sample in both homogenous and inhomogeneous cases were calculated by the proposed method and Eq. (28). The results are listed in Table.5. The stress fields in limit state are shown in Fig. 25. 8

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 14. Equivalent plastic strain evolution for each validation points.

the structure, the load-bearing capacity variation resulting from the strength inhomogeneity altered the limit load of the sample, as is illustrated in Table.5. Therefore, it’s concluded that the proposed method has reproduced the effect of inhomogeneity on the load-bearing capacity reasonably.

and the results in Table.4, the results predicted by the proposed method match well with FVM results in all cases. For EP case, limit load P1 is higher than P2, indicating the inhomogeneous plate has higher strength in x direction. For SH case, the difference between P1 and P2 is neglectable, implying similar plate strength in both x and y direction. Causality between the anisotropic load-bearing capacity and material inhomogeneity can be explained by the limit stress contours. In EP case, it is exhibited that limit stress state shown in Fig. 20(a) has a higher stress level compared with that shown in Fig. 20(c), indicating higher inner force under x directional tensile which can sustain more external loads. Similar to EP case, in SH case, the inner force level under P1 and P2 are close because the ultimate strength distribution shown in Fig. 19(b) allows a sufficient plastic deformation and stress redistribution, therefore, difference between limit loads in x and y direction disappeared in SH case. In example III, Fig. 25(a) and (c) demonstrate that in EP and KH cases, high stress locates at the inhomogeneous zone of the sample where failure occurs. Meanwhile, Fig. 25(b) and (d) exhibit most inhomogeneous zone in the sample reached limit state, i.e. the ratio between Mises stress and local material strength σe/σs(σu) is approximated as 0.9 ~ 1. According to the mechanism of plastic collapse, the plastic hinge produced by these high stress zone leads to the plastic failure of

5.3. Inhomogeneous kinematic hardening behavior Impact of kinematic hardening behavior on the failure of structure is complicated due to Bauschinger effect. According to the discussion in Section 5.1, in a degenerated case (homogeneous case), the proposed method predicts the kinematic hardening effect in shakedown analysis precisely. So kinematic hardening behavior in inhomogeneous cases can be further discussed. In example I, the shakedown domain of inhomogeneous pipe shown in Fig. 13 exhibits a similar profile compared with the homogeneous case shown in Fig. 9, that is, kinematic hardening enhances the ratchet resistance but have neglectable influence on the alternating plasticity failure boundary. This character has been intensively studied and reported by Refs. [10–12,20,23] and has been reproduced by BSIM. Furthermore, based on the results of validation points shown in Fig. 14, the rationality of shakedown domain shown in Fig. 13 can be further 9

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

In example II, Fig. 17(a) has proven that kinematic hardening enlarges the limit load domain. In inhomogeneous case, according to the results in Table.4, the enhancement of limit load has also been reproduced by the proposed method. Besides, from Figs. 20 and 21 one can found the proposed method and FVM suggests identical failure zone and similar stress contours in all cases, whereas the stress contour in EP case and KH case exhibit differences, implying a further stress redistribution resulting from kinematic hardening behavior. Comparing Fig. 21(a) & (c) with Fig. 21(b) & (d), one can conclude that in strain hardening (SH) case, the proposed method predicts same failure zone and stress distribution compared with results from FVM, implying the present method predicts the stress field and failure point reasonably. The kinematic hardening induced stress redistribution further changed the limit load of the plate. As is shown in Table 4, in SH case, limit load under P1 and P2 are identical. The reason is that the pre-defined ultimate strength shown in Fig. 18(b) exhibits the same strength in zone A and zone B. Additionally, limit load in KH case suffer non-uniform enhancement compared with that in EP case, which means that the effect of kinematic hardening depends on local kinematic hardening ratio σuGu(xi)/σsGs(xi). In example III, the effect of kinematic hardening is revealed by the increasing of smooth sample’s limit load in KH case. Based on the limit state distribution shown in Fig. 25(c), the proposed method successfully predict the limit stress field, where most of the inhomogeneous zone reached limit state, see Fig. 25(d). Comparison between Fig. 25(c) and (d) demonstrate the stress distribution pattern changing from EP case to KH case, which leads to load-bearing capacity change. As a result, it’s concluded that the proposed method predicts the characters of kinematic hardening behavior in inhomogeneous limit and shakedown analysis reasonably.

Fig. 15. Dimensions and loads of the plate with central hole.

5.4. Closure Limit and shakedown load of three typical examples were calculated and discussed intensively. The characters of the examples are listed in Table.6. It’s clear that the examples with various complexity can cover most engineering cases in modern industries. Based on aforementioned discussions, for all three examples, the proposed method BSIM has predicted the limit and shakedown loads, the stress field, the characters of kinematic hardening behavior and the failure process of inhomogeneous structure precisely. As a result, it’s concluded that the proposed method is of validity and can predict the limit and shakedown load of kinematic hardening structure made of arbitrary inhomogeneous material reasonably. Fig. 16. Finite element model of the plate with central hole.

6. Conclusions verified. It’s demonstrated in Fig. 14(a) that equivalent plastic strain of the pipe under validation point A & E increase firstly and then kept constant, exhibiting an elastic shakedown response. Meanwhile, load combinations B & F result in unconfined incremental plastic deformation, implying a ratchet failure. Similar conclusion can be draw for KH ratchet boundary from Fig. 14(b), indicating the proposed method predicts EP & KH ratchet reasonably. Fig. 14(c) & (d) demonstrate that load combinations M & I & O & K remain elastic shakedown response while load combinations N & P exhibit typical alternating plasticity response. Therefore, the realistic alternating boundary locates between the validation points M & I & O & K and N & P. As a result, proposed method can be considered as reasonable in this case. What’s more, for load combinations J & L, structure’s response exhibits both the character of ratchet and alternating plasticity, implying these load combinations locate near the branch point of the shakedown boundary. According to Fig. 13, this branch point locates near load combinations I(K) & J(L), which is another evidence for the validity of the proposed method.

In the present paper, a numerical technique named Back Stress Iteration Method (BSIM) is proposed for limit and shakedown analysis with arbitrary inhomogeneity. The proposed method is developed based on basis reduction method and a novel back stress construction algorithm. Three typical examples were carried out to validate BSIM. Main remarks are concluded as following. (1) The examples proved that BSIM can predict limit and shakedown load of kinematic hardening structure with arbitrary inhomogeneity reasonably. (2) BSIM extends the advantage of basis reduction method to inhomogeneous kinematic hardening shakedown problem, with which the computational burdens of shakedown analysis can be significantly reduced. (3) BSIM is compatible with arbitrary material properties distribution pattern, therefore its potential in composite structure engineering is huge. 10

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 17. Limit and shakedown load of the homogenous plate.

Fig. 18. Stress field of homogeneous plate under different load cases. 11

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 19. Material strength distribution in the plate.

Fig. 20. Limit stress contour of inhomogeneous EP plate. 12

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 21. Limit stress contour of inhomogeneous SH plate.

Table 4 Limit load of the inhomogeneous plate. P1: P2

Present method EP

FVM EP

Present method KH

FVM IH

1:0 1:1 0:1

189.409 191.367 173.622

189.443 189.467 173.125

250.797 280.480 250.010

258.309 286.605 258.309

Acknowledgement This work was supported by the National Science Foundation of China (Grand No. 51905173). The first author would thank his PhD Tutor Prof. Zhiping Chen for everything he did for him during last five years.

Fig. 23. Finite element model of smooth tensile sample.

Fig. 22. Dimensions and loads of the smooth tensile sample. 13

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Fig. 24. Strength distribution in the tensile sample. Table 5 Limit load of the smooth tensile sample.

Table 6 Complexity of examples.

Material state

Limit load (MPa)

No. of Example

Load

Stress state

Inhomogeneity

Homogeneous case EP Homogeneous case KH theoretical solution EP theoretical solution KH Inhomogeneous EP Inhomogeneous KH Homogenized EP Homogenized KH

49.9815 74.8654 50 75 43.6153 51.1325 41.6598 54.0555

I

Mechanical & Thermal Mechanical Mechanical

3-dimentional

1-dimentional

2-dimentional 1-dimentional

2-dimentional 3-dimentional random

II III

Declaration of Competing Interest The authors declare that they have no conflict of interest. Data availability. The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Fig. 25. limit Stress contour in the inhomogeneous tensile sample. 14

Composite Structures xxx (xxxx) xxxx

S. Huang, et al.

Appendix A. Supplementary data

[10] Stein E, Zhang G, König J. Shakedown with nonlinear strain-hardening including structural computation using finite element method. Int J Plast 1992;8:1–31. [11] Heitzer M, Pop G, Staat M. Basis reduction for the shakedown problem for bounded kinematic hardening material. J Global Optim 2000;17:185–200. [12] Huang S, Wan F, Jiao P, Chen Z. A modified basis reduction method for limited kinematic hardening shakedown analysis under complex loads. Mech Based Des Struct Mach 2018;46:85–100. [13] Simon J, Kreimeier M, Weichert D. A selective strategy for shakedown analysis of engineering structures. Int J Numer Meth Eng 2013;94:985–1014. [14] Magoariec H, Bourgeois S, Débordes O. Elastic plastic shakedown of 3D periodic heterogeneous media: a direct numerical approach. Int J Plast 2004;20:1655–75. [15] Zhang H, Liu Y, Xu B. Plastic limit analysis of ductile composite structures from micro- to macro-mechanical analysis. Acta Mech Solida Sin 2009;22:73–84. [16] You JH, Kim BY, Miskiewicz M. Shakedown analysis of fibre-reinforced copper matrix composites by direct and incremental approaches. Mech Mater 2009;41:857–67. [17] Chen G, Bezold A, Broeckmann C, Weichert D. On the statistical determination of strength of random heterogeneous materials. Compos Struct 2016;149:220–30. [18] Tornabene F, Fantuzzi N, Bacciocchi M. Linear static behavior of damaged laminated composite plates and shells. Materials 2017;10:811. [19] Wang B, Mai Y, Noda N. Fracture mechanics analysis model for functionally graded materials with arbitrarily distributed properties. Int J Fract 2002;116:161–77. [20] Simon J. Direct evaluation of the limit states of engineering structures exhibiting limited, nonlinear kinematical hardening. Int J Plast 2013;42:141–67. [21] Kanno Y. A note on formulations of static shakedown analysis with bounded kinematic hardening. Mech Res Commun 2016;74:57–9. [22] Liu Y. Lower bound shakedown analysis by the symmetric Galerkin boundary element method. Int J Plast 2005;21:21–42. [23] Chen S, Liu Y, Li J, Cen Z. Performance of the MLPG method for static shakedown analysis for bounded kinematic hardening structures. Eur J Mech A Solids 2011;30:183–94. [24] François A, Abdelkader H, Said M, et al. Application of lower bound direct method to engineering structures. J Global Optim 2007;37:609–30.

Supplementary data to this article can be found online at https:// doi.org/10.1016/j.compstruct.2019.111641. References [1] Huang S, Chen Z, Li Y, Zhang D. Impact of hot bending on the high-temperature performance and hydrogen damage of 2.25Cr-1Mo-0.25V Steel. J Mater Eng Perform 2019;28:567–77. [2] Garbatov Y, Tekgoz M, Soares CG. Experimental and numerical strength assessment of stiffened plates subjected to severe non-uniform corrosion degradation and compressive load. Ships Offshore Struct 2017;12:461–73. [3] Neugebauer R, Schieck F, Rautenstrauch A, Bach M. Hot sheet metal forming: the formulation of graded component characteristics based on strategic temperature management for tool-based and incremental forming operations. CIRP J Manuf Sci Technol 2011;4:180–8. [4] Reza M, Mohsen D, Ali N, Hamid D. Thick-walled functionally graded material cylinder under thermo-mechanical loading. 2018 9th International Conference on Mechanical and Aerospace Engineering (ICMAE). 2018. p. 505–10. [5] Friend CM. The effect of matrix properties on reinforcement in short alumina fibrealuminium metal matrix composites. J Mater Sci 1987;22:3005–10. [6] Pierard O, LLorca J, Segurado J, Doghri I. Micromechanics of particle-reinforced elasto-viscoplastic composites: finite element simulations versus affine homogenization. Int J Plast 2007;23:1041–60. [7] Melan E. Zur Plastizität des räumlichen Kontinuums. Ingenieur-Archiv. 1938;9:116–26. [8] Koiter W. A new general theorem on shake-down of elastic-plastic structures. Proceeding of eighth international congress of applied mechanics. 1956. p. 220–30. [9] Weichert D, Grossweege J. The numerical assessment of elastic plastic sheets under variable mechanical and thermal loads using a simplified 2-surface yield condition. Int J Mech Sci 1988;30:757–67.

15