The material non linear analysis of 2D strutures using a radial point interpolation method

FACULDADE DE ENGENHARIA DA UNIVERSIDADE DO PORTO The Material Non Linear Analysis of 2D Structures Using a Radial Point Interpolation Method Henrique Manuel Sousa Duarte Mestrado Integrado em Engenharia Mecânica Supervisor: Jorge Belinha (Post-PhD Researcher and Invited Auxiliar Professor) Co-Supervisor: Renato M. Natal Jorge (Associated Professor) Co-Supervisor: Lúcia M. J. S. Dinis (Associated Professor) July 2014 The Material Non Linear Analysis of 2D Structures Using a Radial Point Interpolation Method Henrique Manuel Sousa Duarte Mestrado Integrado em Engenharia Mecânica July 2014 Agradecimentos Em primeiro lugar, gostaria de agradecer ao Professor Jorge Belinha, orientador. Obrigado por ter despertado o meu interesse pelo estudo dos métodos sem malha. Obrigado pela disponibilidade que demonstrou ao longo da realização desta tese, assim como ao longo dos trabalhos e apresentações que a precederam. Ao Professor Renato Natal Jorge e Professora Lúcia Dinis, agradeço a oportunidade que me foi dada para entrar no mundo da investigação. O convite pessoal e a visita guiada aos diferentes projetos em desenvolvimento ficarão para sempre gravados na minha memória como um dos momentos altos do meu percurso académico. Aos colegas de curso e amigos José Andrade, Rafael Tavares, João Ferreira, Francisco Rua, Tiago Ramos, entre muitos outros, agradeço todo o tempo passado em grupo, quer em ambiente sério e de trabalho, quer nos momentos de descontração. À minha família, em particular aos meus pais e padrinhos, o meu profundo agradecimento pelo apoio incondicional ao longo deste trabalho e de todo o percurso até ele. Obrigado por criarem todas as condições necessárias para que fosse bem sucedido. Obrigado pelo incentivo e ânimo constantes. Finalmente, gostaria de agradecer a todos os que, de uma forma ou de outra, me ajudaram na realização deste trabalho. i ii Abstract In this work, the Radial Point Interpolation Method (RPIM) is used to study several solid mechanics problems, starting with elastic benchmark examples, then elastoplastic problems, with monotonic and non monotonic loading conditions. Additionally, examples of trabecular and cortical bone are analysed. Meshless methods have been the focus of interest in the past few years, mainly because of the advantages they have when compared to the well known and established Finite Element Method. Some of these advantages are the production of smoother stress fields and the more flexible discretization technique. In the RPIM, the nodal connectivity is imposed using the concept of neighbour nodes, which are radially searched around the interest node. This creates variable sized influence domains that contain the neighbour nodes of the interest node, and that are used to create the interpolation functions. The RPIM is not a trully meshless method, since a background mesh is required to perform the numerical integration of the interpolation functions. This mesh can be fitted to the domain of the problem or blindly regular, with the subsequent exclusion of integration points outside the domain. The interpolation functions used in the RPIM contain a polynomial basis and a radial basis. The Radial Basis Function (RBF) used in this work is the Multiquadric RBF. Since the interpolation functions possess the delta Kronecker property, the imposition of natural and essencial boundary conditions is simple and direct. Regarding the non-linear analysis, small deformations and elastoplastic material behaviour are considered. In order to solve the non-linear equation system, an incremental and iterative method is used, a modified Newton-Raphson method. In the initial stiffnes method (KT0), the stiffness matrix is calculated only once at the start of the analysis, which reduced the computational cost, but increases the number of iterations. The yield criterion used for the non-linear analysis is the Von Mises criterion, and the forward-Euler procedure is used to return the stresses to the yield surface. Several linear elastic benchmark examples are solved, in order to assess the quality of the meshless method. After that, non-linear examples are analysed, first with monotonic loading and then with non monotonic loading. The final application is the study of the trabecular and cortical bone using elastoplastic material behaviour. The obtained results suggest that the RPIM is an accurate and reliable meshless method. iii iv Resumo Neste trabalho, o "Radial Point Interpolation Method" (RPIM) é usado no estudo de diversos problemas de mecânica dos sólidos, começando com exemplos lineares de referência, passando depois para problemas não lineares, com cargas monotónicas e não monotónicas. Álem disso, exemplos de osso trabecular e cortical são analisados. Os métodos sem malha têm sido o foco da atenção nos últimos anos, principalmente devido às vantagens que apresenta quando comparado com o conhecido e estabelecido Método dos Elementos Finitos. Algumas destas vantangens são a produção de campos de tensões suaves e maior flexibilidade de discretização. No RPIM, a conetividade nodal é imposta usando o conceito de nós vizinhos, que são radialmente procurados à volta do nó de interesse. Assim são criados domínios de influência de tamanho variável que contêm os nós vizinhos do nó de interesse, usados para a criação das funções de interpolação. O RPIM não é um método sem malha puro, uma vez que é necessária a criação de uma malha de fundo para a integração numérica das funções interpoladoras. A malha pode ser ajustada ao modelo ou regular, sendo necessário neste caso a remoção posterior dos pontos de integração fora do domínio. As funções de interpolação usadas no RPIM são formadas por uma base polinomial e por uma base radial. A base radial usada neste trabalho é a multiquadrica. Uma vez que as funções de forma possuem a propriedade de delta Kronecker, a imposição das condições de fronteira naturais e essenciais é simples e direta. No que diz respeito à análise não linear, foram consideradas pequenas deformações e comportamento material elastoplástico. De forma a resolver o sistema de equações não linear, é usado um método incremental e iterativo, o método de Newton-Raphson modificado. No método de rigidez inicial (KT0), a matriz de rigidez é calculada apenas uma vez no início da análise, o que reduz o custo computacional mas aumenta o número de iterações. O critério de cedência usado para a análise não linear é o critério de Von Mises, e o processo "forward-Euler" é utilizado para retornar as tensões à superfície de cedência. Vários exemplos lineares de referência são resolvidos, de forma a avaliar a qualidade do método sem malha. Depois disso, são analisados exemplos não lineares, primeiro com cargas monotónicas e depois com cargas não monotónicas. A última aplicação é o estudo do osso trabecular e cortical utilizando comportamento material elastoplástico. Os resultados obtidos sugerem que o RPIM é um método sem malha preciso e fiável. v Introduction In the last few years, many meshless methods using interpolation functions have been developed, with the main goal of making the imposition of the natural and essential boundary conditions easier. The main ones are the Point Interpolation Method (PIM) [13,14], the Radial Point Interpolation Method (RPIM) [15,16], the Natural Neighbour Finite Element Method (NNFEM) [17,18], the Meshless Finite Element Method (MFEM) [19] and more recently the Radial Natural Element Method (NREM) [20,21]. 1.2 Objective The main objective of this work is the analysis of elastoplastic materials. For this, an existing MATLAB code [22] was studied, and additional code was created in order to perform the non linear analysis. The code includes the KT0 algorithm, the stress return algorithm and the analysis of non monotonic loads. 1.3 Thesis Structure The structure of this work is described as follows. In Chapter 2, the meshless method used in this work, the RPIM, is presented. Chapter 3presents the elastic and elastoplastic mechanical fundamentals. The nonlinear solution algorithm and the stress return algorithm used to solve elastoplastic problems are described in Chapter 4. Chapter 5presents the analysis of linear elastic problems, in order to evaluate the performance of the meshless method used. In Chapter 6, the analysis of some benchmark nonlinear problems is presented, along with the study of some non monotonic loading cases. In addition, the example of the trabecular bone and polymer analogue material is analysed, as well as the cortical bone case. Finally, in Chapter 7the main conclusions of this work are presented, along with perspectives for future works. 2 Chapter 2 Meshless Method In the past few years, meshless methods have come into focus of interest, specially in the engineering community. Their development is related to the need to overcome some limitations found in the finite element methods (FEM). In the meshless methods, the nodes can be arbitrarily distributed inside the domain, since the field functions (such as displacement and stresses) are approximated within an influence domain, rather than an element [2,4,23]. These influence domains may and must overlap, as opposed to FEM, in which there is a no-overlap rule between elements [24,25]. In this chapter, the Radial Point Interpolation Method (RPIM) is presented. 2.1 Radial Point Interpolators The Radial Point Interpolators (RPI) have their origin in the Point Interpolation Method (PIM) [13]. In this method, polynomial interpolants possessing the delta Kronecker property are constructed, based only on a group of arbitrarily distributed points. This technique, however, has several numerical problems, for example the perfect alignment of the nodes produces a singular solution in the interpolation functions. To overcome these problems, a radial basis function was introduced in the construction process of the interpolation functions to stabilize the method, and the PIM evolved to the Radial Point Interpolation Method (RPIM) [15]. The RPIM was initially developed to perform data surface fitting. Later, Kansa’s developments [26] to the RBF allowed it to be used for diferential equation solving. Unlike Kansa’s algorithm, which uses the concept "global domain", the RPIM uses the "influence domain". This generates sparse and banded stiffness matrices, more adequate to complex geometry problems. 2.2 Influence Domain and Nodal Connectivity Consider a problem domain Ω, bounded by Γ, discretized in a set of randomly placed nodes N={n0,n1,...,nN}∈R2, as shown in figure 2.1. 3 Meshless Method Γ Ω Figure 2.1: Discretization of the Problem Domain in Several Randomly Distributed Nodes In early works, the nodal connectivity in the RPIM is obtained by overlaping the influence domains of the nodes [13,14]. These domains are constructed by searching nodes inside a fixed area (or a fixed volume). This concept is very simple, however the irregular boundaries or different node densities inside the model can lead to unbalanced influence domains. It would be ideal if the influence domains in the problem contained the same number of nodes [22]. The solution to this problem is the use of a variable influence domain. Instead of a fixed size, this influence domain has a fixed number of nodes n, while the size may vary, as shown in figure 2.2. For each interest point, a radial search around the point itself is performed, and the nclosest nodes are defined. xi xj ni=5 nj=5 di6=dj di dj Figure 2.2: Variable Influence Domain Previous work [27] suggests that the number of nodes inside each influence domain should vary between n∈[16,25]. 2.3 Numerical Integration In the present work, the Gauss-Legendre integration scheme is considered. The solid domain is divided in a regular grid, and each grid cell is filled with integration points, respecting the Gauss-Legendre quadrature rule. Figure 2.3 shows an example for a single cell. 4 Meshless Method The initial quadrilateral cell (figure 2.3(a)) is transformed in an isoparametric square (figure 2.3(b)). The Gauss-Legendre quadrature points are placed inside the isoparametric square. In this case, as in the whole work, a 2×2 quadrature is used. The cartesian coordinates of the integration points are obtained by using the isoparametric interpolation functions (figure 2.3(c)). The weight of each integration point is the product of the isoparametric weight and the determinant of the Jacobian matrix’s inverse for the respective cell [22]. y x n1 n2 n3 n4 (a) η ξ n1 n2n3 n4 (b) y x n1 n2 n3 n4 (c) Figure 2.3: Gauss-Legendre Integration [22]: (a) Initial Cell; (b) Transformation of the Initial Cell in an Isoparametric Square and Application of 2×2 Quadrature Rule; (c) Return to the Initial Cell. If the solid domain is fairly regular, a regular integration mesh that fits the domain can be constructed, and no additional post-treatment is required (figure 2.4(a)). Although, if the domain is irregular, a regular mesh can be created, as figure 2.4(b) shows, and after the integration points are placed, the ones that are outside the domain are removed. (a) (b) Figure 2.4: Background Integration Mesh [22]: (a) Fitted Mesh; (b) Regular Blind Mesh. Consider the function F F F(x)defined in the domain Ω. The numerical integration can be expressed by, Z Ω F F F(x)dΩ= ng ∑ i=1 _ wiF F F(xi)(2.1) where _ wiis the weight of the integration point xi. 5 Meshless Method 2.4 Interpolation Functions In this work, the punctual radial type functions - RPI (radial point interpolators) are used [11,15]. Consider a function u u u(x x x), defined in the domain Ω, discretized by a set of Nnodes. It is assumed that only the nodes inside the influence domain of the interest point x x xIhave effect on the function u u u(x x x). Using a radial basis function, the function u u u(x x x)passes through all the nodes of the influence domain. The value of the function u u u(x x x)at the interest point x x xIis obtained by, u u u(x x xI) = n ∑ i=1 Ri(x x xI)ai(x x xI)+ m ∑ j=1 pj(x x xI)bj(x x xI) = R R RT(x x xI),p p pT(x x xI)(a a a b b b)(2.2) where Ri(x x xI)is the RBF and nis the number of nodes inside the influence domain of the interest point x x xI. The coefficients ai(x x xI)and bi(x x xI)are, respectively, non constant coefficients of Ri(x x xI)and pj(x x xI). The monomials of the polynomial basis are defined by pj(x x xI)and mis the basis monomial number. The vectors of equation 2.2 are defined as, R R RT(x x xI) = {R1(x x xI),R2(x x xI),...,Rn(x x xI)}(2.3) p p pT(x x xI) = {p1(x x xI),p2(x x xI),...,pm(x x xI)}(2.4) a a aT(x x xI) = {a1(x x xI),a2(x x xI),...,an(x x xI)}(2.5) b b bT(x x xI) = {b1(x x xI),b2(x x xI),...,bm(x x xI)}(2.6) Several known RBF’s are well studied and developed [15,28]. In this work, the Multiquadric (MQ) function is used, initially proposed by Hardy [29]. The MQ-RBF is defined as, R(rIi) = r2 Ii +c2p(2.7) where cand pare two shape parameters that require an optimization study, since the variation of these parameters greatly affects the performance of the RBF. Previous work suggests that c∼ =0 and p∼ =1 [11,30]. The variable of the RBF is the Euclidian norm rIi, which defines the distance between the interest point x x xIand the neighbour node x x xi, ri j =q(xI−xi)2+(yI−yi)2(2.8) 6 Meshless Method The polinomial basis added must be complete to assure that the interpolation of the RBF matrix is invertible [11]. The polinomial basis that can be used for the two-dimensional space are, Null Basis - x x xT={x,y};p p pT(x x x) = {0};m=0,(2.9) Constant Basis - x x xT={x,y};p p pT(x x x) = {1};m=1,(2.10) Linear Basis - x x xT={x,y};p p pT(x x x) = {1,x,y};m=3,(2.11) Quadratic Basis - x x xT={x,y};p p pT(x x x) = 1,x,y,x2,xy,y2;m=6.(2.12) In this work, the constant polinomial basis is used. There is an additional requirement that the polinomial basis must satisfy in order to obtain a unique solution [11], n ∑ i=1 pj(xi)ai(xi) = 0,j=1,2,...,m(2.13) Therefore, a new equation system can be defined, (u u us 0 0 0)=G G G(a a a b b b)(2.14) where G G Gis a matrix defined by, G G G="R R RQP P Pm P P PT m0 0 0#(2.15) being R R RQthe moment matrix of the RBF, R R RQ=      R(r11)R(r12)··· R(r1n) R(r21)R(r22)··· R(r2n) . . .. . ..... . . R(rn1)R(rn2)··· R(rnn)       (2.16) and P P Pmthe moment matrix of the polynomial basis, P P Pm=      P1(x1)P2(x1)··· Pm(x1) P1(x2)P2(x2)··· Pm(x2) . . .. . ..... . . P1(xn)P2(xn)··· Pm(xn)       (2.17) 7 Meshless Method Solving equation 2.14, (a a a b b b)=G G G−1(u u us 0 0 0)(2.18) Substituting the previous equation in the interpolation function (equation 2.2), u u u(x x xI) = R R RT(x x xI),p p pT(x x xI)G G G−1(u u us 0 0 0)=ϕ(x x xI)u u us(2.19) where ϕ(x x x)is the interpolation function defined as, ϕ(x x xI) = R R RT(x x xI),p p pT(x x xI)G G G−1={ϕ1(x x xI),ϕ2(x x xI),...,ϕn(x x xI)}(2.20) Figure 2.5 shows a schematic representation of the interpolation functions for a one-dimensional space [31]. It is visible that the functions have a null value on all nodes of the influence domain except the interest point x x xI. The interpolation function has a unit value for this point. Therefore, the RPIM function possesses the delta Kronecker property. Its derivatives are easily obtainable [11,30,32]. (a) (b) (c) xI xI xI xI ϕ(xI) ϕ(xI) ϕ(xI) Figure 2.5: Interpolation Functions for a One-Dimensional Space: (a) 2 Node Influence Domain; (b) 4 Node Influence Domain; (c) 6 Node Influence Domain [31] 8 Chapter 3 Mechanical Fundamentals The science of solid mechanics is the foundation for the modelling and designing of structural systems. It defines the relations between stresses (which are used to verify how far the structure is from yielding), and strains (caused by the loading conditions of the structure) [33]. In relation to mechanical properties, materials can be divided in two groups. Isotropic materials have the same mechanical properties in all directions, as oposed to anisotropic materials, in which the mechanical properties vary depending on the direction. The first group can be defined by two independent material constants (usually the Young’s modulus and the Poisson’s ratio). For the second group, multiple constants may be required, depending on the level of anisotropy (transversely isotropic, orthotropic) [34]. Regarding the relationship between deformation and force aplied to a solid, materials can be defined as elastic (the deformation disappears when the solid is unloaded) or plastic (the deformation in the solid can’t be fully recovered after unloading). In this work, both cases are considered. The boundary conditions are an important consideration in solid mechanics [33,34]. These can be applied through forces or imposed displacements. In this work, only static (or quasistatic) forces are considered, which means the stress, strain and displacement will not be a function of time. 3.1 Problem Formulation 3.1.1 Displacement Field Consider a two dimenionsal problem, where the domain Ωis bounded by Γ. In case of the application of an external force, the points in domain Ωwill change to the domain Ω0. Figure 3.1 shows the initial particle Pmoving to position P0[35]. 9 Mechanical Fundamentals The vector u u u(u,v)represented in figure 3.1 is a function of the coordinates (x,y), as shown in equation 3.1. It represents the displacement field for the two dimensional body. u u u(u,v) = (u(x,y) v(x,y))(3.1) x y Γ Ω Γ0 Ω0 P P0 u u u= (u,v) Figure 3.1: Displacement Field 3.1.2 Strain Field Strain is a measure of deformation, which can be caused by external loads, body forces or internal forces, such as the ones caused by temperature gradients. It can be defined as the change of displacement per unit length. x y A0 D0 C0 B0 A DC B u v ∆x ∆y u+∂u ∂y∆y v+∂v ∂y∆y u+∂u ∂x∆x v+∂v ∂x∆x Figure 3.2: Cartesian Components of Strain for an Infinitesimal Material Element 10 Mechanical Fundamentals Strains can be divided into normal strains (related to dilatations or contractions) or shear strains (related to distortions). Figure 3.2 represents the strain field components (normal and shear strain) for the two dimensional material element. The strain matrix can be defined by equation 3.2, ε ε ε="εxx εxy εyx εyy #(3.2) where the normal strains are given by, εxx =∂u ∂x;εyy =∂v ∂y(3.3) and the shear strain by (note that γi j =2×εi j), γxy =γyx =∂u ∂y+∂v ∂x(3.4) Since the strain matrix is symmetric, the Voigt notation can be implemented, turning equation 3.2 into, ε ε ε=     εxx εyy εxy      (3.5) Replacing equations 3.3 and 3.4 in the previous one gives, ε ε ε=     εxx εyy εxy      =     ∂u ∂x ∂v ∂y ∂u ∂y+∂v ∂x      (3.6) Equation 3.6 can be represented as a product of a partial differential operator L L Land the displacement vector u u u, ε ε ε=L L L·u u u=   ∂ ∂x0 0∂ ∂y ∂ ∂y ∂ ∂x   ·(u v)(3.7) 3.1.3 Constitutive Equations When dealing with elastic materials, the constitutive equation that expresses the relationship between stresses and strains is the Hooke’s Law, given by, σ σ σ=c c c·ε ε ε(3.8) 11 Mechanical Fundamentals 3.5 Yield Criterion The yield criterion mathematically defines the stress values that match the elastic limit behaviour of the material [40]. It is a function of the stress state and the loading history, and can be generally expressed by, F(σ σ σ,ε ε εp,κ) = f(σ σ σ,ε ε εp,κ)−σY(κ) = 0 (3.39) where f(σ σ σ,ε ε εp,κ)is the yield function, dependent on the stress state ε ε εp, the plastic strain ε ε εp and a hardening parameter κ. The stress σY(κ)represents the yield stress of the material. Equation 3.39 can be simplified in the case of an isotropic material with an yield stress independent of κ[41], F(σ σ σ) = f(σ σ σ)−σY=0 (3.40) In this work, Von Mises yield criterion is used, which is commonly applied to metals. Figure 3.4 shows a representation of the Von Mises and Tresca yield criteria in the Westergaard space. It is visible that the Tresca yield criterion presents edges in the surface, which are hard to deal with in terms of computation (due to discontinuities and singularities). The Von Mises criterion has a smooth surface, therefore the computation is easier when compared with the Tresca criterion. (a) σ3 σ1σ2 Tresca Yield Surface Von Mises Yield Surface q2 3σY (b) Figure 3.4: Von Mises and Tresca Yield Criterion: (a) General View [31]; (b) Deviatoric Plane. Von Mises suggested that the yielding occured when the second deviatoric stress invariant (I2) reached a critical value [42], defined as, I2=σY √3(3.41) The previous equation defines the surface of a cilinder circumscribed in the Tresca hexagon, and whose interception with the πplane is the circle represented on figure 3.4(b). 18 Mechanical Fundamentals The Von Mises criterion can also be expressed as a function of the effective stress ¯ σ, ¯ σ2=σ2 Y=3I2(3.42) The second deviatoric stress invariant I2can be defined in terms of the cartesian components of stress, ¯ σ=r1 2h(σxx −σyy)2+(σyy −σzz)2+(σzz −σxx)2+6σ2 yz +σ2 zx +σ2 xyi(3.43) Substituing the previous equation in equation 3.40, the following equation is obtained, F(σ) = r1 2h(σxx −σyy)2+(σyy −σzz)2+(σzz −σxx)2+6σ2 yz +σ2 zx +σ2 xyi−σY=0 (3.44) The yield function can be defined by many different expressions, each one having a distinct geometrical representation in the Westergaard space. The material to be simulated has to be taken into account when deciding the yield function to use [43]. 3.6 Plastic Flow By observing equation 3.40, it is possible to conclude that if f<σY, the material is in the elastic domain, and f=σYrepresents the limit of the elastic domain and start of the plastic behaviour. Once this point is reached, the subsequent behaviour of the material is conditioned by the value of the variation of the yield function fwith the stress σ σ σ. This variation can be defined as, d f =∂f ∂σ σ σT dσ σ σ(3.45) where ∂f/∂σ σ σis the gradient of f, and consequently an orthogonal vector to the yield surface for a given stress σ σ σ, as shown in figure 3.5. The outcome of equation 3.45 can be divided into three situations. If d f <0, the stress point is inside the yield surface, and elastic unloading has occured. The relation between stress and strain is linear. If d f =0, the stress point is on the yield surface. This can represent a perfectly plastic state, in case the material has no hardening parameter κ, or otherwise the start of the plastification state. Finally, d f >0 corresponds to the plastic loading state, with the stress remaining in an expanding yield surface. At this point, the material is in a plastic flow state [41]. Introducing the concept of the plastic potential function, g(σ σ σ,ε ε εp,κ), the flow rule can be defined as, dε ε εp=dλ∂g ∂σ σ σ(3.46) 19 Mechanical Fundamentals σ2 σ1 Elastic Behaviour Plastic Behaviour Yield Surface d f dσ d f dσ1 d f dσ2 Figure 3.5: Orthogonality Condition in a Two-Dimensional Stress Space where dλis the plastic strain-rate multiplier [43]. The plastic multiplier defines the magnitude of the plastic strain increment vector, while the gradient of the plastic potential defines the direction. In this work, an associative flow rule is used. This means that the plastic potential function is actually the yield function f(σ σ σ,ε ε εp,κ), and so equation 3.46 can be rewritten as, dε ε εp=dλ∂f ∂σ σ σ(3.47) As mentioned earlier, the vector ∂f ∂σ σ σis normal to the yield surface. Therefore the plastic strain increment vector is also normal to the yield surface. This orthogonality ensures a unique solution [39]. 3.7 Hardening Rule The hardening rule has the objective of describing the conditions for a new plastic flow, given that the plastic phase of the material has been achieved. These conditions are essential, since the yield surface may change its size and shape with the increasing plastic deformation. Two types of hardening rules may be defined, depending on how the hardening parameter κis calculated. These are strain hardening and work hardening. Figure 3.6 presents some strain hardening models. A perfect plastic material has no hardening parameter κand its yield surface never changes size or shape. 20 Mechanical Fundamentals σ1 σ2 1 2 (a) αi j σ1 σ2 1 2 (b) αi j σ1 σ2 1 2 (c) Figure 3.6: Hardening Rules: (a) Isotropic Hardening; (b) Kinematic Hardening; (c) Mixed Hardening. In the isotropic hardening rule, figure3.6(a), it is assumed that the yield surface expands uniformly without distortion or movement. Figure 3.6(b) represents the Kinematic Hardening Rule. In this case, the yield surface evolves by moving in relation to its initial origin, without changing its size. With the different types of movements (expansion, translation and even rotation), many hardening models can be defined [39], as is the example of the mixed hardening rule, represented in figure 3.6(c). Generically, the yield criterion can be written as, F(σ σ σ,ε ε εp,α α α,κ) = f(σ σ σ,ε ε εp,α α α,κ)−k2(κ) = 0 (3.48) being k2the yield surface size, f(σ σ σ,ε ε εp,α α α,κ)the yield surface shape and α α αthe centre of the yield surface. In this work, the isotropic hardening rule is used. Therefore, equation 3.48 can be simplified, f(σ σ σ) = k2(κ)(3.49) By applying the Von Mises criterion, the previous equation becomes, F(σ σ σ,k) = 3I2−σ2 Y(κ)(3.50) The hardening parameter κcan be obtained by the effective plastic strain ¯ εp, ¯ εp=Zd¯ εp=Zr2 3dε ε εpdε ε εp(3.51) Therefore, the κparameter is simply defined as, κ=¯ εp(3.52) 21 Mechanical Fundamentals 22 Chapter 4 Nonlinear Solution Algorithms The elastoplastic analysis of structures requires the solution of a set of non linear equations. These equations can’t be solved by direct methods. Therefore, two types of numerical algorithms can be employed to obtain a solution, the incremental method and the incremental and iterative method, also know as the Newton-Rapson method. The incremental algorithm is very simple and intuitive. The non linear problem is simplified to a series of linear solutions, where the loads are applied incrementally until the desire value is reached. However, the method usually has a bad performance for considerable non linearities, leading to great errors [42]. The Newton-Rapson method overcomes the problems of the previous method, since in every load increment, there is an iteration process that reduces the transition error between load increments to an insignificant value. Figure 4.1(a) represents the incremental and iterative Newton-Rapson method. Even though this method generates very good results in relation to the convergence to the final solution, this method is a very slow numerical procedure, since the stiffness matrix and its inverse must be recalculated in every iteration process. f u ∆f ∆f Prediction Iterations Predictor Iterations u0u1u0 2u1 2u2 2u2 (a) f u ∆f ∆f Prediction Iterations Predictor Iterations (b) f u ∆f ∆f (c) Figure 4.1: Newton-Rapson Methods: (a) Classic Newton-Rapson; (b) Initial Increment Stiffness Variation (KT1); (c) Initial Stiffness Variation (KT0). 23 Nonlinear Solution Algorithms In order to reduce the computational cost, several modified versions of the Newton-Rapson were developed. The first variant is the initial increment stiffness method (KT1), represented in figure 4.1(b). In this variation, the stiffness matrix is calculated at the start of each load increment, and remains constant throughout the entire iteration process. This results in a lower computational cost, however it has slower convergence. Further simplifications can be made to the KT1 method. Figure 4.1(c) represents the initial stiffness method, in which the stiffness matrix is only calculated in the first increment and used throughout the whole process. When compared with the KT1 method, the KT0 has lower computational cost, since there is only one calculation and inversion of the stiffness matrix. However, more iterations are required to achieve the value of the incremental load. 4.1 KT0 Algorithm As mentioned previously, the KT0 algorithm is used in this work, which is described in this section. The main characteristic of this method is the single calculation of the stiffness matrix and its inverse at the start of the first load increment. Before the KT0 is ready to initiate, initial data must be introduced, such as: problem dimensions and discretization information, material properties, maximum load, essencial boundary conditions and the number of increments and iterations of the algorithm. With this information, the nodal mesh is constructed, along with the integration mesh. The influence domains are defined and the interpolation functions are calculated. With the material properties and the nodal connectivity, the initial stiffness matrix can be determined. With the natural and essential boundary conditions, the final stiffness matrix K K K0is obtained. After this pre-processing phase, the KT0 algorithm can start. The incremental load f f fiis defined as, f f fi=F F F inc (4.1) where F F Fis the maximum load and inc is the number of increments. The displacement field can be obtained by, u u ui=K K K−1 0f f fi(4.2) and consequently the stress field, σ σ σi=c c cB B Bu u ui(4.3) 24 Nonlinear Solution Algorithms The stress state is verified for every Gauss point. If the stress passes the yield surface, the return algorithm is applied. In the end of this process, the actualized stress vector of the increment, ∆ σ σ σi=σ σ σi−1+∆σ σ σ(4.4) where σ σ σi−1is the total stress vector of the previous increment. If any Gauss point passed the yield surface in the present increment or iteration, the residual forces must be calculated, using the following equation, f f fres i=f f fi−Z Ω B B BT∆σ σ σdΩ(4.5) The value of the residual forces is zero if no Gauss point passed the yield surface. In order to check whether the residual forces are significant or not, the following equation is calculated,   n ∑ m=1f f fres i m n  j ∑ k=1  n ∑ m=1f f fres i m n  <toler (4.6) where nis the number of nodes and jis the number of iterations. This serves as a stop criterion, where toler is the value introduced to define how significant the residual forces can be. If the equation is verified, the algorithm goes to the next load increment. Otherwise, the residual forces are still too big, and they must be applied to the structure. To do this, the process described in this section is repeated, replacing the incremental force defined in equation 4.1 for the residual force. When all the load increments are applied, the algorithm stops, and the load-displacement curves, as well as the stress and displacement fields can be obtained. 4.2 Mathematical Plasticity Since this work as the objective of analysing elastoplastic problems, it makes sense to present the concepts shown in chapter 3in a matricial form [42]. The elastic strain vector can be defined using Hooke’s law, dε ε ε=c c c−1dσ σ σ(4.7) Once the material elastic limit is achieved, the previous equation is no longer valid, since the total deformation now has a plastic and irreversible component. 25 Nonlinear Solution Algorithms 4.2.1 Uniaxial Yielding Considering an uniaxial test for a elastoplastic material, which generates the stress/strain curve presented in figure 3.3. It is visible that the material has an elastic behaviour, characterized by the elastic modulus, until the yield stress is achieved. After this point, the stress/strain relation is given by the tangent modulus. Having into account figure 3.3, the hardening parameter H0can be defined by, H0(ε ε εp) = ∂¯ σ σ σ ∂¯ ε ε εp =dσ σ σ dε ε ε−dε ε εe =1 dε ε ε dσ σ σ−dε ε εe dσ σ σ =1 1 ET−1 E =ET 1−ET E (4.8) Therefore, H0can be determined by a simple uniaxial test. 4.2.2 Plastic Flow Rule The expression of the Von Mises yield criterion is again presented, F(σ σ σ,ε ε εp,κ) = f(σ σ σ,ε ε εp,κ)−σY(κ) = 0 (4.9) Differentiating the previous equation, the following equation is obtained, dF =∂f ∂σ σ σT dσ σ σ−∂σY ∂κ dκ=0 (4.10) which can be also expressed as, dF =a a aTdσ σ σ−Adλ=0 (4.11) where a a ais the flow vector, normal to the yield surface, dλis the plastic strain-rate multiplier and Ais defined as, A=1 dλ ∂σY ∂κ dκ(4.12) The relation between stress and strain variation can be expressed by, dσ σ σ=c c c·dε ε εe=c c c·dε ε ε−dλ·c c c·a a a(4.13) The plastic strain-rate multiplier dλcan be defined by applying equation 4.13 to equation 4.11, dλ=a a aTc c cdε ε ε a a aTc c ca a a+A(4.14) The plastic deformation can now be obtained by, ε ε εp p p=dλa a a(4.15) 26 Nonlinear Solution Algorithms Strain hardening is used in this work, and by applying equation 3.52, dκ=d¯ εp=r2 3dε ε εT pdε ε εp(4.16) Taking into account equation 3.47 and the definition of the flow vector a a a=∂f ∂σ σ σ, the previous equation becomes, dκ=dλr2 3a a aTa a a(4.17) The Aparameter can be obtained from equation 4.12, A=σY dk r2 3a a aTa a a=σY dk ¯a(4.18) A linear variation is considered for the yield stress, σY(κ=¯ ε ε εp) = σ0 Y+H0¯ ε ε εp(4.19) Therefore, the parameter Abecomes, A=H0¯a(4.20) When the Von Mises criterion is used, ¯a=1, and so, A=H0(4.21) 4.2.3 Stress/Strain Relation As mentioned previously, Hooke’s law is not valid once the plastic phase is achieved. Hence, the stress variation can be expressed as a function of the strain variation by substituting equation 4.14 in equation 4.13, dσ σ σ=c c cepdε ε ε(4.22) where c c cep is defined by, cep =c c c−c c ca a aa a aTc c c A+a a aTc c ca a a(4.23) Defining d d dD=c c ca a a, the previous equation becomes, c c cep =c c c−d d dDd d dT D A+d d dT Da a a(4.24) 27 Linear Applications 0 0.2 0.40.60.8 1 −100 −80 −60 −40 −20 0 y (m) σxx (Pa) RPIM Analytic (a) 0 0.2 0.40.60.8 1 −100 −80 −60 −40 −20 0 y (m) σxy (Pa) RPIM Analytic (b) Figure 5.4: Square Plate Under Parabolic Stress: (a) σxx Stress Along x=0; (b) σxy Stress Along x=L/2. The displacement field and stress field distributions for the complete domain are displayed in figure 5.1. Is is possible to observe that the variable fields produced are extremely smooth. −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 (a) u −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 (b) v 0.02 0.04 0.06 0.08 0.1 0.12 (c) d −100 −80 −60 −40 −20 0 20 40 60 80 100 (d) σxx −80 −60 −40 −20 0 20 40 60 80 100 (e) σyy −180 −160 −140 −120 −100 −80 −60 −40 −20 0 (f) σxy Square Plate Under Parabolic Stress: Displacement Field and Stress Field Distributions 34 Linear Applications 5.2 Cantilever Beam The second example is a cantilever beam, represented in Figure 5.6(a). The stress field applied is given by, σxx =σ0x2 L2−y2 D2 σyy =σ0 x2L2−2D2 L2·D2+y2 L2!(5.5) σxy =−σ02x.y L2 The geometry of the model, the boundary and loading conditions, as well as the mechanical properties are presented in figure 5.6(a). Both regular and irregular meshes were analysed. Figure 5.6(b) and figure 5.6(c) represent a regular and an irregular mesh of 561 nodes. x y L=2m D=1mσxy σxy σxx A E=1000Pa ν=0.3 σ0=10Pa (a) (b) (c) Figure 5.6: Cantilever Beam: (a) Material and Geometric Conditions; (b) Regular 561 Node Mesh; (c) Irregular 561 Node Mesh. The analytical displacement solution for this model is given by [44], u=σ0 E"x3 3L2−x·y2 D2−ν x3L2−2D2 3L2·D2+x.y2 L2!# v=σ0 E"x2.yL2−2D2 L2.D2+y3 3L2−νx2.y L2−y3 3D2#(5.6) 35 Linear Applications Figure 5.7 shows the convergence study for the vertical displacement of point A (represented in figure 5.6(a)). The horizontal displacement is not presented, because the analytical solution is 0, and the erros obtained with RPIM are of the order of computer precision (10−15). It is possible to observe that the converged result for the vertical displacement is very close to the analytical solution (vA=0.375), whether regular or irregular meshes are used. 101102103104 0.32 0.34 0.36 0.38 Nodes vA regular mesh irregular mesh Analytic Figure 5.7: Cantilever Beam: vADisplacement Figure 5.8(a) shows the results for the convergence study of the medium displacement error. The medium error was calculated using the same equation 5.3 that was used for the first example. It is visible that the influence of the type of mesh in the medium displacement error is negligible, and the converged error is very low (0.5%). The same logic was used to analyse the convergence of the medium stress error (equation 5.4). The results obtained with regular meshes are presented in figure 5.8(b). The converged medium stress error is low (5%), and the irregular meshes produce very similar results. 101102103104 10−2 10−1 Nodes Medium Displacement Error regular mesh irregular mesh (a) 101102103104 10−1 100 Nodes Medium Stress Error σxx σxy (b) Figure 5.8: The Cantilever Beam: (a) Medium Displacement Error; (b) Medium Stress Error. 36 Linear Applications 0 0.2 0.40.60.8 1 −50 0 50 y (m) σxx (Pa) RPIM Analytic (a) 0 0.2 0.40.60.8 1 −1 −0.5 0 0.5 1 y (m) σyy (Pa) RPIM Analytic (b) 0 0.2 0.40.60.8 1 0 5 10 15 y (m) σxy (Pa) RPIM Analytic (c) Figure 5.9: Cantilever Beam: (a) σxx Stress Along x=L/2; (b) σyy Stress Along x=L/2; (c) σxy Stress Along x=L/2. Following the line of thought of the first example, the next analysis refers to the stress distribution along some interest lines of the solid. The following distributions were obtained using a regular mesh of 129×65 =8385 nodes. Figure 5.9 shows the evolutio of the three stresses (σxx in 5.9(a),σyy in 5.9(b) and σxy in 5.9(c)) along the line x=L/2. Apart from a small oscillation in the σyy stress, the RPIM and the analytical solution are practically the same. The displacement field and stress field distributions for the complete domain are displayed in figure 5.10. It is visible that the variable fields produced are very smooth. 37 Linear Applications −0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1 (a) u 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 (b) v 0.05 0.1 0.15 0.2 0.25 0.3 0.35 (c) d −100 −80 −60 −40 −20 0 20 40 60 80 100 (d) σxx 0 2 4 6 8 10 12 14 (e) σxy Figure 5.10: Cantilever Beam: Displacement Field and Stress Field Distributions. 5.3 Square Plate with a Circular Hole The final elastic example is a square plate with a circular hole. Figure 5.11(a) represents the simplified model analysed, since there are two planes of symmetry. The stress field applied is given by [44], σxx =σ01−a2 r23 2cos(2θ)+cos(4θ)+3a4 2r4cos(4θ) σyy =σ0−a2 r21 2cos(2θ)−cos(4θ)−3a4 2r4cos(4θ)(5.7) σxy =σ0−a2 r21 2sin(2θ)+sin(4θ)+3a4 2r4sin(4θ) The boundary conditions and the mechanical properties are presented in figure 5.11(a), while figures 5.11(b) and 5.11(c) represent a regular and an irregular mesh of 276 nodes each. 38 Linear Applications x y L=1m D=1m r=0.25m E=1000Pa ν=0.3 σ0=100Pa σxy σxy σxx σyy A B C D E F (a) (b) (c) Figure 5.11: Square Plate with a Circular Hole: (a) Material and Geometric Conditions; (b) Regular 276 Node Mesh; (c) Irregular 276 Node Mesh. The exact solution for the displacement is given by, u=σ0(1+¯ ν) ¯ E(1−¯ ν)r+2a2 rcosθ+a2 2r−a4 2r3cos(3θ) v=σ0(1+¯ ν) ¯ E−¯ νr−(1−2¯ ν)a2 rsinθ+a2 2r−a4 2r3sin(3θ)(5.8) where ¯ νand ¯ Eare defined as, ¯ ν=ν(1+ν) ¯ E=E1−¯ ν2(5.9) Figure 5.12 shows the convergence results for the displacement values of points A to F (represented in figure 5.11(a)). By analysing figure 5.12, it is possible to observe that the converged results are close to the analytical solution, and in general, the type of mesh has low impact on the solution. The convergence study for the medium displacement error is presented in figure 5.13(a). It is visible that the converged errors are very low, although the difference between the use of regular or irregular meshes is slightly bigger when compared to the previous examples. Figure 5.13(b) shows the the results for the convergence study of the medium σxx stress error, using a regular mesh. The results for irregular meshes are similar to the ones presented, and the converged error is low (below 2%). 39 Linear Applications 101102103104 0.101 0.102 0.103 0.104 0.105 Nodes uA regular mesh irregular mesh Analytic (a) 101102103104 −3 −2.8 −2.6 −2.4·10−2 Nodes vA regular mesh irregular mesh Analytic (b) 101102103104 0.116 0.117 0.118 Nodes uB regular mesh irregular mesh Analytic (c) 101102103104 7.4 7.5 7.6 7.7 7.8·10−2 Nodes uC regular mesh irregular mesh Analytic (d) 101102103104 5.2 5.25 5.3 5.35 5.4 5.45 ·10−2 Nodes uD regular mesh irregular mesh Analytic (e) 101102103104 −1.8 −1.7 −1.6 −1.5 ·10−2 Nodes vD regular mesh irregular mesh Analytic (f) 101102103104 −2.4 −2.2 −2 −1.8·10−2 Nodes vE regular mesh irregular mesh Analytic (g) 101102103104 −3.86 −3.84 −3.82 −3.8 −3.78 ·10−2 Nodes vF regular mesh irregular mesh Analytic (h) Figure 5.12: Square Plate with a Circular Hole: (a) uADisplacement; (b) vADisplacement; (c) uB Displacement; (d) uCDisplacement; (e) uDDisplacement; (f) vDDisplacement; (g) vEDisplacement; (h) vFDisplacement. 40 Linear Applications 102103104 10−3 10−2 10−1 Nodes Medium Displacement Error regular mesh irregular mesh (a) 102103104 10−2 10−1 Nodes Medium Stress Error σxx (b) Figure 5.13: Square Plate with a Circular Hole: (a) Medium Displacement Error; (b) Medium Stress Error. The study of the stress distribution along interest lines of the solid focused on the two edges with essencial boundary conditions. Therefore, figures 5.14(a) and 5.14(b) represent the evolution of the σxx and σyy stresses along the line x=0, while figures 5.14(c) and 5.14(d) represent the evolution of the same stresses along the line y=0. It is clear that the RPIM solution replicates the analytical solution with great precision. Figure 5.15 displays the displacement and stress fields across the whole domain. Even though the displacement fields are very smooth, the stress fields display small zones where the stress variations are slightly more abrupt. 0 0.1 0.2 0.3 0.40.5 0.60.7 0.8 0.9 1 100 150 200 250 300 y (m) σxx (Pa) RPIM Analytic (a) 0 0.1 0.2 0.3 0.40.5 0.60.7 0.8 0.9 1 0 10 20 30 40 y (m) σyy (Pa) RPIM Analytic (b) 0 0.1 0.2 0.3 0.40.5 0.60.7 0.8 0.9 1 0 20 40 60 80 x (m) σxx (Pa) RPIM Analytic (c) 0 0.1 0.2 0.3 0.40.5 0.60.7 0.8 0.9 1 −100 −50 0 x (m) σyy (Pa) RPIM Analytic (d) Figure 5.14: Square Plate with a Circular Hole: (a) σxx Stress Along x=0; (b) σyy Stress Along x=0; (c) σxx Stress Along y=0; (d) σxx Stress Along y=0. 41 Linear Applications 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 (a) u −0.035 −0.03 −0.025 −0.02 −0.015 −0.01 −0.005 0 (b) v 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 (c) d 0 50 100 150 200 250 300 (d) σxx −80 −60 −40 −20 0 20 40 (e) σyy −80 −70 −60 −50 −40 −30 −20 −10 0 10 (f) σxy Figure 5.15: Square Plate with a Circular Hole: Displacement Field and Stress Field Distributions. 42 Chapter 6 Non Linear Applications In this chapter several two dimensonal examples are analysed, in order to assess the performance of the nonlinear algorythm (KT0) with the RPIM. First, two well-known examples (cantilever beam and cook’s membrane) are studied using a monotonic load. After that, the cantilever beam is analysed with two different kinds of cyclic load. The results obtained with the RPIM are compared with FEM solutions (either ABAQUS or ANSYS). The last example is the study of the trabecular bone and a polymer analogue. 6.1 Cantilever Beam The cantilever beam problem represented in figure 6.1(a) is analysed. The domain was discretized in a regular mesh of 561 nodes, represented in figure 6.1(b). The dimensions and mechanical properties are also presented in figure 6.1(a). y x L=8m D=4mAP E=3×107Pa ET=3×106Pa σy=3×104Pa ν=0.3 (a) (b) Figure 6.1: Cantilever Beam: (a) Material and Geometric Conditions; (b) Regular 561 Node Mesh. Figure 6.2 presents the vertical displacement evolution of point A (represented in figure 6.1(a)) with the increasing load P. It is possible to observe that the RPIM fits the ANSYS solution (obtained from [31]) perfectly. 43 Non Linear Applications 6.4 Trabecular Bone In this section the trabecular bone is analysed. Figure 6.11 shows the experimental results obtained by [45] for the confined square test. The geometry, load and boundary conditions used are presented in figure 6.12. By observing the experimental data, it is visible that the trabecular bone’s mechanical behaviour can be represented by a bilinear material model, at least until densification occurs, as shown in figure 6.11. 0 0.1 0.2 0.3 0.40.5 0.6 0 20 40 60 Strain P(MPa) Experimental Bilinear Figure 6.11: Trabecular Bone: Experimental Results and Bilinear Model The mechanical properties of the bilinear model are presented in figure 6.12 (the Poisson’s ratio is obtained from the literature [45]). For the RPIM analysis, a 289 node mesh was used, similar to the one represented in figure 5.1(a). x y L=1cm D=1cm P ATrabecular Bone E=350.5MPa ET=18.6MPa σY=13.4MPa ν=0.16 Polymer Analogue E=147.6MPa ET=3.6MPa σY=5.6MPa ν=0.20 Figure 6.12: Trabecular Bone and Polymer Analogue: Material and Geometry Conditions The strain of point Ais presented in figure 6.13. It is visible that, even though the RPIM presents a slightly stiffer behaviour than the bilinear model, the relation between the two is fairly good. 50 Non Linear Applications 05·10−20.10.15 0.20.25 0.30.35 0.4 0 10 20 30 Strain P(MPa) RPIM Experimental Bilinear Figure 6.13: Trabecular Bone: εAStrain A trabecular bone analogue material is also studied. Figure 6.14 shows the experimental results obtained for the celular rigid polyurethane (PU) foam. The same boundary and load conditions are applied. Once again, a bilinear model is defined to simulate the foam material until densification. The material properties for the bilinear model are presented in figure 6.12. 0 0.1 0.2 0.3 0.40.5 0.6 0 5 10 Strain P(MPa) Experimental Bilinear Figure 6.14: Polymer Foam: Experimental Results and Bilinear Model Figure 6.15 shows the strain of point Afor the PU foam analysis. The results are very similar to the trabecular bone study. Although the RPIM solution is slightly stiffer, the relation between the RPIM solution and the bilinear model is fairly good. 51 Non Linear Applications 05·10−20.10.15 0.20.25 0.30.35 0.4 0 2 4 6 8 10 Strain P(MPa) RPIM PU Model Figure 6.15: Polymer Foam: εAStrain 6.5 Cortical Bone The final elastoplastic example is the study of the trabecular bone. Figure 6.17 shows the bilinear models of the cortical bone. These were obtained through the analysis of experimental tests [46]. The tests were performed at different speed levels, and figure 6.17 presents the ones performed with the highest and lowest speeds. The geometry, boundary and load conditions, as well as the mechanical properties of the bilinear models are presented in figure 6.16. The model has the same dimensions as the model for the trabecular bone, but in this case, the base is clamped and the sides are free. x y L=1cm D=1cm P AHigh Speed E=18.5GPa ET=2.2GPa σY=128.9MPa ν=0.36 Low Speed E=11.6GPa ET=0.95GPa σY=79.3MPa ν=0.36 Figure 6.16: Cortical Bone: Material and Geometry Conditions The RPIM solution for the strain of point Ais presented in figure 6.17. It is visible that the relation between the RPIM solution and the bilinear model is very close, for both the slow and high speed models. 52 Non Linear Applications 01·10−22·10−2 0 50 100 150 200 Strain P(MPa) RPIM Low Speed High Speed Figure 6.17: Cortical Bone: εAStrain for High and Low Speed 53 Non Linear Applications 54 Chapter 7 Conclusions and Future Works In this work, the RPIM was applied to the study of non-linear problems, with both monotonic and non monotonic loads. First, several linear elastic examples were analysed. With the results obtained from these tests, it is possible to conclude that the RPIM solutions are very close to the exact problem solutions, even with relatively small meshes. The displacement and stress fields obtained are accurate and very smooth. Also, the use of regular or irregular meshes has little or no influence in the final solution. The RPIM was then applied to elastoplastic problems. Regarding the monotonic load cases, the results obtained for the displacement are extremely close to the FEM solution, with an equal or smaller mesh. The stress results are also very close to the FEM solution, however in some cases the relation is slighty different. From the analysis of the non monotonic load cases, it was concluded that, even though the model yields a little sooner than it is supposed to when the structure is reloaded in the same direction, the RPIM solution is very close to the ABAQUS solution. The final application was the trabecular and cortical bone study. For the trabecular bone (and polymer analogue), the RPIM solution is slightly stiffer than the bilinear model, although the relation between the two is fairly good. The results obtained for the cortical bone are very close to the bilinear model, for both the high speed and low speed cases. In general, the elastoplastic results show that the non-linear solution algorithm, as well as the procedure for the stress returning were successfully implemented in the RPIM. 55 Conclusions and Future Works Even though the RPIM has good accuracy and produces smooth stress fields, its main disadvantage when compared to the FEM is the high computational cost. Several additions could enhance this work, such as: •The inclusion of large deformations; •The inclusion of dynamic analysis; •The inclusion of new hardening rules, namely the kinematic hardening rule, in order to reproduce the Bauschinger effect with non monotonic loads; •The study of other modified Newton-Raphson methods (for example KT1), as well as the study of new non-linear algorithms, such as the arc-length method, to reproduce structural instabilities, namely the "snap-through" and the "snap-back" phenomena; 56 References [1] B. Nayroles, G. Touzot, and P. Villon, “Generalizing the finite element method: Diffuse approximation and diffuse elements,” Computational Mechanics, vol. 10, pp. 307–318, 1992. [2] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, “Meshless methods: an overview and recent developments,” Computer Methods in Applied Mechanics and Engineering, vol. 139, pp. 3–47, 1996. [3] P. W. Randles and L. D. Libersky, “Smoothed particle hydrodynamics: Some recent improvements and applications,” Computer Methods in Applied Mechanics and Engineering, vol. 139, pp. 375–408, 1996. [4] V. P. Nguyen, T. Rabczuk, S. Bordas, , and M. Duflot, “Meshless methods: A review and computer implementation aspects,” Mathematics and Computers in Simulation, vol. 79, pp. 763–813, 2008. [5] J. J. Monaghan, “Smoothed particle hydrodynamics: Theory and applications to nonspherical stars,” Monthly Notices of the Astronomical Society, vol. 181, pp. 375–389, 1977. [6] L. D. Libersky and A. G. Petschek, “Smoothed particle hydro-dynamics with strength of materials,” in The Next Free Lagrange Conference, pp. 248–257, 1991. [7] T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free galerkin methods,” International Journal for Numeric Methods in Engineering, vol. 37, pp. 229–256, 1994. [8] P. Lancaster and K. Salkauskas, “Surfaces generation by moving least squares methods,” Mathematics of Computation, vol. 37, pp. 141–158, 1981. [9] W. K. Liu, S. Jun, S. Li, J. Adee, and T. Belytschko, “Reproducing kernel particle methods for structural dynamics,” International Journal for Numeric Methods in Engineering, vol. 38, pp. 1655–1679, 1995. [10] Z. T. Atluri, “A new meshless local petrov–galerkin (mlpg) approach in computational mechanics,” Computational Mechanics, vol. 22(2), pp. 117–127, 1998. [11] L. M. J. S. Dinis, R. M. N. Jorge, and J. Belinha, “Analysis of 3d solids using the natural neighbour radial point interpolation method,” Computer Methods in Applied Mechanics and Engineering, vol. 196, pp. 2009–2028, 2007. [12] K. M. Liew, X. Zhao, and A. J. M. Ferreira, “A review of meshless methods for laminated and functionally graded plates and shells,” Composite Structures, vol. 93, pp. 2031–2041, 2011. 57 REFERENCES [13] G. Liu and Y. Gu, “A point interpolation method for two-dimensional solids,” International Journal for Numeric Methods in Engineering, vol. 50, pp. 937–951, 2001. [14] J. Wang, G. Liu, and Y. Wu, “A point interpolation method for simulating dissipation process of consolidation,” Computational Methods in Applied Mechanics and Engineering, vol. 190, pp. 5907–5922, 2001. [15] J. Wang and G. Liu, “A point interpolation meshless method based on radial basis functions,” International Journal for Numeric Methods in Engineering, vol. 54, pp. 1623–1648, 2002. [16] J. Wang and G. Liu, “On the optimal shape parameters of radial basis functions used for 2-d meshless methods,” Computational Methods in Applied Mechanics and Engineering, vol. 191, pp. 2611–2630, 2002. [17] J. Braun and M. Sambridge, “A numerical method for solving partial differential equations on highly irregular evolving grids,” Nature, vol. 376, pp. 655–660, 1995. [18] N. Sukumar, B. Moran, A. Semenov, and V. Belikov, “Natural neighbour galerkin methods,” International Journal for Numeric Methods in Engineering, vol. 50, pp. 1–27, 2001. [19] R. Sergio, S. Idelsohn, E. Onate, N. Calvo, and F. D. Pin, “The meshless finite element method,” International Journal for Numeric Methods in Engineering, vol. 58, pp. 893–912, 2003. [20] J. Belinha, L. M. J. S. Dinis, and R. M. N. Jorge, “The natural radial element method,” International Journal for Numerical Methods in Engineering, vol. 93, pp. 1286–1313, 2013. [21] J. Belinha, L. M. J. S. Dinis, and R. M. N. Jorge, “Composite laminated plate analysis using the natural radial element method,” Composite Structures, vol. 103, pp. 50–67, 2013. [22] J. Belinha, “The radial point interpolation meshless method introduction course.” Institute of Mechanical Engineering - FEUP, September 2011. [23] Y. T. Gu, “Meshfree methods and their comparisons,” International Journal of Computational Methods, vol. 2(4), pp. 477–515, 2005. [24] O. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, 1989. [25] K. J. Bathe, Finite Element Procedures. Prentice-Hall: Englewood Cliffs, 1996. [26] E. J. Kansa, “A scattered data approximation scheme with applications to computational fluid-dynamics,” Computers and Mathematics with Applications, vol. 19, pp. 127–161, 1990. [27] J. Belinha, “Elasto-plastic analysis considering the element free galerkin method,” Master’s thesis, Faculdade de Engenharia da Universidade do Porto, 2004. [28] G. R. Liu, “A point assembly method for stress analysis for two-dimensional solids,” International Journal of Solid and Structures, vol. 39, pp. 261–276, 2002. [29] R. L. Hardy, “Theory and applications of the multiquadrics - biharmonic method (20 years of discovery 1968-1988),” Computers and Mathematics with Applications, vol. 19, pp. 163– 208, 1990. 58 REFERENCES [30] J. Belinha, R. M. N. Jorge, and L. M. J. S. Dinis, “Bone tissue remodelling analysis considering a radial point interpolator meshless method,” Engineering Analysis with Boundary Elements, vol. 36, pp. 1660–1670, 2012. [31] S. F. Moreira, “Elastoplastic analysis using the natural neighbour radial point interpolation method,” Master’s thesis, Faculdade de Engenharia da Universidade do Porto, 2013. [32] S. Moreira, J. Belinha, L. M. J. S. Dinis, and R. M. N. Jorge, “Analysis of laminated beams using the natural neighbour radial point interpolation method (nnrpim),” Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, p. in press, 2013. [33] G. R. Liu and S. S. Quek, The Finite Element Method: A Pratical Course. ButterworthHeinemann, 2003. [34] S. Timoshenko and J. N. Goodier, Theory of Elasticity. McGraw-Hill, 1951. [35] J. F. S. Gomes, Mecânica dos Sólidos e Resistência dos Materiais. Porto: Edições INEGI, 2009. [36] G. R. Liu, Meshfree Methods: Moving Beyond the Finite Element Method. CRC Press, New York, 2000. [37] J. N. Reddy, Applied Functional Analysis and Variational Methods in Engineering. McGrawHill International Edition, 1986. [38] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method for Solid and Structural Mechanics. Elsevier Butterworth-Heinemann, 2005. [39] J. A. O. P. Belinha, The Natural Neighbour Radial Point Interpolation Method. PhD thesis, Faculdade de Engenharia da Universidade do Porto, 2010. [40] J. Chakrabarty, Theory of Plasticity. McGraw-Hill International Edition, 1987. [41] R. Hill, The Mathematical Theory of Plasticity. Oxford University Press, 1950. [42] D. R. J. Owen and E. Hinton, Finite Elements in Plasticity: Theory and Practice. Pineridge Press Limited, 1980. [43] M. A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures. John Wiley Sons, 2000. [44] J. Belinha, L. M. J. S. Dinis, and R. M. N. Jorge, “The natural radial element method,” International Journal for Numerical Methods in Engineering, p. Published online in Wiley Online Library, 2012. [45] N. Kelly and J. P. McGarry, “Experimental and numerical characterisation of the elastoplastic properties of bovine trabecular bone and a trabecular bone analogue,” Journal of the Mechanical Behavior of Biomedical Materials, vol. 9, pp. 184–197, 2012. [46] A. N. Natali, E. L. Carniel, and P. G. Pavan, “Constitutive modelling of inelastic behaviour of cortical bone,” Medical Engineering Physics, vol. 30, pp. 905–912, 2008. 59