PSI - Issue 33

I.J. Sánchez-Arce et al. / Procedia Structural Integrity 33 (2021) 149–158 Sánchez-Arce et al. / Structural Integrity Procedia 00 (2019) 000–000

153

5

respect to each principal stretch was calculated, giving the Second Piola-Kirchhoff stress tensor ( S ). Then, the Kirchhoff stress (  ), or engineering stress, is T  FSF  , which divided by J ( det( ) J  F ), gives the Cauchy stress (  ). For the Ogden model, assuming incompressibility and a plane stress case,  would be:

N

OGD  τ i 

 1 2          p p p  1

 

; 1, 2 , i 

(12)



1

p



while for plane strain,  would be:





N

1 3





  

  

p

OGD 

  ln ;









τ

1 ,

(13)

J

K J



p 

i 



p   

p    p

3 



3

p

1

2

3

i

1

p



being K the material’s bulk modulus. Equations (12) and (13) were implemented into the MATLAB scripts corresponding to the meshless method. Therefore, for each increment of the non-linear solution, the updated stress state was calculated. Then, the non-linear solution algorithm was that described by Farahani et al., (2019) for elastic plastic analyses. 4. Validation of the methodology The first testing case corresponds to a square plate with a side size of 2 m subject to a displacement of 0.5 m on its right vertical edge, symmetry conditions around the XZ plane were considered (Figure 1). The material was considered as a Mooney Rivlin with C 10 = 349 kN/m 2 ; C 01 = 87 kN/m 2 , and D 1 = 2.2935 x 10 7 m 2 /N, as described by Khosrowpour et al., (2019); the relationship between D 1 and K is 1 2 K D  . The geometry was divided into an array of 21 x 9 nodes (20 x 8 divisions) along the X and Y dimensions, respectively (Figure 1b). This nodal distribution, pre-processed with the NNRPIM, gave a distribution of 1280 integration points. The left vertical edge was constrained in both X and Y directions (Ux = Uy = 0). The lower horizontal edge was only constrained in the vertical direction to represent the symmetry condition (Uy = 0). Finally, the right vertical edge was subject to a displacement (Ux = 0.5 m) while left free to displace in the vertical direction; thus, shrinking of this edge is expected (Figure 1a). The problem was modelled as a plane stress case, as suggested by Khosrowpour et al., (2019). The models were run using custom-written MATLAB scripts, and so the initial Young’s modulus, E 0 , and shear modulus, G 0 , were calculated from the Mooney-Rivlin material constants (Equation (7));  was calculated using G 0 and D 1 . Although this procedure is not necessary for a mono-material hyper-elastic analysis, the scripts had been written to consider multi-material models involving elastic-plastic and hyper-elastic behaviours, and so the method can be expanded for further analyses. In addition, to simulate the Mooney-Rivlin behaviour with the Ogden model, N ODG =2. The parameters  1 and  2 for the Ogden model were calculated from the recommendations from de Souza Neto et al., (1995), while  1 and  2 were 2 and -2, respectively (de Souza Neto et al., 1995). The methodology was also validated by comparing the solution from NNRPIM with a solution from the ABAQUS software (ABAQUS 6.17. Dassault Systèmes, Inc); however, the nodal distribution in the ABAQUS model was finer (40 x 16 divisions) to reduce stress concentrations due to mesh. The reaction force at the left vertical edge (Figure 1) was calculated as the sum of the horizontal forces on all the nodes along that edge. The displacement corresponds to the horizontal translation of the right vertical edge’s lowest node. This procedure was performed with both NNRPIM and FEM (using ABAQUS’) solutions. Normal and shear stresses were determined at x=1, vertical displacement of the upper node at the higher right node was also tracked through the solution.