Numeric Computational Scheme for Adaptive Tearing Model

Numeric Computational Scheme for Adaptive Tearing Model

Astris Dyah Perwita

Informatics Engineering Master Program Institut Teknologi Sepuluh Nopember Surabaya Indonesia

Rully Soelaiman, Darlis Herumurti Department of Informatics Engineering Institut Teknologi Sepuluh Nopember Surabaya, Indonesia

AbstractTo produce details and component of tearing and cracking is complicated. The detail and component for every object would also produce certain type of tearing line. Meanwhile using traditional finite element method will need such a complex computation caused by updated value for every crack made. Such model needed to predict magnitude and angle of tearing or cracking depend on the variables value.

To generate tearing model that satisfy the requirements and produce tearing line with small computational cost, numerical scheme will be needed. This research will be focused on thin plates especially paper. Modification of finite element method will be proposed. Every step of process will be done adaptively depending on the existing results each step.

Experiments generate a model which stiffness matrix convergent for each time function and each variation of elastic variables such as Young modulus and Poissons ratio. Experiments shows that computation time for tearing model without inclusion component works on average 21.17 seconds faster than with inclusion with average 30.57 seconds. Computation time for tearing model with circular inclusion work faster dan linear inclusion.

Keywords Graphic, Stiffness Matrix, Tearing Propagation, FEM Modification

1. INTRODUCTION

   Dynamic tearing for various objects has been developed for years in computer graphic field. However, to produce details and tearing variation for various objects still considered challenging [1]. Meanwhile using traditional finite element method for tearing and cracking will need a complex computation. Therefore, this research will proposed on simulation technique which focused on applying proper method and components values with modification of finite element method so that the model produce tearing line, magnitude, and angle precisely [2].

   Research by Terzopoulos and Fleischer represent viscoelastic and plastic deformation model. This method applied with three fundamental metric tensors to represent changes in energy function which calculate the deformation on face and volume of object. The behaviors are viscoelasticity, plasticity, and fracture. The deformation will be done locally based on viscous and plastic processes within the models. The simple fracture mechanic process depends on yield and creep relationship that be affected by force and/or instantaneous deformation. This fracture introduces local discontinuities as a function of the instantaneous deformation through the model [3].

   When objects deform by forces, thin sheets of material such as paper tend to crumple. A crumple defines as distinctive patterns characterized by networks of sharp folds and cone singularities. These crumples form due to interaction between low bending resistance and high in-plane stiffness. Framework by Narain allows efficient modeling sharp features and avoids bend locking that would be otherwise caused by stiff in-plane behavior. Shape diffusion caused by remeshing prevented by using an explicit plastic embedding space [4].

   In particular case, fractures caused by tearing are different form fractures caused by in-plane motions. Stress relaxation method applied to avoid shattering artifacts after generating each fracture cut. To adaptively provide some fracture details, a fracture-aware remeshing scheme based on constrained Delaunay triangulation formulated. The method will provide efficient and realistic multi-layered fractures of thin plates [5].

   Displacement-based approximation caused by discontinuities with the finite element method is chalenging due to the need for update the mesh to match the discontinuities of the geometry. Some research provide new technique to model cracks and crack growth without remeshing [6].

   Deformation technique without remeshing applied by using the unity properties of finite element method. This concept will allows local enrichment elements to be combined into finite element approach. This approach aplied in such a way that degree of freedom from enrichment elements connect with the mesh. Enrichment elements on tearing and cracking will be focused on crack tip area and calculate discontinuities function to represent displacement along the tearing lines [7].

   To model strong discontinuities which represent tearing area, two types of enrichment elements needed. The first element is Heaviside elements which represent discontinuities along tearing line. The second element is crack tip element which considered as determinant of the next crack. Tearing line represented by mesh independently by enrichment element so that the geometry of tearing could be updated without adding or rebuild the mesh in the domain. Material interface represent by enrichment element will be processed using absolute value enrichment which gives the shortest signed distance from a given point to the interface between two materials [8].

   To produce accurate and efficient computational stages, such numerical schemes which satisfy the requirements will be needed. Based on previous research, we proposed tearing models with specific value of variables so that the computational process works realistically and efficiently

   against displacement by the time function. Approach on

   (), 0 0

   parameters value expected to optimize and reduce

   computational cost.

   0 ( ), 0

   0 0 (),

2. THEORETICAL BACKGROUND

   =

   (6)

   0 ( ), (),

   1. XFEM

      (), 0 ,

      Finite Element method used to solve problem related to

      [(), ( ), 0 ]

      discontinuity. Discontinuity will happened around the tearing tip nodes and will cause displacement of the nodes position on domain. The discontinuity value and displacement will be calculated with Finite Element method. The point change on

3. METHOD

   domain is calculated using (1).

   () = () [ + ()]

   (1)

   Tearing model performed on Matlab with certain required

   stages of method. The method stages explained as pseudocode on Fig. 1.

   Input consists of five main elements. The five elements are

   Domain which represents the objects size, MAT which represents objects variables, Crack which represents the initial

   Variable in Equation 1 is domain while is the discontinued domain. () is traditional finite elements,

   () is enrichment elements discontinuities function, while and is degree of freedom for traditional and enrichment elements.

   Inclusion elements represented as formulation using (2).

   () = ( )2+( )2 (2)

   Variable and are the i-th coordinate meanwhile and are center coordinates for composite material, and is the radius of the inclusion.

   On discrete formulation of finite element method, Heaviside and tearing tip elements approach represents as (3) and (4)

   crack of the object, INC which represents the inclusion elements, Grow which represents the amount of tearing that must be established, and Force which represents the amount of force that must be distributed on the object. These inputs proceeded in stages as follows

   earing Model Pseudocode

   1. Stiffness matrix

   ((), ) 0 (3)

   1. Input

   2. If isempty(GROW) = 1

   3. iter 1

   4. elseif length(GROW) = 2

   5. iter GROW(1)

   6. end

   7. for i = 1 to iter

   8. if i = 1

   9. connectivity

   10. pHDOF []

   11. else

   12. clear status on NODES

   13. end

   14. omega levelSet(i)

   15. [DOF,DISP] calcDOF

   16. [updElem, IElem] = enrElem(i,pHDOF)

   17. if i = 1

   18. globalK stiffnessMatrix

   19. else

   20. globalK updateStiffness

   21. end

   22. globalF forceVector

   23. freeDOF boundaryCond

   24. DISP(freeDOF,:) globalK \ globalF

   25. if isempty(CRACK) != 0

   26. pHDOF 2*max(NODES(:,2))

   27. [KI,KII] JIntegral(omega,DISP) 28. Gopt -(KI(1)^2+KII(1)^2)*(1-

   MAT(2)^2)/MAT(1);

   1. exit growCrack(KI,KII,omega)

   2. if exit = YES

   3. plot result

   4. break

   5. end

   6. print time

   7. end

   8. end

   ((), ) = 0 (4)

   Displacement calculated by (1) happen numeric computation field. To apply the displacement into mesh domain and interface domain, the displacement value needs to be adapted into its nodal form. Interface domain could be represented as matrix. This matrix also include stress and strain values which declare force within domain.

   To obtain stress and strain value, stiffness matrix which represent derivative of initial shape function needs to be calculated. For traditional stiffness matrix without enrichment elements, the matrix could be represents as (5). And for stiffness matrix for enrichment elements could be represents as (6)

   , 0 0

   0 , 0

   0 0

   =

   ,

   (5)

   0 , ,

   , 0 ,

   [, , 0 ]

   Figure 1. Tearing Model Pseudocode

   1. Connectiviy Formation

      Connectivity is formed by Domain input. The size that must be established is expressed in the form of a grid and represented as a matrix of connectivity. Connectivity matrix stores coordinates of the nodes, the index element and its members, as well as index element enrichment information at any nodes. This stage only executed at initialization or first iteration. For the next iteration neutralizing index element enrichment information by returning the value to null is performed.

   2. Level Set Formation

      Level Set used as part of discretization method. Level Set is a function that used for stating the shape function of elements that exist in the domain so that the value of each element may be calculated as information on the connectivity matrix.

      For each node, the calculated value (x) as a shape function of inclusions elements in accordance with Equation

      2. Declaration of initial tearing tip node on the matrix domain is calculated. In accordance with Equation 3 and 4, the points are eligible to be categorized as Heaviside elements and endpoints tear.

   3. Degree of Freedom Calculation

      Degree of Freedom is calculated in corresponding to each node in each element. For Heaviside element calculated in accordance with the terms Equation 7 and tip tear element is calculated by Equation 8.

4. RESULT AND ANALYSIS

The main experiment was conducted with two scenarios, the first one is without inclusion and the second one is using inclusion. From both of the main scenarios, there are several sub scenarios that are affected by Youngs modulus variable and Poissons ratio based on user input. In the scenario with inclusion, there are two scenarios that are determined by the inclusion that are appliedcircular and linear.

The results of the experiment without inclusion are shown in the form of nodal stress values as shown in Fig. 2. Meanwhile, the results of the trial with inclusion are shown in the contour of nodal stress values as shown in Fig. 3. Nodal stress value is a stress values that are influenced by Gausss value at a point. In this contour image shown that every point tear is depended on stiffness matrix computed at any time function. The contour image also shows that the distribution of the stiffness matrix at each iteration are still associated with the stress distribution matrix in the previous iteration.

The experiment result in terms of stiffness matrix assembly model without inclusion is shown in Fig 4. While the experiment results in terms of stiffness matrix assembly model with the inclusion is shown in Fig 5. Both of the implementation graphs of stiffness matrix assembly against the time function indicates the properties of convergence model against matrix stiffness calculation for each iteration. The experiment results show that the model began to converge on average at fifth iteration. The trial results also showed that the convergence rate will be applicable on both

= {

( ) 1 , above tearing line

1 , below tearing line

(7)

models with or without the inclusion as well as with the variation of elasticity variable.

=14 [

(), = , , , ] (8)

2 2 2 2

1. Stiffness Matrix Calculation

   The Stifness matrix is calculated by calculate and as in (5) and (6). Matrix value of and will later be a global stiffness matrix which used in the next step to calculate stress and strain value.

   The updated stiffness matrix for the second iteration until the appointed total iteration will be done without calculating inclusion nodes matrix because it doesnt altered by time function.

2. Force Calculation

   The force distributed in nodes is in corresponding with Degree of Freedom which calculated before. Force magnitude value obtained from users input.

3. Tearing Tip Determination

   The first thing to do to determine the second and further tearing is calculate stress intensity factor by J-Integral method. Stress intensity factor is used to determine tearing direction. Further tearing coordinates determined by force from users input and the angle which calculated from stress intensity factor.

   From the experiment results obtained an average computation time that are required to form a tear in the object. The computing time can be seen entirely in Table 1. It can be seen that in table 1 the form of a tear in material without inclusion is faster than the one with inclusion. Meanwhile the time required for the circle inclusion is faster than linear inclusion.

   Table 1. Experiments Computational Time
   No	Modulus	Rasio	Inclusion	Computational
   Young Poisson time (seconds)
   1	10E9	0.3	–	23.67
   2	5E9	0.3	–	21.34
   3	10E9	0.1	–	24.49
   4	10E9	0.3	Circular	24.51
   5	5E9	0.3	Circular	25.29
   6	10E9	0.1	Circular	24.46
   7	10E9	0.3	Linear	33.68
   8	5E9	0.3	Linear	36.44
   9	10E9	0.3	Linear	39.05

   (a) (b) (c) (d)

   Figure 2. Stress distribution for tearing model without inclusion (a) 1st iteration (b) 10th iteration (c) 20th iteration (d) 3th iteration

   (a) (b) (c) (d)

   Figure 3 Stress distribution for tearing model with inclusion (a) 1st iteration (b) 10th iteration (c) 20th iteration (d) 30th iteration

   Figure 4. Stiffness matrix assembly time without inclusion model

   Figure 5. Stiffness matrix assembly time with inclusion model

   REFERENCES

   1. T. Pfaff, R. Narain, J. M. de Joya and J. E. O'brien, "Adaptive Tearing and Cracking of Thin Sheets," ACM Transaction on Graphics, vol. 33, 2014.

   2. M. J. Pais, N.-H. Kim and T. Davis, "Reanalysis of the Extended Finite Element Method for Crack Initiation Propagation," Florida, 2010.

   3. D. Terzopoulos and K. Fleischer, "Modeling Inelastic Deformation: Viscoelasticity, Plasticity, Fracture," Computer Graphics, vol. 22, no. 4, pp. 269-278, 1988.

   4. R. Narain, T. Pfaff and J. F. O'brien , "Folding and Crumpling Adaptive Sheets," California, 2013.

   5. O. Busaryev, T. K. Dey and H. Wang, "Adaptive fracture simulation of multi-layered thin plates," ACM Trans. Graph, vol. 4, no. 32, 2013.

   6. N. Moes, J. Dolbow and T. Belytschko, "A Finite Element Method For Crack Growth Without Remeshing," International Journal For Numerical Methods in Engineering, vol. 46, pp. 131-150, 1999.

   7. J. M. Melenk and I. Babuska, "The partition of unity finite element method: Basic theory and applications," Computer methods in applied mechanics and engineering, vol. 139, pp. 289-314, 1996.

   8. N. Sukumar, D. Chopp and T. Belytschko, "Modeling holes and inclusions by level sets in the extended finite-element method," Computer methods in applied mechanics and engineering, vol. 190, pp. 6183-6200, 2001.