Numerical simulations of fabric shields for bulletproof vest: determining the number of sheets to withstand ballistic impacts

The main contribution of this research is to develop a FE model for evaluating impact loading against fabrics. Although the fabric model has been studied by some authors, this article contributes with a FE model with a new material formulation which takes into account the stress-strain limit to simulate the fabric rupture. As a case study, the author uses a bulletproof vest impact analysis and determine the number of sheets necessary to resist the ballistic impact. A sophisticated numerical algorithm is used. New constitutive equations are developed capable to take into account the material damage, yield stresses and rupture limits. This is based on elastoplasticity theory, with an adaptation of the function related to the hardening of the material. With this, it is possible to control the limits of the permissible stresses on the yarn and coating of the structural fabric and accurately represent the physical behavior of the materials used. The elastoplasticity and damage theories are developed and, finally, a numerical algorithm is presented for the determination of stresses, displacements and material deformations.


INTRODUCTION
The particle method can be used in several fields, especially into researches about structural fabrics. It consists in studying the armor systems modeling, protections against explosions and lightweight structural fabrics submitted to the impact of a projectile. These lines of research use the concepts of the theory of the particles, see Zohdi and Powell (2006), Zohdi (2003), Zohdi and Wriggers (2001a) and Bandeira and Zohdi (2019).
Fabric shields have been studied by several researchers. Specially, because there exist a wide range of applications for lightweight ballistic fabric shields, such as the protection of critical structural components in transport systems and the human body. Important works in this field were presented by Zohdi and Powell (2006), Kato, Yoshino and Minami (1999), Pargana (2007) and Zohdi (2014c). In these researches, the authors affirm that some deficiencies are the susceptibility to being abruptly severed by sharp objects, which completely eliminates the fabric's ability to stretch and absorb incoming kinetic energy and the environmental degradation of the fabric due to moisture, heat and sunlight, which is of growing concern, since many new fabrics have multiple purposes, such as electrical and chemical sensing, in addition to being part of a protective system.
Basically, Zohdi developed a computational framework using a coated network model in order to capture the basic characteristics of such systems that can provide qualitative information to guide and reduce time-consuming experiments. Fabric shields are solved by a method called network modeling, which although there are no particles in this method, the algorithm of solution (i.e., the integration of the equations in time) ends up being relatively similar to the particle method. The main difference, but not the only one, is that the numerical algorithm doesn't have to search for contact and geometric entities between the nodes, because they aren't identified as particles. The numerical algorithm treats these nodes as pseudo-particles, which means that only the forces obtained by the finite element method are considered in these nodes, excluding the contact between eachother. Obviously, the contact forces developed between the rigid sphere and the nodes of the mesh are considered. This is a technique that allows the simulation of breakage, tearing and material separation (segregation) due to the impact at high rates of speed. Performing this modeling by using the classical finite element method procedures becomes extremely difficult, maybe impossible.
Nowadays exists a wide variety of materials that can be used to fabricate fabric shields such as Kevlar, zylon and other aramid-based materials. Zylon is a polymeric material produced by the Toyobo Corporation, see Taylor and Vinson (1990), and has a multiscale structure constructed from polybenzoxale (PBO) microscale fibrils, which are bundled to form yarn, which are then tightly woven into sheets. Experimental ballistic data on zylon has recently became available, see Zohdi and Powell (2006) for more details.
Experimental ballistic test for the evaluation of such materials are extremely expensive and time-consuming. For these reasons is necessary to develop numerical models to analyze fabric shields.
This research contributes to the development of a FE model with a new constitutive equation for fabric shields. It is the result of a combination of the theory of damage with the theory of elastoplasticity under large deformations. The hardening parameter is developed and subsequently modified to faithfully represent the physical and mechanical behavior of the fabric shields. This demonstration is concluded studying the physical behavior of this constitutive equation, as well as the calibration of the variables involved in this formulation. The advance in this type of approach, is to establish and respect maximum limits in the stresses values. Within this, there is a very positive advance in the numerical results obtained to represent the phenomenon of ballistic impact in fabric shields reliably to the reality. The rupture of the fabric shield and the damage of the yarn and coating must be similar when compared with an experimental result. In the numerical example a methodology is presented to design bulletproof vests. All fabric shield theory is minutely detailed and demonstrated in this article. An incremental layer analysis is performed through numerical simulations, seeking to determine the minimum number of sheets necessary to completely stop the projectile and, to investigate the rupture shapes of the bulletproof vest. Using the comercial software LS-DYNA this research sought to compare the results obtained by the program developed in this research, verifying the rupture configuration along the time. The final configuration of the ballistic impact experiment, presented by Zohdi (2014c) as a photography, are also used to check the final shape configuration.
Latin American Journal of Solids and Structures, 2020, 17(4), e272 3/42 The yarn is comprised of one-dimensional microscale fibrils, see Figure 1. The yarn has a coating to increase the resistance of the material. The fabric shields are made by yarn with coating arranged to form two-dimensional planar sheets. The fabric shield should be able to support a ballistic impact as shown in Figure 1.

DESCRIPTION OF THE COMPUTATIONAL NETWORK MODEL TO SIMULATE FABRIC SHIELDS
Basically, if the properties of the fibrils are known, the structural scale properties can be constructed. For structural fabric, a typical quantity of interest is the global tensile force-deflection response. The compressive response is usually of little interest.
The numerical model consists to consider a network of coated yarn that captures stretching of interconnected yearn networks with coating contributions and interaction with impact projectiles, incorporating contact with friction. The deformation of the fabric is calculated by solving the coupled system of differential equations for the motion of lumped masses.
A real fabric shield is presented in Figure 2.

DYNAMICS OF A NETWORK MODEL OF FABRIC AND PROJECTILES
The dynamics of the lumped masses are given by where = 1, 2, . . . , , is the number of lumped masses, is the mass of a single lumped mass (i.e., the total fabric mass divided by the total number of masses), is the position of node , is the total force of node , represents the normal contact forces contribution to node (contact between the fabric node and the projectile), represents the friction contact force contribution to node , and represent the axial contributions of the four yarns and coatings intersecting at mass , respectively, is the external damping, that can be used to simulate the effects of the environment and is the gravity force. The forces from the ℎ surrounding yarn-segment (there are four of them for the type of rectangular weaving pattern considered, see Figure 1) acting on the ℎ lumped mass is . Clearly, is a function of the mass positions ( ), which are all coupled together, leading to a system of equations. The dynamics of the projectiles are given by where = 1, 2, . . . , , is the number of particles, is the number of fabric nodes in contact with the projectile, is the mass of a single projectile, is the position of projectile , is the total force of the projectile. Note that the contact forces and presented in equation (2) is the same used in equation (1) for the correspondent lumped mass in contact with the projectile .

Mechanical forces of a yarn-segment formulation under finite elastic deformation
To develop the constitutive equation for yarn-segments the following conditions are assumed: (1) it has uniaxial-stress condition, i.e., the forces only act along the length of the yarn-segments; (2) the yarn-segments remain straight, undergoing a homogeneous (axial) stress state, (3) the compressive response of a yarn-segment is insignificant and (4) the buckling phenomena is ignored.
In terms of one-dimensional constitutive laws, the stored energy of a single fibril is given by where is the Young's modulus and is the Green-Lagrange strain. Using equation (125), the second Piola-Kirchhoff stress is given by It is convenient to work with quantities expressed in terms of the stretch ratio , that can be obtained as = , (5) where is the deformed length of the fibril and is its original length. The Green-Lagrange strain is obtained from (119) as following For this relaxed one-dimensional model, the Jacobian of is = and the stresses defined in equation (116) can be rewritten by When ≤ 1 the yarn has compression and when > 1 it has tensile stress. The compression is neglected, so, when ≤ 1 the first Piola-Kirchhoff stresses is enforced to be zero, i.e., = 0. The equations (118), (120) and (7) results in = = .
Latin American Journal of Solids and Structures, 2020, 17(4), e272 5/42 The first Piola-Kirchhoff stress can be expressed by where is the force of the yarn on the undeformed cross-sectional area , in the reference configuration. Using the relations presented in equations (5), (8), (9) and (10) the force carried in the yarn-segment results in Substituting (4) and (6) into (11) lies to So, the yarn force presented in equation (1) can be defined by + denotes the position vector of the endpoint connected to the lumped mass, − denotes the endpoint that is connected to its neighboring mass and ‖•‖ indicates the Euclidean norm in ℛ 3 . See the geometrical interpretation of the yarn force in Figure 3. Note that the deformed length of the fibril is calculated by

Mechanical forces due to the yarn and coating segments under finite elastic deformation
In this section the mechanical forces developed in the yarn-segment with coating is defined. Because of the networking model's kinematic, the following relations are valid Observing the yarn force presented in equation (13), the mechanical forces for the yarn and coating can be defined, respectively, by A simple estimate of the increase in stiffness due to the coating can be determined by summing the forces As explained in (17), for the coating contribution, the enforcement of zero compressive stress is not necessary. In practice, for ballistic applications, the enforcement of this condition makes little difference since the material is nearly always in tension. The cross-section area of a microfibril is a circle of radii . Let be the number of microfibril in a yarn-segment, the radii of initial cross section of the yarn-segment and the thickness of the coating. So, the cross-section area of the yarn and coating are, respectively, For Zylon each yarn approximately contains = 350 microfibrils.

Damage evaluation in a yarn-segment network
The damage phenomenon was initially incorporated in fabric shields by Zohdi (2002b). The idea is considering that the rupture of most fibers used as shielding is not sudden, but rather gradual. For Zohdi Zohdi (2002b), this is attributed to inhomogeneous rupture of microscale filaments which comprise a single fiber, as well as internal microscale friction, which impedes micro-filament pull-out between the individual micro-filaments. Basically, on the continuum scale, this progressive degradation in each fiber comprising the fabric net can be represented by a damager parameter such that it can construct a new Young's modulus in each time as following where and 0 is the Young modulus in the initial of the analysis, i.e., 0 = ( = 0). The damage variable for each yarn-segment typically has an evolution law associated with it, which represents progressive stretch-induced damage. Specifically, for a yarn that is undamaged, = 1, while for a yarn that is completely damaged, = 0. Hence, in the beginning of the analysis ( = 0) = 1. Note that when the parameter goes to zero the physical material parameter also goes to zero. The damage parameter proposed by Zohdi (2002b) is governed by the following over-stress evolution equation is the Cauchy stress, is the Cauchy stress limit to starts the damage of the fibers and is a rate parameter. It's important to mention that this formulation considers that the Cauchy stress is positive in tensile case and zero in compressive case, i.e., if ≤ 1 then = 0.
For Zylon, according to the manufacturer, for a fibril, = 1.03 and the Young's modulus is = 180 . From equation (121)  . The rupture of the material Zylon is presented by Zohdi and Powell (2006) and is supposed to occur when the stretch is around ≅ 1.06, see Figure 4. The damage parameter was also proposed by Zohdi (2014c) as shown in the following equation where the parameter assumes the following values The equation (29) can also be rewritten as where ∆ is the increment of time. The equations (27) and (31) are valid when the stretch is equal or higher than . The inconvenience of using this kind of progressive damage is the control of the maximum stresses and stretches. It means that the rupture configuration can be significantly different from the final experimental shape.
Another equation to simulate damage was proposed by Zohdi and Powell (2006) in the following manner where ( ) is the stretch of the yarn-segment at time , , is the critical stretch of the yarn-segment and ≥ 0 is a damage decay parameter which is related to the amount of fibril misalignment that exists in the yarn. The authors in Zohdi and Powell (2006) suggests that a best fit value is ≈ 124. The equation (32) has a numerical instability problem to represent the case of no damage generated, → 0, and sudden rupture, → ∞.
Basically, a parameter ( ) is created to describe the failure of the yarn-segment by checking whether a critical stretch has been attained or exceeded, i.e., ( ) ≥ and additionally, this parameter will track the progressive damage with a simple damage (with isotropic material).
This parameter, also called damage representation, was also proposed by Zohdi (2014c) in the following manner Analyzing equation (33) is concluded that when → ∞ the type of failure tends towards sudden rupture and when → 0 there is no damage generated. This formulation is better than the one presented in (32) because it doesn't have any division by zero and all values of damage are positive, i.e., attend equation (26). This work suggests that a best fit value for equation (33) is ≈ 135. This value gives a curve similar to the theoretical equation presented in (32), that is approximately equal to the experimental curve obtained by Zohdi and Powell (2006). The curves behavior of equations (32) and (33) are presented in Figure 4. The parameters used in equation (32) and (33) are = 124 and = 135, respectively. It is important to mention that in practical procedures, is suggested to consider the total damage when approaches 0.02, see Figure 4. The equations (27), (32) and (33) indicate that damage is irreversible, i.e., is a monotonically decreasing function.
The same formulations are applied for the yarn and coating.

Material behavior
It is important to mention that the yarn of fabric shield doesn't have the standard plasticity with hardness, like the Von Mises criteria for example. The plasticity can be considered due to the coating of the yarn-segments.
To incorporate the elastoplastic formulation in the yarn-segment with coating, a plastic stretch strain model is used. A complete review of elastoplasticity in large deformation is presented in the end of this article. Consider one dimension and the relation presented in equation (9), is possible to write the following decomposition where and is sometimes referred as an elastic and plastic stretch strain, respectively. The Green-Lagrange strain is defined by where and are the elastic and plastic part of Green-Lagrange strain, respectively. The elastic part of the second Piola-Kirchhoff stress presented in equation (4) can be rewritten as = . (36) Introducing equation (121) into (37), result in Latin American Journal of Solids and Structures, 2020, 17(4), e272 9/42 The above equation can be rewritten as From equations (8) and (9) is possible to determine the Cauchy stress tensor The equation (40) can also be written as where = ( ) −1 . The variation of the plastic stretch strain can be computed by the following relation is a rate parameter, is the yield stress (connected to stretch sensitivity), is the initial yield stress and ℋ is a hardening variable.
It is important to mention that when | | ≤ the plastic strain is zero, i.e., = 0. In this case, the plastic stretch tensor is = 1 and the second Piola-Kirchhoff stress tensor is = , where = . These relations can be verified in equation (39). The update of the plastic part of the stretch tensor can be computed from (42) as following where and ∆ is the size of the time increment. Analyzing the equations (40), (44) and (45) is possible to conclude that to compute the variation of the stretch tensor ∆ is necessary to know the Cauchy stress tensor , that is unknown in the current time step and need to be calculated by equation (40). So, to solve it, the value of to be used in equation (40) is chosen to be the same of the previous step and its update is done in the end of the time step by using the equations (44) and (45) with the value of the Cauchy stress tensor obtained in equation (40). In the beginning of the analysis or in the elastic range the value of the plastic stretch tensor is = 1.

Integration scheme to obtain the stresses in plasticity
The present section introduces one of the contributions to the field. Basically, the author uses an adapted algorithm from elastoplasticity and includes damage to it, to make the control of the maximum stretch and stress. Moreover, the parameters are calibrated to an abrupt rupture, not allowing for hardening.
The integration scheme of the elastoplastic constitutive equation is shown in Algorithm 1. The formulation is presented in Appendix B and Appendix C.

3)
Compute the elastic part of the stretch:

4)
Compute the elastic part of the Green-Lagrange strain:

8)
Compute the variation of the plastic part of the stretch tensor ∆ from equation (45).

9)
Update the plastic stretch tensor:

10)
Compute the elastic part of the stretch:

11)
Compute the elastic part of the Green-Lagrange strain:
To illustrate the behavior between stretch and stress by using the elastoplastic algorithm, is considered a simple numerical example of a Zylon yarn segment, with length = 0.00254 , Young modulus = 180 and critical stretch = 1.03. As explained before, the Cauchy stress limit lies to = 5645.43 and the rupture occurs when ≅ 1.06. In order to generate a tensile force, in each time increment of ∆ = 10 −6 the yarn segment was increased its length in ∆ = 2.03 • 10 −7 . The case that have no damage in the yarn segment ( = 0) and have only elastoplasticity with parameter = 10 4 , the numerical result obtained is shown in Figure 5. Observe that this constitutive law doesn't give values of stresses higher than the limit established, i.e., ≤ .

Mechanical forces due to the yarn and coating segments considering plasticity and damage
In this section the real force developed by the yarn and coating segments considering plasticity and damage effects is presented. Substituting equation (39) into (11) and using the relations presented in equations (16) where , ( ) and ( ) are defined in (14), (22) and (23), respectively. The Cauchy stress tensor for the yarn and coating are respectively Note that and are calculated independently by the Algorithm 1. To illustrate the behavior of the curve correlating stretch and stress taking into account the elastoplasticity and damage, the same example of the yarn segment presented in section 4.4.1 is considered. The equation (33) is used to calculate the damage. In case that have a rather gradual damage in the yarn segment ( = 135) and elastoplasticity with parameter = 10 4 , the numerical result obtained is shown in Figure 6. For the plastic part of the graphic, the damage curve is similar to the one presented in Figure 4. It's a gradual and exponential loss of stiffness. It is important to mention that the Figure 6 can be also obtained with the parameters = 135 and = 0. This formulation presented in section 4.4 is more stable than the others ones presented. To simulate a sudden rupture, the parameter → ∞. In this numerical example it should be = 10 4 , see Figure 7. The material stress-stretch relations presented in Figure 6 and Figure 7, obtained by the constitutive formulation developed in this research, are in agreement with the ones presented in Zohdi and Powell (2006). This formulation presented in this section will be used in the numerical examples brought at the end of this article to simulate the yarn-segment and coating. It is important to mention that each material has its damage parameter as well as its physical parameters. This constitutive equation means a considerable contribution to the study of fabric shields. It is important to note that the hardening function has been adapted and that the plasticity algorithm has the main function of controlling the tensions acting on the material, preventing them from being higher than the yield stresses. The abrupt rupture constitutive model perfectly simulates the zylon behavior as its curve shows no yielding.

MECHANICAL FORCES DEVELOPED BY THE CONTACT BETWEEN YARN-SEGMENTS AND PROJECTILES
The fabric/rigid-body interaction was developed by Zohdi (2010a) and Zohdi (2014c) based on classical contact algorithms. A variety of contact algorithms is presented by Wriggers (2002).
The contact between a rigid projectile with a fabric shield is performed as a contact point-to-point, i.e., between the fabric node and the projectile, that is described by the position of its center of mass. Basically, according to Zohdi (2010a) and Zohdi (2014c), if penetration of a fabric node into the projectile is detected ("contact"), the following procedures are enacted: (1) move the fabric node to the closest point on the rigid body's surface and (2) if there is friction, then a stick condition is assumed, i.e, assuming a common position and velocity for the point on the surface of the projectile and the fabric node. The friction force is then checked against the static limit, which if violated, enacts a sliding friction force. See Figure 8 for a geometrical representation.

Geometrical interpretation for a fabric node
Assuming interpenetration for a node , see Figure 8, the contact force on a node can be calculate as where is the lamped mass of the node of the fabric. From above equation, the corrected contact position is given by where the predicted contact position is defined by where , ( + ∆ ) is the velocity of the node in time step + ∆ , ( ) is the velocity of the converged value from the previous time step ( ) and is the trapezoidal time-stepping parameter. The velocity of the predicted contact position is calculated by where , ℎ ( ) are all forces other than the contact forces. In this work, the forces other than contact are considered The equation (52) can be rewritten as Substituting (53) into (55), lies to To updated each fabric node , basically, the projectile is frozen and a new position of the fabric node is computed taking into account the contact conditions, i.e., if occurs penetration into the projectile, the position of the fabric node must be correct, see Figure 9.

Procedure to calculate the Fabric motion
The procedure used in this research is described as following. In the beginning of each time increment, the first step consists to calculate the predicted contact position , ( + ∆ ) by using equation (52). If the contact occurs, see Figure 9, the fabric node must move to the closest point on the rigid body's surface. For this purpose, first is assumed that the projectile-fabric surface "sticks". So, the fabric node is corrected to where , ( + ∆ ) is the closest point on the surface of the contacting object. For an in-depth review of these types of methods see Wriggers (2002). Additionally, values to the velocity must be assigned. Considering "stick" case, is assigned the velocity of the fabric node , ( + ∆ ) to that of the closest point on the body , ( + ∆ ), for node , as Latin American Journal of Solids and Structures, 2020, 17(4), e272 14/42 The friction force is calculated using where the contact force is computed by equation (50) and ̇ can be considered to be the acceleration of the projectile, since it is considered "stick" case. The acceleration ̇ can be computed by equation (83).
If the friction force exceeds the theoretical limit then the assumption that there is "stick" is incorrect, and sliding must occur, resulting in a nonmatching velocity of with the (frictional) projection where the velocity of the node of the fabric at the closest contact point , is assumed to be

Procedure to calculate the projectile motion
After the calculations culminating in equation (62) have been completed, the velocity of the projectile is recalculated using (summing the nodal contributions) where is the number of nodes in contact and In this work, these forces other than contact are considered to be The equation (65) can be rewritten by Substituting (64) into (67) lies to The velocity of any point on the surface of the body can be determined by simply calculating The projectile's angular velocity and rotation are determined in a similar manner by integrating the equations for an angular momentum balance where ̅ , is the mass moment of the projectile, is the angular velocity and is the sum of all moment contributions external to the projectile, around its center of mass. A general rigid body solution algorithm is provided by Zohdi (2010a).
Solving equation (70): Then, based in equation (71), the time discretization is done by the following equation: where ∅ is a scalar value with range 0 < ≤ 1. It is important to mention that when the parameter = 0 the system is solved by explicit method and for 0 < ≤ 1 the system is solved by implicit method. The equation (72) can be rewritten as The moment can also be calculated by where is the friction force between the nodes of the fabric and the projectile. This force is applied at the contact point on the surface of the projectile.
In this work, the rolling of the projectile is taken into account. Therefore, the friction force is modeled by assuming that sliding and rolling may occur between the contact pair, i.e., whether it is pure sliding, sliding with rolling or pure rolling.
The rotation vector of a particle is obtained by The time discretization is done by the following equation: Finally, for a point ′ belonging to a projectile in the initial configuration placed → ′ � =0 from its center, its current position is defined by: Based on equation (64), the dynamics of the projectile is given by It is important to mention that in this work spheres are used to model projectiles. Therefore, the contact point can easily be found and is not necessary to find the current position of any point at the projectile's surface related to its center. Then, it is neglected in the numerical algorithm.

Damping force
The other force that can be taken into account is the damping effect given by where ̇ is the external damping, that can be used to simulate the effects of the environment.

Gravity force
The other force that needs to be taken into account is the gravitational force, which acts in a downwards (− direction) and is given by where is the gravity acceleration and � is the unit vector in the direction.

NUMERICAL SOLUTION SCHEME
In order to describe the overall time-stepping scheme, is considered the dynamics of a single ℎ lumped mass as presented in equation (1) and the dynamics of a single ℎ particle presented in equation (2). The rotation in the node of the fabric is neglected and is only considered in the projectiles. Basically, comparing the fabric shields formulation with the standard contact mechanics formulation used in finite element method, the fabric nodes can be considered to be the "slave node" and the projectiles the "master node (target)". Within this meaning, the algorithm has the basic aim to calculate the converged positions and velocities for all fabric nodes, taking into account that in the end of each iteration the positions and velocities of the projectiles must be corrected. So, in the end of each time step, the values converged to the correct solution.
For this purpose, let the position vector be defined by and the velocity vector defined by where = 1, 2, ⋯ , , is the number of lumped mass, = 1, 2, ⋯ , , is the number of projectiles.

Adaptive time-stepping scheme
The time-stepping scheme used in this work is based on the original works presented by Zohdi and Powell (2006), Zohdi and Wriggers (2008), Zohdi (2010a) and Zohdi (2010b). The only difference is that, in this article, the convergences is verified not only considering displacements, but also velocities. Employing a trapezoidal rule, such that 0 ≤ ≤ 1, the following relations are valid: The position can be computed by applying the mid-point rule, which results in Substituting equation (84) into (85) is obtained The equation (86) can be solved recursively by recasting the relation as and the same for equation (84) as where superscript and + 1 indicates the indices of the current and the future time-step in the discretization, respectively. Basically, the superscript is a time interval counter. This is a generalized one-step scheme, with ∈ (0,1], in which, when = 0 is the explicit Euler integration scheme, while = 1 is the fully implicit Euler integration scheme. The implicit implementation of the position and velocity updates, as is evident from equations (84) and (86), requires an iterative solution for +1 and +1 . The set of equations represented by (87) and (88) can be solved recursively by recasting the relation as where = 1,2,3 … is the index of iteration within time-step + 1. These equations can also be represented in the following form: where, The terms ℛ , and ℛ , are remainder terms that do not depend on the solution. The convergence of such a scheme is dependent on the behavior of and . The form of the error function is presented below as , +1, ≝ +1, − +1 , , +1, ≝ +1, − +1 .
Based in the formulations presented by Zohdi (2014a), a sufficient condition for convergence is the existence of a contraction mapping regarding the displacements. In this work, the author utilizes a contraction mapping regarding both the displacements and velocities, as where 0 ≤ +1, < 1 and 0 ≤ +1, < 1 for each iteration . Then , +1, → and , +1, → for any arbitrary starting value +1, =0 and +1, =0 , as → ∞, which is a contraction condition that is sufficient, but not necessary, for convergence. The convergence of equation (89) is scaled by ∝ ( ∆ ) 2 (see equation (95)) and the convergence of equation (90) is scaled by ∝ ∆ (see equation (97)). Therefore, the contraction constant of and is directly dependent on the magnitude of the interaction forces (‖ ‖), inversely proportional to the masses and directly proportional to ∆ . Thus, decreasing the time-step size improves the convergence. Zohdi (2002a) proposed an approach to maximize the time-step sizes to decrease overall computing time and keep the numerical solution's accuracy. This propose takes into account only the position vector's errors. Here, this propose is extending also for the velocity vector's errors. Basically, these following expressions are assumed: where and are constants. The error resulted within an iteration is according to where , +1,0 = +1, =1 − and , +1,0 = +1, =1 − . These above equations can be rewritten as According to Zohdi (2014a), in the end of the time-step iteration, , +1, = and , +1, = , where and are tolerances for the position and velocity vectors, respectively. So, assuming that and are constants, a new smaller step size can be defined as where is the number of desired iterations. According to Zohdi (2014a), the expression in (109) can also be used for time-step enlargement, if convergence is met in less than iterations, usually chosen to be between 5 to 10 iterations. The numerical solution scheme is shown on Algorithm 2 in Appendix A. This numerical algorithm has a new approach, being presented in a much more complete and detailed way. However, it is adapted fromZohdi (2014b), Zohdi (2014a), Zohdi (2014c) and Zohdi (2010a).

NUMERICAL EXAMPLES
In this section, numerical results using the theory described in this article are presented. The examples are modeled in the currently developed program and in the comercial program LS-DYNA, a module of ANSYS. The examples where simulated in a machine with Intel Core i7 7700HQ (2.80GHz processing), graphic card 4095MB NVIDIA GeForce GTX 1050 Ti, 256 MB storage SSD and 1 TB SSHD-8GB.

Fabric shield composed by one sheet with yarn and coating
The basic aim of this example is to apply the formulations presented in this article to simulate, numerically, structural fabric shields and investigate in details these systems and the behavior of the yarn-segment with coating.
This model consists of one sheet of fabric shield and a spherical projectile with an incoming velocity of = 500 / , see Figure 10. In this model, the boundaries of the fabric shield are completely restricted. It is worth mentioning that the geometry and some material parameters are based on the work developed by Zohdi (2014c). This is done with the objective to compare the numerical results and the final deformed configuration shape. Although damping and gravity effects are considered in this analysis, it is important to remark that their contribution are negligible in the present case study, because of the high projectile speed. The geometry of the system is presented in Figure 10 and the parameters used in the numerical simulation are given as following.
It is observed in Figure 11a that the projectile is in eminence for impact on the fabric shield. In Figure 11b, the projectile penetrates the fabric shield and in Figure 11c, it overtakes it completely. In Figure 11d, the projectile is viewed well away from the fabric, leaving a hole.
The damages of the yarns and coatings at time 3.2 • 10 −4 are presented in Figure 12a and Figure 12b, respectively. It can be seen that the yarns that were not destroyed during the impact, were not damaged, i.e., the damage parameter continues equal to 1. The coating segments have already suffered damage, especially in the region around the hole. The farthest coating segments from the hole were not damaged. The stresses of the yarns and coatings at time 3.2 • 10 −4 are presented in Figure 13. To represent geometrically the physical phenomenon, generated by the ballistic impact between the projectile and the fabric shield, some configurations were selected over time, as shown in the Figure 14. The velocity of the projectile after crossing the fabric shield is = 367.13 / . As shown in the figures, it is observed that the projectile pierced the fabric shield. Analyzing the value of the velocity of the projectile, it is observed that only one layer of shield reduced the velocity by approximately 26.6%. So, as a way to stop the projectile, and consequently prevent full penetration into the fabric shield, it is necessary to add more sheets (layers) of fabric shield. This discussion will be presented in the next section.

Modeling using LS-DYNA
The same example was modelled using LS-DYNA. The projectile was defined as a solid element (as particle) and the structural fabric as a shell element, with a standard Belytschko-Tsay membrane formulation suitable for fabrics and materials where the flexural stiffness is negligible. To reduce the storage cost and the analisys time requirement, shell elements where used instead of bars. For the structural fabric, its material model was chosen from the options with viscoplastic deformation rate ("power law plasticity"). In order to describe the interaction between the bodies during the impact, the automatic type of surface-to-surface contact was used.
The numerical results of this modeling are illustrated in Figure 15 and Figure 16. The rupture configuration shown in Figure 15 has the same shape of the one obtained by the program developed in this research, see Figure 11, Figure 12 and Figure 13.
The lateral view presented in Figure 16, generated by the LS-DYNA analysis, has the same behavior as the one obtained by the numerical program developed in this research, see Figure 14.
It can be observed that the rupture format is very similar between the two kind of analysis, one modelled with shells (LS-DYNA) and the other one by bars (this research). The small difference can be associated with the different way to model the fabric structure. The author affirm that the rupture configuration obtained in this research has more accuracy with the real phenomenon, observed experimentally in Zohdi and Powell (2006).

Fabric shield composed by several sheets with yarn and coating
Since the objective of the analysis is to dimension a bulletproof vest capable of stopping the projectile, one begins to investigate how many layers of shield are necessary to stop the projectile. For this purpose, another layer of fabric shield is added 0.0004 meters away from the previous one. It should be noted that in each layer is considered the same physical and geometric properties of the fabric shield presented in the previous section. In this way, this procedure is done several times until the goal of stopping the projectile is reached. So, twelve kind of reinforcement situations were modeled for the fabric shield, where in each of them another layer of fabric was added keeping the same distance from the previous one.
For the sake of simplicity, in this example, only the numerical results for 12-layer fabric shield will be presented. In that model, the discretization contains 120.000 nodes and 237.600 elements.

Figure 17
Numerical results for the model of twelve layers at time 7.7 • 10 −5 .
The Figure 17 shows the projectile in a configuration where its velocity is almost zero. The Figure 17a illustrates the Cauchy stress on the yarn and Figure 17b illustrates the Cauchy stress on the coating. Figure 17c shows that the yarns were undamaged and Figure 17d illustrates how the coating has been damaged.
To represent geometrically the physical phenomenon, generated by the ballistic impact between the projectile and the fabric shield with twelve sheets, some configurations were selected over time, as shown in the Figure 18. It is concluded in this model that nine layers of fabric shield are required to completely stop the projectile. With twelve sheets, the projectile stops completely and the reaction of the fabric shield change the direction of its velocity. In this way, the structural fabric pushes the projectile to the opposite direction from the initial one. After the velocity reaches the null value, it is inverted and its value begins to increase gradually, as expected for the physical behavior. To better understand the analyzes performed, a comparative study of the velocity of the projectile over time was done for each model studied. Figure 19 describes 4,00E-06  9,00E-06  1,40E-05  1,90E-05  2,40E-05  2,90E-05  3,40E-05  3,90E-05  4,40E-05  4,90E-05  5,40E-05  5,90E-05  6,40E-05  6,90E-05  7,40E-05  7,90E-05  8,40E-05  8,90E-05  9,40E-05  9,90E-05  1,04E-04  1,09E-04  1,14E-04  1,19E-04  1,24E-04  1,29E-04  1,34E-04  1,39E-04  1,44E-04  1,49E-04  1,54E-04  1,59E-04  1,64E-04  1,69E-04  1,74E-04  1,79E- The Figure 20 shows the configuration of the system, composed by the projectile and six layers of fabric shield, at time = 1.15 • 10 −4 . Figure 20a shows that the projectile has penetrated more than half in all fabric layers, and for this reason, its velocity will remain constant and equal to = 131,50 / , since there is no more fabric capable of resisting its movement. It is observed in Figure 20b that the fabric around the projectile is destroyed. The Figure 21 illustrates the exact configuration in which the projectile is completely stopped by the fabric shield composed by nine layers. Analyzing Figure 21a, it is perceived that the projectile penetrated less than half of its diameter. For this reason, the last layers of fabric still offer resistance to the projectile penetration. As explained earlier, once the projectile penetrates beyond its half in the last layer of the fabric shield, practically the projectile will continue without any resistance and with constant velocity, see the time-history of projectile velocities for each model in Figure 19. Another study was carried out observing the projectile in the configuration defined at time 7.7 • 10 −5 . At this specific time, the velocity of the projectile in each model was measured. In this way, a curve was constructed that related the number of layers of fabric shield with the respective speed of the projectile, as shown in Figure 22. It is observed in this curve that the relationship between the velocity ( ) and the number of layers of fabric shield ( ) is non-linear. A study was carried out to define the trend line of the curve, thus obtaining a fifth-degree polynomial, defined by the equation = −0,0123 5 + 0,4375 4 − 5,7882 3 + 37,509 2 − 156,82 + 493,45. With this, it has been proven that there is a highly non-linear relationship when layers are added to the fabric shield. Another pertinent investigation is to observe the shape that the fabric shield ruptures. The hole in the fabric shield (composed by one layer) caused by the ballistic impact obtained through numerical simulation is represented in the Figure 23a. The experimental trial was conducted by Zohdi and Powell (2006) and the photograph of the fabric shield, after their experimental test, is shown in the Figure 23b, being used by the authors as a valid comparison for their work. Comparing the two figures, one obtained by numerical simulation and other by experimental results, it is concluded that both rupture configurations of the of fabric shields are practically the same.
The Figure 23c and Figure 23d illustrate the configuration of the projectile when it is almost stopped.

Figure 23
Results obtained numerically and experimentally.
The shape of the hole left in the structural fabric after the ballistic impact of the projectile was the same in all models. Of course, the velocity of the projectile had a distinct behavior, as shown in Figure 19.

Modeling using LS-DYNA
The same previous example was modelled using LS-DYNA. It´s important to remark that the objective in this section is to investigate the rupture shape of the fabric shield and the number of sheets necessary to completely stop the projectile. The numerical results of this modeling are illustrated in Figure 24.
The analisys from LS-DYNA gives a similar result to the one presented by the program developed in this research. The only diffence is the number of sheets necessary to stop the projectile, in which LS-DYNA required a higher number of sheets, it means that, in this research, the number of sheets necessary to completely stop the projectile is 9, while in LS-Dyna it is 13. Despite this, the same behaviour occurs in both analisys, i.e., the projectile pierces the fabric structure, stops and then, bounces backward.

CONCLUSIONS
The basic aim of this research was to present a methodology to dimension bulletproof vests formed by fabric shields, capable of resisting ballistic impacts. For this purpose, a complete revision was made in the theoretical formulation applied to the fabric shields to obtain a result very close to reality, especially when compared with the LS-DYNA results and the experimental ballistic test performed by Zohdi (2014c). Particular attention was given to the constitutive equation of the material, the yarn and the coating. The reason for this was to respect the maximum stress that the material can withstand, especially at the moment of impact. To this end, the elastoplastic theory was adapted to the constitutive equation of the material. Then, the damage parameter was added in the material constitutive equation. In this sense, this research contributed a lot to the scientific community, in order to more accurately determine the number of sheets in the bulletproof vest to stop a projectile. If tensile stresses were not respected, a smaller layer of fabric shield would probably be required. Thus, the author conclude that the numerical results obtained, after comparison with LS-DYNA, were very satisfactory.
It was observed that, as one layer of fabric is very close to the other, a representative number of layers is necessary to prevent the bullet from overtaking the vest. It suggests for future work the study of a configuration for these structures so that the layers can be further apart from each other, increasing the efficiency of the system and consequently, the resistance of the shield.
Algorithm 2 Numerical solution scheme.

APPENDIX B: OVERHAUL OF SOLID MECHANICS
From the continuum mechanics, the second Piola-Kirchhoff stress tensor is defined by where is the first Piola-Kirchhoff stress tensor that correlates nominal stress per unit reference area and is the deformation gradient defined as In above equation, and are the position vector of the nodes in current configuration and reference configuration, respectively. The Cauchy stress tensor associates the superficial force vector (or stress vector per unit of area in the current configuration) with the unit normal vector at the surface in current configuration defined by = . (112) Using the Nanson theory and considering = (113) and = , is possible to demonstrate the following relation The Kirchhoff-Trefftz stress tensor is defined by Note that the second Piola-Kirchhoff stress tensor is symmetric and doesn't have a physical interpretation, the first Piola-Kirchhoff stress is nonsymmetrical, the Kirchhoff-Trefftz stress tensor is symmetrical, the Cauchy stress tensor is symmetrical and they are energetically conjugate by this relation :̇= :̇. The displacements gradient is defined by where is the displacements vector of the nodes and is the identity matrix. The quadratic stretches tensor or right Cauchy-Green strain is defined by = . (118) The Green-Lagrange strain tensor is calculated by The stretch tensor is calculated by The Green-Lagrange strain tensor can be written also as The following strains belong to the Hill family: where is a integer value and is the stretch. The equation (121) corresponds to (122) when the parameter is equal to 2. The rotation tensor can be calculated by the following polar decompositions where is the left stretch tensor. The rotation tensor is orthogonal, i.e., = +1. For hyperelastic materials, the following relations that correlates stresses with strains are valid: where ( ) and ( ) are the function of specific strain energy. A hyperelastic material is necessarily elastic and reversible. The tensor that represents the tangential stiffness is given by = 2 2 = . (127)

Introduction
In order to incorporate the plasticity in the yarn segments, an overhaul of elastoplastic constitutive equation for large deformation is needed.
Elastic-plastic behavior is characterized by a response from the material, initially elastic, and from a certain level of tension by essentially plastic behavior. The plastic behavior of the material is usually accompanied by an invariance of its volume. The typical behavior of an elastoplastic material submitted to uniaxial tension is shown in Figure 25 When the initial limit tension 0 is reached, the limit tension value may or may not remain constant with increasing deformation. If this value does not depend on the increase of the plastic extension, the material has a perfectly plastic behavior, see Figure 25a. If, on the other hand, the value of the tension increase with the growth of the plastic extension, it is said that the material is undergoing a hardening, see Figure 25b.
Over time, several studies on elastoplasticity have been carried out. One can cite, for example, the works done by Chen and Han (1988), Hill (1950), Kachanov (2004), Lubliner (1984), Pimenta (1992), Simo and Ortiz (1985), Simo and Taylor (1985), Simo (1992) and Wriggers, Eberlein and Reese (1996). In this section, a succinct approach on the formulation of elastoplasticity under large deformation will be presented. Through this formulation, the elastoplasticity algorithm used in this research will be developed. It should be noted, that after obtained the plasticity algorithm, some modifications will be done, i.e., the function corresponding to the hardening parameter will be adapted and the damage parameter will be included, as presented in chapter 4.
Kinematics Consider the composition where is total deformation gradient, is the elastic part of the deformation gradient from the natural to the current configuration and = corresponds to the plastic part from the initial reference configuration to the natural one, see Figure 26. In most cases, one sets (0) = . The elastic and the plastic parts have the following polar decompositions where is the right elastic stretching tensor, is the left elastic stretching tensor, is the elastic rotation tensor, is the right plastic stretching tensor, is the left plastic stretching tensor and is the plastic rotation tensor. As the natural configuration is not fixed, this means that ̇≠ whenever the material yields. So, at each instant, is considered a new fixed reference configuration that is set instantaneously coincident with the current natural one. The deformation gradient from this new reference configuration is now denoted by and given by see again Figure 26 for a geometrical representation. This gradient has the following polar decomposition At each instant, the following identities hold = , = , = .
Considering these relations = =̇ −1 and = =̇ −1 , where =̇ and =̇, the velocity gradient from the instantaneously fixed natural configuration is Changing and ̇ to the same reference configuration, one concludes that Hence, one has Since the equation (137) is valid, one has For the same reason one may instantaneously write ̇=̇.
The elastic Green-Lagrange strain tensor, as shown in equation (121), is defined by Considering the principal axes of the right Cauchy-Green strain tensors , once has the spectral decomposition = ⊗ , where = � , and using the relations present in (121) one may write as the spectral decomposition below where The total strain from the fixed natural configuration is defined by In view of (133), at current time, the following identity holds = .
Isotropic plasticity In plasticity the elastic domain encompasses all admissible stress states and is described by the yield function where : × ℝ ℎ → ℝ and is the ( ℎ × 1) vector of the hardening parameters. does not depend on hardening parameters in the particular case of ideal plasticity. The surface in described by ( , ) = 0 is called yield surface. In isotropic plasticity, first is assumed that elasticity is isotropic and described by ̇=̇, and the yield function (147) is written as displayed below in which � is an isotropic function of its argument and equal to the uniaxial stress in a simple tension test. This means the following where is the unitary vector in the axial direction of the test and is the corresponding uniaxial stress. � is then called equivalent stress. On the other hand is the yield stress measured in a simple tension test and is the isotropic hardening parameter. The plastic dissipation per unit volume in the natural configuration is given by the internal dissipation, as displayed below where ( ) is the strain energy function per unit volume in the natural configuration of an elastic material and the operator is defined as : = ( ), where (•) indicates the trace operator. In order to formulate the constitutive rate equations of plasticity, is stated that the Principle of Maximum Plastic Dissipation is a necessary and sufficient condition for the determination of . Mathematical Programming asserts that (153) can be transformed to the following unconstrained problem where is the Lagrange multiplier. Thus, one has is the normal to the yield surface = 0 in the stress space. Since is an isotropic function of and the elasticity is isotropic, , and are coaxial. Thus, (155) implies that ̇−̇, , and are coaxial. It is easy to demonstrate the following result ̇=̇+̇, which is similar to the additive decomposition of the total strain in small strain plasticity. From and one arrives at and ̇= �̇−̇� .
The above equation leads to the following result Note that for ̇> 0 and introducing the equations (162) and (167) into (161) can be rewritten by ̇= �̇−̇�.
From (175) Introducing ̇= − �̇−̇� in the equation above, one arrives at (155). In metal plasticity, is assumed that the plastic part of the deformation gradient is isochoric. Thus, one has which implies and Σ = Σ .
According to (155), if is deviatoric, then ̇ is deviatoric too. This implies (179). In the case of isochoric plasticity, one has = and = . Stress integration In plasticity, the constitutive equations are rate equations. So, it is necessary to formulate a stress integration algorithm in order to obtain an approximate solution for the quasi-static problem at discrete instants { 0 , 1 , 2 , ⋯ , , +1 , ⋯ }.
Consider the material state at instant , at which , and are supposed to be known. Hence, from (128) one has If the polar decomposition = is performed, the stress tensor at this instant can be computed through The above equation lies to ̇= ̇ .
At instant +1 only the deformation gradient +1 is assumed to be known. The following trial elastic state is defined, assuming that the plastic part of the deformation has not changed, The polar decomposition is then performed and the strain tensor is computed with the aid of The stress tensor in this state is given by If = ( ) < 0, then the trial state is admissible, ∆ = 0, and If ≥ 0, there is an incremental plastic stretch ∆ such that The plastic stretch tensor at the end of the increment, denoted by +1 , is obtained from +1 through The incremental stress integration algorithm is formulated by the following constrained minimization problem, which must be solved for +1 , min * (Δ ) ( +1 , +1 ) = 0, where Δ = +1 − .
Considering that equation (207) can be rewritten as The value of ∆ is obtained from the above equation, just solving the following relation The so-called tensor of algorithmic tangent moduli is defined by Remember that (211) is referenced to the natural configuration at the instant and is valid for the logarithmic stress and strain tensors from this configuration. Plasticity algorithmic The geometrical interpretation of the update procedure is shown in Figure 27.