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

Abstract

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.

Keywords:
Fabric shield; bulletproof vest; ballistic impact

Graphical Abstract

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.

Technological advances in the area of fabric shields are presented by Zohdi and Powell (2006), Kollegal and Sridharan (2000), Lanir (1983), Lim, Shim and Ng (2003), Powell and Zohdi (2009a), Powell and Zohdi (2009b), Tabiei (1999), Zohdi (2002b), Zohdi (2007), Zohdi (2009), Zohdi and Steigmann (2002), Zohdi and Wriggers (2001b), Zohdi and Wriggers (2008) and Atai and Steigmann (1997). Among these, it is important to highlight the work of Zohdi and Powell (2006) that simulate a numerical experimental test of a fabric shield subjected to ballistic impact.

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.

DESCRIPTION OF THE COMPUTATIONAL NETWORK MODEL TO SIMULATE FABRIC SHIELDS

Figure 1
Representation of a network model of fabric (based onZohdi (2014c)).

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.

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.

Figure 2
Photography of fabric shield.

DYNAMICS OF A NETWORK MODEL OF FABRIC AND PROJECTILES

The dynamics of the lumped masses are given by

$m_i{\ddot{\mathit{r}}}_i=\underset{total}{\underbrace{{\mathbf{F}}_i^{tot}}}=\underset{contactforces}{\underbrace{{\mathbf{F}}_i^{con}}}+\underset{frictionforces}{\underbrace{{\mathbf{F}}_i^{fric}}}+\underset{surroundingyarn}{\underbrace{\sum _{I=1}^4{\mathbf{F}}_I^{yarn}}}+\underset{surroundingcoating}{\underbrace{\sum _{I=1}^4{\mathbf{F}}_I^{coating}}}+\underset{damping}{\underbrace{{\mathbf{F}}_i^{damp}}}+\underset{gravity}{\underbrace{{\mathbf{F}}_i^{grav}}},$ (1)

where $i=\mathrm{1,}\mathrm{2,}\mathrm{...,}{N}_f$, $N_f$ is the number of lumped masses, $m_i$ is the mass of a single lumped mass (i.e., the total fabric mass divided by the total number of masses), $r_i$ is the position of node $i$, $F_i^{tot}$ is the total force of node $i$, $F_i^{con}$ represents the normal contact forces contribution to node $i$ (contact between the fabric node $i$ and the projectile), $F_i^{fric}$ represents the friction contact force contribution to node $i$, $F_I^{yarn}$ and $F_I^{coating}$ represent the axial contributions of the four yarns and coatings intersecting at mass $i$, respectively, $F_i^{damp}$ is the external damping, that can be used to simulate the effects of the environment and $F_i^{grav}$ is the gravity force. The forces from the $I^{th}$ surrounding yarn-segment (there are four of them for the type of rectangular weaving pattern considered, see Figure 1) acting on the $i^{th}$ lumped mass is $F_I^{yarn}$. Clearly, $F_I^{yarn}$ is a function of the mass positions ($r_i$), which are all coupled together, leading to a system of equations.

The dynamics of the projectiles are given by

$m_p{\ddot{\mathit{r}}}_p=\underset{total}{\underbrace{{\mathit{F}}_p^{tot}}}=-\underset{contactforces}{\underbrace{\sum _{i=1}^{N_c}{\mathit{F}}_i^{con}}}-\underset{frictionforces}{\underbrace{\sum _{i=1}^{N_c}{\mathit{F}}_i^{fric}}}+\underset{damping}{\underbrace{{\mathbf{F}}_p^{damp}}}+\underset{gravity}{\underbrace{{\mathbf{F}}_p^{grav}}},$ (2)

where $p=\mathrm{1,}\mathrm{2,}\mathrm{...,}{N}_p$, $N_p$ is the number of particles, $N_c$ is the number of fabric nodes in contact with the projectile, $m_p$ is the mass of a single projectile, $r_p$ is the position of projectile $p$, $F_p^{tot}$ is the total force of the projectile. Note that the contact forces $F_i^{con}$ and $F_i^{fric}$ presented in equation (2) is the same used in equation (1) for the correspondent lumped mass $i$ in contact with the projectile $p$.

FORMULATIONS FOR YARN-SEGMENT NETWORK

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

$\psi =\frac{1}{2}\mathbb{D}{E}^2,$ (3)

where $\mathbb{D}$ is the Young’s modulus and $E$ is the Green-Lagrange strain. Using equation (125), the second Piola-Kirchhoff stress $S$ is given by

$S=\frac{\partial \psi }{\partial E}=\mathbb{D}E.$ (4)

It is convenient to work with quantities expressed in terms of the stretch ratio $U$, that can be obtained as

$U=\frac{L}{L^r},$ (5)

where $L$ is the deformed length of the fibril and $L^r$ is its original length. The Green-Lagrange strain is obtained from (119) as following

$E=\frac{1}{2}\left[{\left(\frac{L}{L^r}\right)}^2-1\right].$ (6)

For this relaxed one-dimensional model, the Jacobian of $F$ is

$J=F$ (7)

and the stresses defined in equation (116) can be rewritten by

$\frac{\Sigma }{J}=P=JS=\sigma .$ (8)

When $F\le 1$ the yarn has compression and when $F>1$ it has tensile stress. The compression is neglected, so, when $F\le 1$ the first Piola-Kirchhoff stresses is enforced to be zero, i.e., $P=0$. The equations (118), (120) and (7) results in

$F=J=U.$ (9)

The first Piola-Kirchhoff stress $P$ can be expressed by

$P=\frac{F^{yarn}}{A^r},$ (10)

where $F^{yarn}$ is the force of the yarn on the undeformed cross-sectional area $A^r$, in the reference configuration. Using the relations presented in equations (5), (8), (9) and (10) the force carried in the yarn-segment results in

$F^{yarn}=US{A}^r=\frac{L}{L^r}S{A}^r.$ (11)

Substituting (4) and (6) into (11) lies to

$F^{yarn}=\frac{L}{2L^r}\mathbb{D}{A}^r\left[{\left(\frac{L}{L^r}\right)}^2-1\right].$ (12)

So, the yarn force presented in equation (1) can be defined by

$F_I^{yarn}=\frac{L_I}{2L_I^r}{\mathbb{D}}_I{A}_I^r\left[{\left(\frac{L_I}{L_I^r}\right)}^2-1\right]{\alpha }_I,$ (13)

where

${\alpha }_I=\frac{r_I^{+}-r_I^{-}}{r_I^{+}-r_I^{-}},$ (14)

$r_I^{+}$ denotes the position vector of the endpoint connected to the lumped mass, $r_I^{-}$ denotes the endpoint that is connected to its neighboring mass and $\parallel \cdot \parallel$ indicates the Euclidean norm in ${\mathcal{R}}^3$. See the geometrical interpretation of the yarn force in Figure 3.

Figure 3
Geometrical interpretation of the force at a yarn-segment.

Note that the deformed length of the fibril is calculated by

$L=\parallel r_I^{+}-r_I^{-}\parallel .$ (15)

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

$L^{yarn}=L^{coating}=L,$ (16)

${\left(L^r\right)}^{yarn}={\left(L^r\right)}^{coating}=L^r,$ (17)

$U^{yarn}=U^{coating}=U.$ (18)

Observing the yarn force presented in equation (13), the mechanical forces for the yarn and coating can be defined, respectively, by

$F_I^{yarn}=\frac{L_I}{2L_I^r}{\left({\mathbb{D}}_I\right)}^{yarn}{\left(A_I^r\right)}^{yarn}\left[{\left(\frac{L_I}{L_I^r}\right)}^2-1\right]{\alpha }_I,$ (19)

$F_I^{coating}=\frac{L_I}{2L_I^r}{\left({\mathbb{D}}_I\right)}^{coating}{\left(A_I^r\right)}^{coating}\left[{\left(\frac{L_I}{L_I^r}\right)}^2-1\right]{\alpha }_I.$ (20)

A simple estimate of the increase in stiffness due to the coating can be determined by summing the forces

$F_I^{yarn}+F_I^{coating}=\frac{L_I}{2L_I^r}\left[{\left(\frac{L_I}{L_I^r}\right)}^2-1\right]\left\{{\left({\mathbb{D}}_I\right)}^{yarn}{\left(A_I^r\right)}^{yarn}+{\left({\mathbb{D}}_I\right)}^{coating}{\left(A_I^r\right)}^{coating}\right\}{\alpha }_I.$ (21)

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 $R_{fibril}$. Let be $N_{fibril}$ the number of microfibril in a yarn-segment, $R_{yarn}$ the radii of initial cross section of the yarn-segment and $\mathcal{T}$ the thickness of the coating. So, the cross-section area of the yarn and coating are, respectively,

${\left(A_I^r\right)}^{yarn}=\pi {\left(R_{yarn}\right)}^2=N_{fibril}\pi {\left(R_{fibril}\right)}^2,$ (22)

${\left(A_I^r\right)}^{coating}=\pi \left[{\left(R_{yarn}+\mathcal{T}\right)}^2-{\left(R_{yarn}\right)}^2\right],$ (23)

where

$R_{yarn}=\sqrt{N_{fibril}}{R}_{fibril}.$ (24)

For Zylon each yarn approximately contains $N_{fibril}=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 ${\alpha }_I^{yarn}$ such that it can construct a new Young’s modulus in each time as following

${\mathbb{D}}_I^{yarn}\left(t\right)={\alpha }_I^{yarn}\left(t\right){\mathbb{D}}_{0I}^{yarn},$ (25)

where

$0<{\alpha }_I^{yarn}\le 1$ (26)

and ${\mathbb{D}}_{0I}^{yarn}$ is the Young modulus in the initial of the analysis, i.e., ${\mathbb{D}}_{0I}^{yarn}={\mathbb{D}}_I^{yarn}\left(t=0\right)$. 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, ${\alpha }_I^{yarn}=1$, while for a yarn that is completely damaged, ${\alpha }_I^{yarn}=0$. Hence, in the beginning of the analysis ${\alpha }_I^{yarn}\left(t=0\right)=1$. Note that when the parameter ${\alpha }_I^{yarn}$ goes to zero the physical material parameter ${\mathbb{D}}_I^{yarn}$ also goes to zero.

The damage parameter proposed by Zohdi (2002b) is governed by the following over-stress evolution equation

${\dot{\alpha }}_I^{yarn}=\zeta \left(\frac{{\sigma }_I}{{\sigma }_{crit}}-1\right){\alpha }_I^{yarn},$ (27)

where

$\zeta =\{\begin{array}{cc}0& ,if{\sigma }_I<{\sigma }_{crit}\\ {\zeta }^{*}<0& ,if{\sigma }_I\ge {\sigma }_{crit}\end{array},$ (28)

${\sigma }_I$ is the Cauchy stress, ${\sigma }_{crit}$ is the Cauchy stress limit to starts the damage of the fibers and $\zeta$ 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 $U_I\le 1$ then ${\sigma }_I=0$.

For Zylon, according to the manufacturer, for a fibril, $U_{crit}=1.03$ and the Young’s modulus is ${\mathbb{D}}^{fibril}=180GPa$. From equation (121) the Green-Lagrange strain tensor is $E_{crit}=\frac{1}{2}\left({1.03}^2-1\right)=0.03045$. From equations (4), (8) and (9), the Cauchy stress tensor lies to ${\sigma }_{crit}=U\mathbb{D}E=1.03\cdot 180000\cdot 0.03045=5645.43MPa$. The rupture of the material Zylon is presented by Zohdi and Powell (2006) and is supposed to occur when the stretch is around $U_{rupture}\cong 1.06$, see Figure 4.

The damage parameter was also proposed by Zohdi (2014c) as shown in the following equation

${\dot{\alpha }}_I^{yarn}=\mathcal{P}\left(\frac{\left|{\sigma }_I\right|}{{\sigma }_{crit}}-1\right),$ (29)

where the parameter $\mathcal{P}$ assumes the following values

$\mathcal{P}=\{\begin{array}{cc}0& if\left|{\sigma }_I\right|\le {\sigma }_{crit}\\ {\mathcal{P}}^{*}<0& if\left|{\sigma }_I\right|>{\sigma }_{crit}\end{array}.$ (30)

The equation (29) can also be rewritten as

${\alpha }_I^{yarn}\left(t_{i+1}\right)={\alpha }_I^{yarn}\left(t_i\right)+P\left(\frac{\left|{\sigma }_I\right|}{{\sigma }_{crit}}-1\right)\Delta t,$ (31)

where $\Delta t$ is the increment of time. The equations (27) and (31) are valid when the stretch is equal or higher than $U_{crit}$. 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

${\alpha }_I^{yarn}\left(t\right)=\min \{\begin{array}{c}1\\ {\alpha }_I\left(0\le t^{current}\le t\right)\\ \frac{e^{-\lambda \left[U_I\left(t\right)-U_{I,crit}\right]}-e^{-0.03\lambda }}{1-e^{-0.03\lambda }}\end{array},$ (32)

where $U_I\left(t\right)$ is the stretch of the yarn-segment $I$ at time $t$, $U_{I,crit}$ is the critical stretch of the yarn-segment and $\lambda \ge 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 $\lambda \approx 124$. The equation (32) has a numerical instability problem to represent the case of no damage generated, $\lambda \to 0$, and sudden rupture, $\lambda \to \infty$.

Basically, a parameter ${\alpha }_I^{yarn}\left(t\right)$ is created to describe the failure of the yarn-segment by checking whether a critical stretch has been attained or exceeded, i.e., $U\left(t\right)\ge U_{crit}$ 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

${\alpha }_I^{yarn}\left(t\right)=\min \{\begin{array}{c}1\\ {\alpha }_I\left(0\le t^{current}\le t\right)\\ e^{-\lambda \left[\frac{U_I\left(t\right)-U_{I,crit}}{U_{I,crit}}\right]}\end{array}.$ (33)

Analyzing equation (33) is concluded that when $\lambda \to \infty$ the type of failure tends towards sudden rupture and when $\lambda \to 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 $\lambda \approx 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 $\lambda =124$ and $\lambda =135$, respectively. It is important to mention that in practical procedures, is suggested to consider the total damage when ${\alpha }_I^{yarn}$ approaches $0.02$, see Figure 4.

Figure 4
Behavior of the damage parameter.

The equations (27), (32) and (33) indicate that damage is irreversible, i.e., ${\alpha }_I^{yarn}$ 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

$F=F^e{F}^p\Rightarrow U=U^e{U}^p,$ (34)

where $U^e$ and $U^p$ is sometimes referred as an elastic and plastic stretch strain, respectively.

The Green-Lagrange strain is defined by

$E=E^e+E^p,$ (35)

where $E^e$ and $E^p$ 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

$S^e=\mathbb{D}{E}^e.$ (36)

Substituting (35) into (36) results in

$S^e=\mathbb{D}\left(E-E^p\right).$ (37)

Introducing equation (121) into (37), result in

$S^e=\mathbb{D}\left[\frac{1}{2}\left(U^2-1\right)-\frac{1}{2}\left(U^p{}^2-1\right)\right].$ (38)

The above equation can be rewritten as

$S^e=\frac{1}{2}\mathbb{D}\left(U^2-U^p{}^2\right).$ (39)

From equations (8) and (9) is possible to determine the Cauchy stress tensor

$\sigma =U^e{S}^e=\frac{1}{2}\mathbb{D}{U}^e\left(U^2-U^p{}^2\right).$ (40)

The equation (40) can also be written as

$\sigma =U^e{S}^e=\frac{1}{2}\mathbb{D}{U}^e\left[{\left(U^e\right)}^2-1\right],$ (41)

where $U^e=U{\left(U^p\right)}^{-1}$.

The variation of the plastic stretch strain can be computed by the following relation

${\dot{U}}^p=\{\begin{array}{cc}0& if\left|\sigma \right|\le {\sigma }_y\\ \mathcal{P}\left(\frac{\left|\sigma \right|}{{\sigma }_y}-1\right)& if\left|\sigma \right|>{\sigma }_{y}and\sigma >0\\ -\mathcal{P}\left(\frac{\left|\sigma \right|}{{\sigma }_y}-1\right)& if\left|\sigma \right|>{\sigma }_{y}and\sigma <0\end{array},$ (42)

where

${\sigma }_y={\sigma }_{yo}+\mathcal{H}\left(U^p\right),$ (43)

$\mathcal{P}$ is a rate parameter, ${\sigma }_y$ is the yield stress (connected to stretch sensitivity), ${\sigma }_{yo}$ is the initial yield stress and $\mathcal{H}$ is a hardening variable.

It is important to mention that when $\left|\sigma \right|\le {\sigma }_y$ the plastic strain is zero, i.e., $E^p=0$. In this case, the plastic stretch tensor is $U^p=1$ and the second Piola-Kirchhoff stress tensor is $S=\mathbb{D}E$, where $E=E^e$. 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

$U^p\left(t_{i+1}\right)=U^p\left(t_i\right)+\Delta U^p\left(t_i\right),$ (44)

where

$\Delta U^p\left(t_i\right)=\{\begin{array}{cc}0& if\left|\sigma \right|\le {\sigma }_y\\ \mathcal{P}\left(\frac{\left|\sigma \right|}{{\sigma }_y}-1\right)\Delta t& if\left|\sigma \right|>{\sigma }_{y}and\sigma >0\\ -\mathcal{P}\left(\frac{\left|\sigma \right|}{{\sigma }_y}-1\right)\Delta t& if\left|\sigma \right|>{\sigma }_{y}and\sigma <0\end{array},$ (45)

and $\Delta t$ 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 $\Delta U^p$ is necessary to know the Cauchy stress tensor $\sigma$, that is unknown in the current time step and need to be calculated by equation (40). So, to solve it, the value of $U^p$ 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 $U^p=1$.

1.1.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} .

Algorithm 1
Elastoplasticity with damage.

1)	Initialize the variables: $U_i^p\left(0\right)=1$, ${\alpha }_i\left(0\right)=1$, $L^r$, $\mathbb{D}\left({\alpha }_i\right)$.
2)	Calculate the known variables: $L$, $U_{i+1}$.
3)	Compute the elastic part of the stretch: $U_t^e=U_{i+1}{\left(U_i^p\right)}^{-1}$.
4)	Compute the elastic part of the Green-Lagrange strain: $E_t^e=\frac{1}{2}\left[{\left(U_t^e\right)}^2-1\right]$.
5)	Compute the elastic part of the second Piola-Kirchhoff stress: $S_t^e=\mathbb{D}\left({\alpha }_i\right)E_t^e$.
6)	Compute elastic part of the Cauchy stress: ${\sigma }_t^e=U_t^e{S}_t^e$.
7)	Compute: $F_t={\sigma }_t^e-{\sigma }_y$.
8)	Compute the variation of the plastic part of the stretch tensor $\Delta U^p$ from equation (45).
9)	Update the plastic stretch tensor: IF $F_t<0$ THEN $U_{i+1}^p=U_i^p$ ELSE IF $U_{i+1}^p=U_i^p+\Delta U^p$.
10)	Compute the elastic part of the stretch: $U_t^e=U_{i+1}{\left(U_{i+1}^p\right)}^{-1}$.
11)	Compute the elastic part of the Green-Lagrange strain: $E_t^e=\frac{1}{2}\left[{\left(U_t^e\right)}^2-1\right]$.
12)	Compute the second Piola-Kirchhoff stress: $S_{i+1}=\mathbb{D}\left({\alpha }_i\right)E_t^e$.
13)	Compute the Cauchy stress: $\sigma =U_t^e{S}_{i+1}$.

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 $L=0.00254m$, Young modulus $\mathbb{D}=180GPa$ and critical stretch $U_{crit}=1.03$. As explained before, the Cauchy stress limit lies to $\sigma =5645.43MPa$ and the rupture occurs when $U_{rupture}\cong 1.06$. In order to generate a tensile force, in each time increment of $\Delta t={10}^{-6}s$ the yarn segment was increased its length in $\Delta L=2.03\cdot {10}^{-7}m$. The case that have no damage in the yarn segment $\left(\lambda =0\right)$ and have only elastoplasticity with parameter $\mathcal{P}={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., $\sigma \le {\sigma }_y$.

Figure 5
Graphic of the elastoplastic constitutive law: ($\lambda =0$ and $\mathcal{P}={10}^4$).

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) to (21), lies to

$F_I^{yarn}=\frac{L_I}{2L_I^r}\left[{\left(\frac{L_I}{L_I^r}\right)}^2-{\left(U_I^p{}^{yarn}\right)}^2\right]\left\{{\alpha }_I^{yarn}{\mathbb{E}}_{0I}^{yarn}{\left(A_I^r\right)}^{yarn}\right\}{\alpha }_I$ (46)

$F_I^{coating}=\frac{L_I}{2L_I^r}\left[{\left(\frac{L_I}{L_I^r}\right)}^2-{\left(U_I^p{}^{coating}\right)}^2\right]\left\{{\alpha }_I^{coating}{\left({\mathbb{E}}_{0I}\right)}^{coating}{\left(A_I^r\right)}^{coating}\right\}{\alpha }_I,$ (47)

where ${\alpha }_I$, ${\left(A_I^r\right)}^{yarn}$ and ${\left(A_I^r\right)}^{coating}$ are defined in (14), (22) and (23), respectively.

The Cauchy stress tensor for the yarn and coating are respectively

${\sigma }_I^{yarn}=\frac{L_I}{2L_I^r}\left[{\left(\frac{L_I}{L_I^r}\right)}^2-{\left(U_I^p{}^{yarn}\right)}^2\right]\left\{{\alpha }_I^{yarn}{\left({\mathbb{E}}_{0I}\right)}^{yarn}\right\},$ (48)

${\sigma }_I^{coating}=\frac{L_I}{2L_I^r}\left[{\left(\frac{L_I}{L_I^r}\right)}^2-{\left(U_I^p{}^{coating}\right)}^2\right]\left\{{\alpha }_I^{coating}{\left({\mathbb{E}}_{0I}\right)}^{coating}\right\}.$ (49)

Note that $U_I^p{}^{yarn}$ and $U_I^p{}^{coating}$ 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 $\left(\lambda =135\right)$ and elastoplasticity with parameter $\mathcal{P}={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 $\lambda =135$ and $\mathcal{P}=0$.

Figure 6
Graphic of the elastoplastic constitutive law with rather gradual damage ($\lambda =135$ and $\mathcal{P}={10}^4$ or $\mathcal{P}=0$).

This formulation presented in section 4.4 is more stable than the others ones presented. To simulate a sudden rupture, the parameter $\lambda \to \infty$. In this numerical example it should be $\lambda ={10}^4$, see Figure 7.

Figure 7
Graphic of the elastoplastic constitutive law with sudden damage ($\lambda ={10}^4$ and $\mathcal{P}={10}^4$).

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.

Figure 8
Interpenetration of a fabric node (picture based on Zohdi (2014c) and Zohdi (2010a)).

In Figure 8, $r_i^f\left(t\right)$ is the position of the converged value from the previous time step, $r_i^{f,pred}\left(t+\Delta t\right)$ is the predicted contact position in current time step, $r_i^{f,corr}\left(t+\Delta t\right)$ is the corrected contact position in current time step, $r^{p,cm}\left(t+\Delta t\right)$ is the position of the center of mass of the projectile in current time step, $r^{p,cp}\left(t+\Delta t\right)$ is the position of the closest point on the surface of the contacting object in current time step and $R$ is the radii of the projectile.

Geometrical interpretation for a fabric node

Assuming interpenetration for a node $i$, see Figure 8, the contact force on a node $i$ can be calculate as

$F_i^{con}=\frac{m_i}{\Delta t^2}\left[r_i^{f,corr}\left(t+\Delta t\right)-r_i^{f,pred}\left(t+\Delta t\right)\right],$ (50)

where $m_i$ is the lamped mass of the node $i$ of the fabric. From above equation, the corrected contact position is given by

$r_i^{f,corr}\left(t+\Delta t\right)=r_i^{f,pred}\left(t+\Delta t\right)+\frac{\Delta t^2}{m_i}{F}_i^{con},$ (51)

where the predicted contact position is defined by

$r_i^{f,pred}\left(t+\Delta t\right)=r_i^f\left(t\right)+\Delta t\left[\varphi v_i^{f,pred}\left(t+\Delta t\right)+\left(1-\varphi \right)v_i^f\left(t\right)\right],$ (52)

where $v_i^{f,pred}\left(t+\Delta t\right)$ is the velocity of the node $i$ in time step $t+\Delta t$, $v_i^f\left(t\right)$ is the velocity of the converged value from the previous time step $\left(t\right)$ and $\varphi$ is the trapezoidal time-stepping parameter. The velocity of the predicted contact position is calculated by

$v_i^{f,pred}\left(t+\Delta t\right)=v_i^f\left(t\right)+\frac{\Delta t}{m_i}\left[\varphi F_i^{f,pred,oth}\left(t+\Delta t\right)+\left(1-\varphi \right)F_i^{f,oth}\left(t\right)\right],$ (53)

where $F_i^{f,oth}\left(t\right)$ are all forces other than the contact forces. In this work, the forces other than contact are considered to be

$\underset{\begin{array}{c}otherforces\\ thancontact\end{array}}{\underbrace{{\mathit{F}}_i^{f,oth}}}=\underset{surroundingyarn}{\underbrace{\sum _{I=1}^4{\mathbf{F}}_I^{yarn}}}+\underset{surroundingcoating}{\underbrace{\sum _{I=1}^4{\mathbf{F}}_I^{coating}}}+\underset{damping}{\underbrace{{\mathbf{F}}_i^{damp}}}+\underset{gravity}{\underbrace{{\mathbf{F}}_i^{grav}}}.$ (54)

The equation (52) can be rewritten as

$r_i^{f,pred}\left(t+\Delta t\right)=r_i^f\left(t\right)+\Delta t{v}_i^f\left(t\right)+\varphi \Delta t\left[v_i^{f,pred}\left(t+\Delta t\right)-v_i^f\left(t\right)\right].$ (55)

Substituting (53) into (55), lies to

$r_i^{f,pred}\left(t+\Delta t\right)=r_i^f\left(t\right)+\Delta t{v}_i^f\left(t\right)+\frac{\varphi {\left(\Delta t\right)}^2}{m_i}\left[\varphi F_i^{f,pred,oth}\left(t+\Delta t\right)+\left(1-\varphi \right)F_i^{f,oth}\left(t\right)\right].$ (56)

Procedure to calculate the Fabric motion

Figure 9
Geometrical interpenetration of a fabric node update.

To updated each fabric node $i$, basically, the projectile is frozen and a new position of the fabric node $i$ 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.

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 $r_i^{f,pred}\left(t+\Delta t\right)$ 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

$r_i^{f,corr}\left(t+\Delta t\right)=r_i^{p,cp}\left(t+\Delta t\right),$ (57)

where $r_i^{p,cp}\left(t+\Delta t\right)$ 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 $v_i^{f,corr}\left(t+\Delta t\right)$ to that of the closest point on the body $v_i^{p,cp}\left(t+\Delta t\right)$, for node $i$, as

$v_i^{f,corr}\left(t+\Delta t\right)=v_i^{p,cp}\left(t+\Delta t\right).$ (58)

The friction force is calculated using

$F_i^{fric}=m_i{\dot{v}}_i-F_i^{con}-F_i^{f,oth},$ (59)

where the contact force $F_i^{con}$ is computed by equation (50) and ${\dot{v}}_i$ can be considered to be the acceleration of the projectile, since it is considered “stick” case. The acceleration ${\dot{v}}_i$ can be computed by equation (83).

If the friction force exceeds the theoretical limit

$\Vert F_i^{fric}\Vert >{\mu }_s\Vert F_i^{con}\Vert ,$ (60)

then the assumption that there is “stick” is incorrect, and sliding must occur, resulting in a nonmatching velocity of

$v_i^{f,corr}\left(t+\Delta t\right)=v_i^{f,pred}\left(t+\Delta t\right)+\frac{\Delta t}{m_i}\left[F_i^{con}+F_i^{fric}\right],$ (61)

with the (frictional) projection

$F_i^{fric}={\mu }_d\Vert F_i^{con}\Vert \frac{v_i^{p,cp}-v_i^{f,cp}}{\Vert v_i^{p,cp}-v_i^{f,cp}\Vert },$ (62)

where the velocity of the node $i$ of the fabric at the closest contact point $v_i^{f,cp}$ is assumed to be

$v_i^{f,cp}=v_i^{f,corr}\left(t+\Delta t\right).$ (63)

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)

$v_p\left(t+\Delta t\right)=v_p\left(t\right)-\frac{\Delta t}{m_p}{\displaystyle \sum }_{i=1}^{N_c}\left(F_i^{con}+F_i^{fric}\right)+\frac{\Delta t}{m_p}\left[\varphi F_p^{oth}\left(t+\Delta t\right)+\left(1-\varphi \right)F_p^{oth}\left(t\right)\right],$ (64)

where $N_c$ is the number of nodes in contact and

$r_p\left(t+\Delta t\right)=r_p\left(t\right)+\Delta t\left[\varphi v_p\left(t+\Delta t\right)+\left(1-\varphi \right)v_p\left(t\right)\right].$ (65)

In this work, these forces other than contact are considered to be

$\underset{\begin{array}{c}otherforces\\ thancontact\end{array}}{\underbrace{{\mathit{F}}_p^{oth}}}=\underset{damping}{\underbrace{{\mathbf{F}}_p^{damp}}}+\underset{gravity}{\underbrace{{\mathbf{F}}_p^{grav}}}.$ (66)

The equation (65) can be rewritten by

$r_p\left(t+\Delta t\right)=r_p\left(t\right)+v_p\left(t\right)\Delta t+\Delta t\varphi \left[v_p\left(t+\Delta t\right)-v_p\left(t\right)\right].$ (67)

Substituting (64) into (67) lies to

$r_p\left(t+\Delta t\right)=r_p\left(t\right)+v_p\left(t\right)\Delta t+\frac{\varphi {\left(\Delta t\right)}^2}{m_p}\left\{-{\displaystyle \sum }_{i=1}^{N_c}\left(F_i^{con}+F_i^{fric}\right)+\left[\varphi F_p^{oth}\left(t+\Delta t\right)+\left(1-\varphi \right)F_p^{oth}\left(t\right)\right]\right\}.$ (68)

The velocity of any point $j$ on the surface of the body can be determined by simply calculating

$v_j=v_{p,cm}+{\omega }_p\times r_{p,cm\to j}.$ (69)

The projectile’s angular velocity and rotation are determined in a similar manner by integrating the equations for an angular momentum balance

${\dot{H}}_p=\frac{d}{dt}\left({\overline{I}}_{p,s}\cdot {\omega }_p\right)=M_p^{tot},$ (70)

where ${\overline{I}}_{p,s}$ is the mass moment of the projectile, ${\omega }_p$ is the angular velocity and $M_p^{tot}$ 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):

$\underset{t_1}{\overset{t_2}{{\displaystyle \int }}}d{\omega }_p=\underset{t_1}{\overset{t_2}{{\displaystyle \int }}}\frac{1}{{\overline{I}}_{p,s}}{M}_p^{tot}dt\Rightarrow {\omega }_p\left(t_2\right)={\omega }_p\left(t_1\right)+\frac{\Delta t}{{\overline{I}}_{p,s}}{M}_p^{tot}.$ (71)

Then, based in equation (71), the time discretization is done by the following equation:

${\omega }_p\left(t+\Delta t\right)={\omega }_p\left(t\right)+\frac{\Delta t}{{\overline{I}}_{p,s}}\left\{\varnothing M_p^{tot}\left(t\right)+\left[\varnothing M_p^{tot}\left(t+\Delta t\right)-\varnothing M_p^{tot}\left(t\right)\right]\right\},$ (72)

where $\varnothing$ is a scalar value with range $0<\varphi \le 1$. It is important to mention that when the parameter $\varphi =0$ the system is solved by explicit method and for $0<\varphi \le 1$ the system is solved by implicit method. The equation (72) can be rewritten as

${\omega }_p\left(t+\Delta t\right)={\omega }_p\left(t\right)+\frac{\Delta t}{{\overline{I}}_{p,s}}\left[\varnothing M_p^{tot}\left(t+\Delta t\right)+\left(1-\varnothing \right)M_p^{tot}\left(t\right)\right].$ (73)

The moment $M_p^{tot}$ can also be calculated by

$M_p^{tot}={\displaystyle \sum }_{i=1}^{N_c}{r}_{p\to c}\times F_i^{fric},$ (74)

where $F_i^{fric}$ is the friction force between the nodes $i$ 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 ${\theta }^p$ of a particle is obtained by

${\omega }_p=\frac{d{\theta }_p}{dt}\Rightarrow \underset{t_1}{\overset{t_2}{{\displaystyle \int }}}d{\theta }_p=\underset{t_1}{\overset{t_2}{{\displaystyle \int }}}{\omega }_{p}dt\Rightarrow {\theta }_p\left(t_2\right)={\theta }_p\left(t_1\right)+{\omega }_p\Delta t.$ (75)

The time discretization is done by the following equation:

${\theta }_p\left(t+\Delta t\right)={\theta }_p\left(t\right)+\Delta t\left\{\varnothing {\omega }_p\left(t+\Delta t\right)+\left(1-\varnothing \right){\omega }_p\left(t\right)\right\}.$ (76)

Finally, for a point $P^{\prime }$ belonging to a projectile in the initial configuration placed ${r_{p\to c}^{P^{\prime }}|}_{t=0}$ from its center, its current position is defined by:

$r_p^{P^{\prime }}\left(t+\Delta t\right)=r_p\left(t+\Delta t\right)+{r_{p\to c}^{P^{\prime }}|}_{t=0}\times {\theta }_p\left(t+\Delta t\right).$ (77)

Based on equation (64), the dynamics of the projectile is given by

$m_p{\ddot{r}}_p=F_p^{tot}=-{\displaystyle \sum }_{i=1}^{N_c}\left(F_i^{con}+F_i^{fric}\right)+\left[\varphi F_p^{oth}\left(t+\Delta t\right)+\left(1-\varphi \right)F_p^{oth}\left(t\right)\right].$ (78)

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.

OTHER MECHANICAL FORCES

Damping force

The other force that can be taken into account is the damping effect given by

$F_i^{damp}=c_i{\dot{r}}_i,$ (79)

where $c_i{\dot{r}}_i$ 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 ($-z$ direction) and is given by

$F_i^{grav}=-m_{i}g\widehat{z},$ (80)

where $g$ is the gravity acceleration and $\widehat{z}$ is the unit vector in the $z$ direction.

NUMERICAL SOLUTION SCHEME

In order to describe the overall time-stepping scheme, is considered the dynamics of a single $i^{th}$ lumped mass as presented in equation (1) and the dynamics of a single $p^{th}$ 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

$r={\left[\begin{array}{cccccccc}{r}_1^f& r_2^f& \dots & r_{N_f}^f& r_1^p& r_2^p& \cdots & r_{N_p}^p\end{array}\right]}^T={\left[\begin{array}{cc}{r}_i& r_p\end{array}\right]}^T$ (81)

and the velocity vector defined by

$v={\left[\begin{array}{cccccccc}{v}_1^f& v_2^f& \dots & v_{N_f}^f& v_1^p& v_2^p& \cdots & v_{N_p}^p\end{array}\right]}^T={\left[\begin{array}{cc}{v}_i& v_p\end{array}\right]}^T,$ (82)

where $i=\mathrm{1,}\mathrm{2,}\cdots ,N_f$, $N_f$ is the number of lumped mass, $p=\mathrm{1,}\mathrm{2,}\cdots ,N_p$, $N_p$ 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\le \varphi \le 1$, the following relations are valid:

${\dot{v}}_i\left(t+\varphi \Delta t\right)=\frac{v_i\left(t+\Delta t\right)-v_i\left(t\right)}{\Delta t};$ (83)

$v_i\left(t+\Delta t\right)=v_i\left(t\right)+\frac{1}{m_i}\underset{t}{\overset{t+\Delta t}{{\displaystyle \int }}}{F}_i^{tot}dt=v_i\left(t\right)+\frac{\Delta t}{m_i}\left[\varphi F_i^{tot}\left(t+\Delta t\right)+\left(1-\varphi \right)F_i^{tot}\left(t\right)\right].$ (84)

The position can be computed by applying the mid-point rule, which results in

$r_i\left(t+\Delta t\right)=r_i\left(t\right)+\Delta t{v}_i\left(t+\varphi \Delta t\right)\approx r_i\left(t\right)+\Delta t\left[\varphi v_i\left(t+\Delta t\right)+\left(1-\varphi \right)v_i\left(t\right)\right].$ (85)

Substituting equation (84) into (85) is obtained

$r_i\left(t+\Delta t\right)=r_i\left(t\right)+v_i\left(t\right)\Delta t+\frac{\varphi {\left(\Delta t\right)}^2}{m_i}\left[\varphi F_i^{tot}\left(t+\Delta t\right)+\left(1-\varphi \right)F_i^{tot}\left(t\right)\right].$ (86)

The equation (86) can be solved recursively by recasting the relation as

$r_i^{L+1}=r_i^L+v_i^L\Delta t+\frac{\varphi {\left(\Delta t\right)}^2}{m_i}\left[\varphi F_i^{tot,L+1}+\left(1-\varphi \right)F_i^{tot,L}\right],$ (87)

and the same for equation (84) as

$v_i^{L+1}=v_i^L+\frac{\Delta t}{m_i}\left[\varphi F_i^{tot,L+1}+\left(1-\varphi \right)F_i^{tot,L}\right],$ (88)

where superscript $L$ and $L+1$ indicates the indices of the current and the future time-step in the discretization, respectively. Basically, the superscript $L$ is a time interval counter. This is a generalized one-step scheme, with $\varphi \in \left(\mathrm{0,1}\right]$, in which, when $\varphi =0$ is the explicit Euler integration scheme, while $\varphi =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 $r_i^{L+1}$ and $v_i^{L+1}$.

The set of equations represented by (87) and (88) can be solved recursively by recasting the relation as

$r_i^{L+\mathrm{1,}K}=r_i^L+v_i^L\Delta t+\frac{\varphi {\left(\Delta t\right)}^2}{m_i}\left[\varphi F_i^{tot,L+\mathrm{1,}K-1}+\left(1-\varphi \right)F_i^{tot,L}\right]$ (89)

$v_i^{L+\mathrm{1,}K}=v_i^L+\frac{\Delta t}{m_i}\left[\varphi F_i^{tot,L+\mathrm{1,}K-1}+\left(1-\varphi \right)F_i^{tot,L}\right],$ (90)

where $K=\mathrm{1,2,3}\dots$ is the index of iteration within time-step $L+1$. These equations can also be represented in the following form:

$r_i^{L+\mathrm{1,}K}={\mathcal{G}}_r\left(r_i^{L+\mathrm{1,}K-1}\right)+{\mathcal{R}}_{i,r};$ (91)

$v_i^{L+\mathrm{1,}K}={\mathcal{G}}_v\left(v_i^{L+\mathrm{1,}K-1}\right)+{\mathcal{R}}_{i,v}.$ (92)

where,

${\mathbf{F}}_i^{tot,L+1,K-1}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\mathbf{F}}_i^{tot,L+1,K-1}\left({\mathit{r}}_1^{L+1,K-1},{\mathit{r}}_2^{L+1,K-1},\dots ,{\mathit{r}}_{N_p}^{L+1,K-1}\right),$ (93)

${\mathbf{F}}_i^{tot,L}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\mathbf{F}}_i^{tot,L}\left({\mathit{r}}_1^L,{\mathit{r}}_2^L,\dots ,{\mathit{r}}_{N_p}^L\right),$ (94)

${\mathcal{G}}_r\left(r_i^{L+\mathrm{1,}K-1}\right)=\frac{{\left(\varphi \Delta t\right)}^2}{m_i}{F}_i^{tot,L+\mathrm{1,}K-1},$ (95)

${\mathcal{R}}_{i,r}=r_i^L+v_i^L\Delta t+\frac{\varphi {\left(\Delta t\right)}^2}{m_i}\left(1-\varphi \right)F_i^{tot,L},$ (96)

${\mathcal{G}}_v\left(v_i^{L+\mathrm{1,}K-1}\right)=\frac{\varphi \Delta t}{m_i}{F}_i^{tot,L+\mathrm{1,}K-1},$ (97)

${\mathcal{R}}_{i,v}=v_i^L+\frac{\Delta t}{m_i}\left(1-\varphi \right)F_i^{tot,L}.$ (98)

The terms ${\mathcal{R}}_{i,r}$ and ${\mathcal{R}}_{i,v}$ are remainder terms that do not depend on the solution. The convergence of such a scheme is dependent on the behavior of ${\mathcal{G}}_r$ and ${\mathcal{G}}_v$. The form of the error function is presented below as

${\mathit{\varpi }}_{i,r}^{L+1,K}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\mathit{r}}_i^{L+1,K}-{\mathit{r}}_i^{L+1},$ (99)

${\mathit{\varpi }}_{i,v}^{L+1,K}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\mathit{v}}_i^{L+1,K}-{\mathit{v}}_i^{L+1}.$ (100)

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

$\Vert r_i^{L+\mathrm{1,}K}-r_i^{L+1}\Vert =\Vert {\mathcal{G}}_r\left(r_i^{L+\mathrm{1,}K-1}\right)-{\mathcal{G}}_r\left(r_i^{L+1}\right)\Vert \le {\eta }_r^{L+\mathrm{1,}K}\Vert r_i^{L+\mathrm{1,}K-1}-r_i^{L+1}\Vert ,$ (101)

$\Vert v_i^{L+\mathrm{1,}K}-v_i^{L+1}\Vert =\Vert {\mathcal{G}}_v\left(v_i^{L+\mathrm{1,}K-1}\right)-{\mathcal{G}}_v\left(v_i^{L+1}\right)\Vert \le {\eta }_v^{L+\mathrm{1,}K}\Vert v_i^{L+\mathrm{1,}K-1}-v_i^{L+1}\Vert ,$ (102)

where $0\le {\eta }_r^{L+\mathrm{1,}K}<1$ and $0\le {\eta }_v^{L+\mathrm{1,}K}<1$ for each iteration $K$. Then ${\varpi }_{i,r}^{L+\mathrm{1,}K}\to 0$ and ${\varpi }_{i,v}^{L+\mathrm{1,}K}\to 0$ for any arbitrary starting value $r_i^{L+\mathrm{1,}K=0}$ and $v_i^{L+\mathrm{1,}K=0}$, as $K\to \infty$, which is a contraction condition that is sufficient, but not necessary, for convergence. The convergence of equation (89) is scaled by ${\eta }_r\propto \frac{{\left(\varphi \Delta t\right)}^2}{m_i}$ (see equation (95)) and the convergence of equation (90) is scaled by ${\eta }_v\propto \frac{\varphi \Delta t}{m_i}$ (see equation (97)). Therefore, the contraction constant of ${\mathcal{G}}_r$ and ${\mathcal{G}}_v$ is directly dependent on the magnitude of the interaction forces ($\Vert F_i\Vert$), inversely proportional to the masses $m_i$ and directly proportional to $\Delta t$. 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:

${\eta }_r^{L+\mathrm{1,}K}\approx S_r{\left(\Delta t\right)}^2$ (103)

${\eta }_v^{L+\mathrm{1,}K}\approx S_v\Delta t,$ (104)

where $S_r$ and $S_v$ are constants. The error resulted within an iteration is according to

${\left[S_r{\left(\Delta t_r\right)}^2\right]}^K{\varpi }_{i,r}^{L+\mathrm{1,0}}={\varpi }_{i,r}^{L+\mathrm{1,}K}$ (105)

${\left(S_v\Delta t_v\right)}^K{\varpi }_{i,v}^{L+\mathrm{1,0}}={\varpi }_{i,v}^{L+\mathrm{1,}K}.$ (106)

where ${\varpi }_{i,r}^{L+\mathrm{1,0}}=r_i^{L+\mathrm{1,}K=1}-r_i^L$ and ${\varpi }_{i,v}^{L+\mathrm{1,0}}=v_i^{L+\mathrm{1,}K=1}-v_i^L$. These above equations can be rewritten as

$\Delta t_r=\frac{1}{S_r{}^{\frac{1}{2}}}{\left(\frac{{\varpi }_{i,r}^{L+\mathrm{1,}K}}{{\varpi }_{i,r}^{L+\mathrm{1,0}}}\right)}^{\frac{1}{2K}},$ (107)

$\Delta t_v=\frac{1}{S_v}{\left(\frac{{\varpi }_{i,v}^{L+\mathrm{1,}K}}{{\varpi }_{i,v}^{L+\mathrm{1,0}}}\right)}^{\frac{1}{K}}.$ (108)

According to Zohdi (2014a), in the end of the time-step iteration, ${\varpi }_{i,r}^{L+\mathrm{1,}K}=TO{L}_r$ and ${\varpi }_{i,v}^{L+\mathrm{1,}K}=TO{L}_v$, where $TO{L}_r$ and $TO{L}_v$ are tolerances for the position and velocity vectors, respectively. So, assuming that $S_r$ and $S_v$ are constants, a new smaller step size can be defined as

${\mathrm{\Delta }t}_{tol}=\mathrm{\Delta }t\bullet min\left\{\underset{\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\mathrm{\Lambda }}_{K,r}}{\underbrace{\left[\frac{{\left(\frac{{TOL}_r}{{\mathit{\varpi }}_{i,r}^{L+\mathrm{1,0}}}\right)}^{\frac{1}{2K_d}}}{{\left(\frac{{\mathit{\varpi }}_{i,r}^{L+1,K}}{{\mathit{\varpi }}_{i,r}^{L+\mathrm{1,0}}}\right)}^{\frac{1}{2K}}}\right]}};\underset{\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\mathrm{\Lambda }}_{K,v}}{\underbrace{\left[\frac{{\left(\frac{{TOL}_v}{{\mathit{\varpi }}_{i,v}^{L+\mathrm{1,0}}}\right)}^{\frac{1}{2K_d}}}{{\left(\frac{{\mathit{\varpi }}_{i,v}^{L+1,K}}{{\mathit{\varpi }}_{i,v}^{L+\mathrm{1,0}}}\right)}^{\frac{1}{2K}}}\right]}}\right\}=\mathrm{\Delta }t\bullet \underset{{\mathrm{\Lambda }}_K}{\underbrace{min\left({\mathrm{\Lambda }}_{K,r};{\mathrm{\Lambda }}_{K,v}\right)}}=\mathrm{\Delta }t\bullet {\mathrm{\Lambda }}_K,$ (109)

where $K_d$ 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 $K_d$ 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 $v_p=500m/s$, 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.

Figure 10
Numerical model geometry: impact of a spherical projectile with an incoming velocity of $500m/s$ into a fabric shield (one sheet).

The geometry of the system is presented in Figure 10 and the parameters used in the numerical simulation are given as following.

• Numerical solution scheme:

• ¬ trapezoidal time-stepping parameter $\varphi =0.5$ (mid-point rule);

• ¬ initial time-step size of $\Delta t={10}^{-10}\sec onds$;

• ¬ iterative limit per time-step $K_d=20$;

• ¬ iterative tolerances per time-step $TO{L}_r={10}^{-6}$, $TO{L}_v={10}^{-6}$, $TO{L}_w={10}^{-6}$;

• ¬ gravity acceleration of $g=\left(0;0;-9.81\right)m/s$;

• ¬ damping coefficient $c_i=0.01N\cdot s/m$.

• Fabric shield:

• ¬ size of the sheet: $0.254m\times 0.254m$;

• ¬ lumped masses: $100\times 100$ (100 nodes in each direction);

• ¬ initial length $L^r=0.00254m$ (stretched, but without tension).

• Yarn properties:

• ¬ the yarn-segment is made by Zylon;

• ¬ Young’s modulus of a Zylon microfibril ${\mathbb{D}}^{yarn}=180GPa$;

• ¬ radius of a Zylon microfibril $R_{fibril}=3.9555\cdot {10}^{-5}m$;

• ¬ yarn with approximately $N_{fibril}=350$ microfibrils;

• ¬ radius of the yarn-segment $R_{yarn}=7.4\cdot {10}^{-4}m$;

• ¬ critical stretch to failure of a Zylon micro-fibril of $U_{crit}=1.034$;

• ¬ rupture around $U_{rupture}=1.06$;

• ¬ maximum Cauchy stress ${\sigma }_y=6435.66MPa$;

• ¬ density of a Zylon microfibril ${\rho }^{yarn}=1540kg/m^3$;

• ¬ damage rate for the yarn-segment $\lambda =135$;

• ¬ plastic rate for the yarn-segment $\mathcal{P}=0$.

• Coating properties:

• ¬ the coating-segment is made by a metal material;

• ¬ Young’s modulus of the coating ${\mathbb{D}}^{yarn}=50GPa$;

• ¬ thickness of the coating $\mathcal{T}=1.85\cdot {10}^{-4}m$ ($25\%$ of the radius of the yarn-segment);

• ¬ critical stretch to failure $U_{crit}=1.002$;

• ¬ rupture around $U_{rupture}=1.04$;

• ¬ maximum Cauchy stress ${\sigma }_y=100.30MPa$;

• ¬ density of the coating ${\rho }^{coating}=5000kg/m^3$;

• ¬ damage rate for the coating-segment $\lambda =135$;

• ¬ plastic rate for the coating-segment $\mathcal{P}=0$;

• Projectile:

• ¬ Young’s modulus of ${\mathbb{D}}^p=210GPa$;

• ¬ specific mass of ${\rho }^p=4522.55kg/m^3$;

• ¬ diameter of the spherical projectile $d_p=0.025m$;

• ¬ mass of the projectile $m_p=0.037kg$;

• ¬ dynamic viscosity of the air $c^e=0.00265N\cdot s/m^2$;

• ¬ local average velocity of the external interstitial medium of $v^e=0m/s$;

• ¬ projectile velocity of $v_p=500m/s$.

Figure 11
Motion configurations over time.

After the numerical simulation, some configurations were selected to illustrate the motion of the projectile and the behavior of the fabric shield, see Figure 11. Each sheet of fabric shield contains 10.000 nodes and 19.800 elements.

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\cdot {10}^{-4}s$ 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.

Figure 12
Damage of yarn and coating at time $3.2\cdot {10}^{-4}s$.

The stresses of the yarns and coatings at time $3.2\cdot {10}^{-4}s$ are presented in Figure 13.

Figure 13
Cauchy stress of yarn and coating at time $3.2\cdot {10}^{-4}s$.

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.

Figure 14
Particle and fabric shield configurations over time (with one layer).

The velocity of the projectile after crossing the fabric shield is $v_p=367.13m/s$. 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.

Figure 15
Front view of the impact of the projectile against a sheet of Zylon.

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).

Figure 16
Side view of the impact of the projectile against a sheet of Zylon.

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\cdot {10}^{-5}s$.

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.

Figure 18
Particle and fabric shield configurations over time (with twelve layers).

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 the velocity of the projectile over the time for the twelve models realized. The values of the velocities were plotted from the time $t=1.15\cdot {10}^{-4}s$ to $v_p=\mathrm{131,50}m/s$. It is noticed that with the addition of fabric shield layers a significant drop in speed occurs. From the ninth layer this speed drop becomes small.

Figure 19
Relationship between velocity and time for each model from 1 sheet to 12 sheets.

The Figure 20 shows the configuration of the system, composed by the projectile and six layers of fabric shield, at time $1.15\cdot {10}^{-4}s$. 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 $1.82\cdot {10}^{-4}s$, 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.

Figure 20
Structural fabric configuration at time $v$ for the model with 6 sheets.

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.

Figure 21
Structural fabric configuration at time $1.82\cdot {10}^{-4}s$ for the model with 9 sheets.

Another study was carried out observing the projectile in the configuration defined at time $7.7\cdot {10}^{-5}s$. 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 ($v$) and the number of layers of fabric shield ($x$) 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 $v=-\mathrm{0,0123}{x}^5+\mathrm{0,4375}{x}^4-\mathrm{5,7882}{x}^3+\mathrm{37,509}{x}^2-\mathrm{156,82}x+\mathrm{493,45}$. With this, it has been proven that there is a highly non-linear relationship when layers are added to the fabric shield.

Figure 22
Relationship between velocity and number of sheets at time $7.7\cdot {10}^{-5}s$ and the tendency curve.

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.

Figure 24
Front view of the impact of the projectile against multiple sheets of Zylon.

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.

ACKNOWLEDGEMENTS

The author wishes to express sincere appreciation to the independent public foundation Capes (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior, Ministry of Education of Brazil, Brasília-DF, Brazil), which mission is to foster research and the scientific and technological development. This work was supported by Capes, under the grants BEX 2565/15-3. The authors would like to thank the CRML (Computational Materials Research Laboratory, University of California at Berkeley, USA), specially Tarek Ismail Zohdi for the collaboration and support during my post-doctorate.

RFERENCES

• Zohdi, T. I.; Powell, D. (2006). Multiscale construction and large-scale simulation of structural fabric undergoing ballistic impact. Comput. Methods Appl. Mech. Engrg., v. 195 (1-3): 94-109.
• Zohdi, T. I. (Mar. 26, 2003). Genetic design of solids possessing a random-particulate microstructure. The Royal Society.
• Zohdi, T. I.; Wriggers, P. (2001b). Modeling and simulation of the decohesion of particulate aggregates in a binding matrix. Engineering Computations, vol. 18 : 79-95.
• Bandeira, A. A.; Zohdi, T. I. (2019). 3D numerical simulations of granular materials using DEM models considering rolling phenomena. Computational Particle Mechanics, vol. 6 : 97-131.
• Kato, S.; Yoshino, T.; Minami, H. (1999). Formulation of constitutive equations for fabric membranes based on the concept of fabric lattice model. Engineering Structures, vol. 21 : 691-708.
• Pargana, J. B. (2007). Advanced material model for coated fabrics used in tensioned fabric structures. Engineering Structures, vol. 29 : 1323-1336.
• Zohdi, T. I. (2014c). Impact and penetration resistance of network models of coated lightweight fabric shielding. GAMM-Mitt., v. 37, n. 1: 124 - 150.
• Kollegal, M. G.; Sridharan, S. (2000). Strength prediction of plain woven fabrics. J. Compos. Mater., v. 34: 240-257.
• Lanir, Y. (1983). Constitutive equations for fibrous connective tissues. J. Biomech, v. 16 (1): 1-12.
• Lim, C. T.; Shim, V. P.; Ng, Y. H. (2003). Finite element modeling of the ballistic impact of fabric armour. Int. J. Impact Engrg., v. 28, n. 1: 13-31.
• Powell, D.; Zohdi, T. I. (2009a). A note on flaw-induced integrity reduction of structural fabric. The International Journal of Fracture/Letters in Micromechanics, v. 158: L89-L96.
• Powell, D.; Zohdi, T. I. (2009b). Attachment mode performance of network-modeled ballistic fabric shielding. Composites Part B: Engineering, v. 40, Issue 6: 451-460.
• Tabiei, Y. J. (1999). Woven fabric composite material model with material nonlinearity for finite element simulation. Int. J. Solids Struct., v. 36: 2757-2771.
• Zohdi, T. I. (2002b). Modeling and simulation of progressive penetration of multilayered ballistic fabric shielding. Comput. Mech., v. 29: 61-67.
• Zohdi, T. I. (2007). A computational framework for network modeling of fibrous biological tissue deformation and rupture. Computer Methods in Applied Mechanics and Engineering, v. 196: 2972-2980.
• Zohdi, T. I. (2009). Microfibril-based estimates of the ballistic limit of multilayered fabric shielding. The International Journal of Fracture/Letters in Micromechanics, v. 158: L81-L88.
• Zohdi, T. I.; Steigmann, D. J. (2002). The toughening effect of microscopic filament misalignment on macroscopic fabric response. Int. J. Fract., v. 115: L9-L14.
• Zohdi, T. I.; Wriggers, P. (2001a). Computational micro-macro material testing. Arch. Comp. Meth. Eng., v. 9: 131-229.
• Zohdi, T. I.; Wriggers, P. (2008). Introduction to computational micromechanics. Second Reprinting. ed. Berlin: Springer,
• Atai, A. A.; Steigmann, D. J. (1997). On the nonlinear mechanics of discrete networks. Arch. Appl. Mech., v. 67: 303-319.
• Taylor, W.J., Jr.; Vinson, J. R. (1990). Modelling ballistic impact into flexible materials. AIAA J., vol. 28 : 2098-2103.
• Zohdi, T. I. (2010a). High-speed impact of electromagnetically sensitive fabric and induced projectile spin. Comput Mech, vol. 46 : 399-415.
• Wriggers, P. (2002). Computational Contact Mechanics. [S.l.]: John Wiley & Sons Ltd, England,
• Zohdi, T. I. (2010b). On the dynamics of charged electromagnetic particulate jets. Arch. Comput. Methods Eng., vol. 17 : 109-135.
• Zohdi, T. I. (2014a). A direct particle-based computational framework for electrically-enhanced thermo-mechanical sintering of powdered materials. Math. Mech. Solids : 1-21.
• Zohdi, T. I. (2002a). An adaptive-recursive staggering strategy for simulating multifield coupled processes in microheterogeneous solids. Int. J. Numer. Methods Engrg., v. 53: 1511-1532.
• Zohdi, T. I. (2014b). Additive particle deposition and selective laser processing-a computational manufacturing framework.Comput. Mech., vol. 54 : 171-191.
• Chen, W. F.; Han, D. J. (1988). Plasticity for Structural Engineers. New York: Springer-Verlag,
• Hill, R. (1950). The mathematical theory of plasticity. [S.l.]: Oxford University Press, Oxford (UK),
• Kachanov, L. M. (2004). Fundamentals of the Theory of Plasticity. [S.l.]: Dover Publications,
• Lubliner, J. (1984). A maximum-dissipation principle in generalized plasticity. Acta Mechanica, vol. 52 : 225-237.
• Pimenta, P. M. (1992). Finite Deformation Soil Plasticity on Principal Axes. III Complas. Barcelona, Espanha: [s.n.].
• Simo, J. C.; Ortiz, M. (1985). A unified approach to finite deformation elastoplasticity based on the use of hyperelastic constitutive equations. Comp. Meth. Appl. Mech. Engrg., vol. 49 : 221-245.
• Simo, J. C.; Taylor, R. L. (1985). Consistent Tangent Operators for Rate-Independant Elasto-plasticity. Computer Methods in Applied Mechanics and Engineering, vol. 48 : 101-18.
• Simo, J. C. (1992). Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Computer Methods in Applied Mechanics and Engineering, North-Holland, vol. 99 : 61- 112.
• Wriggers, P.; Eberlein, R.; Reese, S. (1996). A comparison of three-dimensional continuum and shell elements for finite plasticity. Int. J. Solids and Structures, vol. 33 : 3309-3326.

APPENDIX A: NUMERICAL SOLUTION SCHEME ALGORITHM

Algorithm 2
Numerical solution scheme.

ITERATIVE IMPLICIT SOLUTION METHOD
Step 1	Start a global fixed iteration and set the iteration counter $K=0$. Store the variables: $r_i^L\left(r_i^{L+\mathrm{1,0}}\right)$, $v_i^L\left(v_i^{L+\mathrm{1,0}}\right)$, $r_p^L\left(r_p^{L+\mathrm{1,0}}\right)$, $v_p^L\left(v_p^{L+\mathrm{1,0}}\right)$, ${\theta }_p^L\left({\theta }_p^{L+\mathrm{1,0}}\right)$ and ${\omega }_p^L\left({\omega }_p^{L+\mathrm{1,0}}\right)$. Compute other forces: $F_i^{f,oth,L}\left(r_i^L,v_i^L,r_p^L,v_p^L,{\omega }_p^L\right)$ and $F_p^{oth,L}\left(r_i^L,v_i^L,r_p^L,v_p^L,{\omega }_p^L\right)$. Compute the total force: $F_i^{tot,L}\left(r_i^L,v_i^L,r_p^L,v_p^L,{\omega }_p^L\right)$ and $F_p^{tot,L}\left(r_i^L,v_i^L,r_p^L,v_p^L,{\omega }_p^L\right)$. Compute the total momentum: $M_p^{tot,L}\left(r_i^L,v_i^L,r_p^L,v_p^L,{\omega }_p^L\right)$. Update the iteration counter $K=K+1$.
Step 2	Make a loop from $i=1$ to $i=N_f$ and compute for each fabric node $i$: a) Compute forces: $F_i^{f,pred,oth\mathrm{,,}L+\mathrm{1,}K-1}\left(r_i^{L+\mathrm{1,}K-1},v_i^{L+\mathrm{1,}K-1},r_p^{L+\mathrm{1,}K-1},v_p^{L+\mathrm{1,}K-1},{\omega }_p^{L+\mathrm{1,}K-1}\right)$. b) Compute mass predicted position and predicted velocity: $r_i^{f,pred,L+\mathrm{1,}K}=r_i^{f,L}+v_i^{f,L}\Delta t+\frac{\varphi {\left(\Delta t\right)}^2}{m_i}\left[\varphi F_i^{f,pred,oth,L+\mathrm{1,}K-1}+\left(1-\varphi \right)F_i^{f,oth,L}\right]$, $v_i^{f,pred,L+\mathrm{1,}K}=v_i^{f,L}+\frac{\Delta t}{m_i}\left[\varphi F_i^{f,pred,oth,L+\mathrm{1,}K-1}+\left(1-\varphi \right)F_i^{f,oth,L}\right]$. c) Enforce contact if necessary (if contact occurs do the following procedures): Consider “stick case”. Make $r_i^{f,corr,L+\mathrm{1,}K}=r_i^{p,cp,L+\mathrm{1,}K}$ and $v_i^{f,corr,L+\mathrm{1,}K}=v_i^{p,cp,L+\mathrm{1,}K}$. Compute $F_i^{con}=\frac{m_i}{\Delta t^2}\left[r_i^{f,corr,L+\mathrm{1,}K}-r_i^{f,pred,L+\mathrm{1,}K}\right]$ and $F_i^{fric}=\frac{m_i}{\Delta t}\left(v_i^{f,corr,L+\mathrm{1,}K}-v_i^{f,L}\right)-F_i^{con}-F_i^{f,oth,L}$. IF $F_i^{fric}\le {\mu }_s{F}_i^{con}$ THEN “stick case”. Leave like that. The contact condition is satisfied. ELSE “slip case”. Update the velocity $v_i^{f,corr,L+\mathrm{1,}K}=v_i^{f,pred,L+\mathrm{1,}K}+\frac{\Delta t}{m_i}\left[F_i^{con}+F_i^{fric}\right]$ and the friction force $F_i^{fric}={\mu }_d{F}_i^{con}\frac{v_i^{p,cp,L+\mathrm{1,}K}-v_i^{f,corr,L+\mathrm{1,}K}}{v_i^{p,cp,L+\mathrm{1,}K}-v_i^{f,corr,L+\mathrm{1,}K}}$.
Step 3	Make a loop from $p=1$ to $p=N_p$ and compute for each projectile $p$: a) Compute forces: $F_p^{oth,L+\mathrm{1,}K-1}\left(r_i^{L+\mathrm{1,}K},v_i^{L+\mathrm{1,}K},r_p^{L+\mathrm{1,}K-1},v_p^{L+\mathrm{1,}K-1},{\omega }_p^{L+\mathrm{1,}K-1}\right)$. b) Compute momentum: $M_p^{tot,L+\mathrm{1,}K-1}\left(r_i^{L+\mathrm{1,}K},v_i^{L+\mathrm{1,}K},r_p^{L+\mathrm{1,}K-1},v_p^{L+\mathrm{1,}K-1},{\omega }_p^{L+\mathrm{1,}K-1}\right)$. c) Compute position, velocity, angular velocity and rotation: $r_p^{L+\mathrm{1,}K}=r_p^L+v_p^L\Delta t+\frac{\varphi {\left(\Delta t\right)}^2}{m_p}\left\{-{\displaystyle \sum }_{i=1}^{N_c}\left(F_i^{con}+F_i^{fric}\right)+\left[\varphi F_p^{oth,L+\mathrm{1,}K-1}+\left(1-\varphi \right)F_p^{oth,L}\right]\right\}$, $v_p^{L+\mathrm{1,}K}=v_p^L-\frac{\Delta t}{m_p}{\displaystyle \sum }_{i=1}^{N_c}\left(F_i^{con}+F_i^{fric}\right)+\frac{\Delta t}{m_p}\left[\varphi F_p^{oth,L+\mathrm{1,}K-1}+\left(1-\varphi \right)F_p^{oth,L}\right]$. ${\omega }_p^{L+\mathrm{1,}K}={\omega }_p^L+\frac{\Delta t}{{\overline{I}}_{i,s}}\left[\varnothing M_p^{tot,L+\mathrm{1,}K-1}+\left(1-\varnothing \right)M_p^{tot,L}\right]$. ${\theta }_p^{L+\mathrm{1,}K}={\theta }_p^L+\Delta t\left[\varnothing {\omega }_p^{L+\mathrm{1,}K-1}+\left(1-\varnothing \right){\omega }_p^L\right]$.
Step 4	Compute/update forces: $F_i^{tot,L+\mathrm{1,}K}$ (similar to step 2).
Step 5	Compute/update fabric/coating damage: ${\alpha }_I^{yarn}$ and ${\alpha }_I^{coating}$.
Step 6	Compute/update fabric/coating plastic stretch: $U_{p,I}^{yarn}$ and $U_{p,I}^{coating}$.
Step 7	Measure normalized error quantities: a) ${\varpi }_{K,r}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{\sum _{i=1}^{N_f+N_p}‖{\mathit{r}}_i^{L+1,K}-{\mathit{r}}_i^{L+1,K-1}‖}{\sum _{i=1}^{N_f+N_p}‖{\mathit{r}}_i^{L+1,K}-{\mathit{r}}_i^L‖}$ , ${\varpi }_{K,v}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{\sum _{i=1}^{N_f+N_p}‖{\mathit{v}}_i^{L+1,K}-{\mathit{v}}_i^{L+1,K-1}‖}{\sum _{i=1}^{N_f+N_p}‖{\mathit{v}}_i^{L+1,K}-{\mathit{v}}_i^L‖}$ and ${\varpi }_{K,\omega }\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{\sum _{p=1}^{N_p}‖{\mathit{\omega }}_p^{L+1,K}-{\mathit{\omega }}_p^{L+1,K-1}‖}{\sum _{p=1}^{N_p}‖{\mathit{\omega }}_p^{L+1,K}-{\mathit{\omega }}_p^L‖}$. b) $Z_{k,r}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{{\varpi }_{K,r}}{{TOL}_r}$ , .$Z_{k,v}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{{\varpi }_{K,v}}{{TOL}_v}$ and $Z_{k,\omega }\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{{\varpi }_{K,\omega }}{{TOL}_{\omega }}$ . c) ${\mathrm{\Lambda }}_{K,r}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{{\left(\frac{{TOL}_r}{{\varpi }_{0,r}}\right)}^{\frac{1}{2K_d}}}{{\left(\frac{{\varpi }_{K,r}}{{\varpi }_{0,r}}\right)}^{\frac{1}{2K}}}$ , ${\mathrm{\Lambda }}_{K,v}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{{\left(\frac{{TOL}_v}{{\varpi }_{0,v}}\right)}^{\frac{1}{2K_d}}}{{\left(\frac{{\varpi }_{K,v}}{{\varpi }_{0,v}}\right)}^{\frac{1}{2K}}}$ and ${\mathrm{\Lambda }}_{K,\omega }\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{{\left(\frac{{TOL}_{\omega }}{{\varpi }_{0,\omega }}\right)}^{\frac{1}{2K_d}}}{{\left(\frac{{\varpi }_{K,\omega }}{{\varpi }_{0,\omega }}\right)}^{\frac{1}{2K}}}$.
Step 8	IF $\left(Z_{k,r}\le 1AND{Z}_{k,v}\le 1AND{Z}_{k,\omega }\le 1\right)ANDK\le K_d$ THEN Increment time: $t=t+\Delta t$, Construct the next time step: ${\left(\Delta t\right)}^{new}={\Lambda }_K{\left(\Delta t\right)}^{old}$, where ${\Lambda }_K=\min \left[{\Lambda }_{K,r},{\Lambda }_{K,v},{\Lambda }_{K,\omega }\right]$, Select the minimum size: $\Delta t=\min \left[{\left(\Delta t\right)}^{\lim },{\left(\Delta t\right)}^{new}\right]$, Go to step 1. ELSE Go to step 9.
Step 9	IF $\left(Z_{k,r}>1OR{Z}_{k,v}>1OR{Z}_{k,\omega }>1\right)ANDK<K_d$ THEN Update the iteration counter: $K=K+1$, Go to step 2. ELSE Go to step 10.
Step 10	IF $\left(Z_{k,r}>1OR{Z}_{k,v}>1OR{Z}_{k,\omega }>1\right)ANDK=K_d$ THEN Construct a next time step: ${\left(\Delta t\right)}^{new}={\Lambda }_K{\left(\Delta t\right)}^{old}$, where ${\Lambda }_K=\min \left[{\Lambda }_{K,r},{\Lambda }_{K,v},{\Lambda }_{K,\omega }\right]$, Select the minimum size: $\Delta t=\min \left[{\left(\Delta t\right)}^{\lim },{\left(\Delta t\right)}^{new}\right]$, Restart at time $t$: Set $r_i^{L+\mathrm{1,0}}=r_i^L$, $v_i^{L+\mathrm{1,0}}=v_i^L$, $r_p^{L+\mathrm{1,0}}=r_p^L$, $v_p^{L+\mathrm{1,0}}=v_p^L$ and ${\omega }_p^{L+\mathrm{1,0}}={\omega }_p^L$, Go to step 1. ELSE Set a larger number for the variable $K_d$, Restart at time $t$: Set $r_i^{L+\mathrm{1,0}}=r_i^L$, $v_i^{L+\mathrm{1,0}}=v_i^L$, $r_p^{L+\mathrm{1,0}}=r_p^L$, $v_p^{L+\mathrm{1,0}}=v_p^L$ and ${\omega }_p^{L+\mathrm{1,0}}={\omega }_p^L$, Go to step 1.

APPENDIX B: OVERHAUL OF SOLID MECHANICS

From the continuum mechanics, the second Piola-Kirchhoff stress tensor $S$ is defined by

$S=F^{-1}P,$ (110)

where $P$ is the first Piola-Kirchhoff stress tensor that correlates nominal stress per unit reference area and $F$ is the deformation gradient defined as

$F=\nabla x=\frac{\partial x}{\partial x^r}.$ (111)

In above equation, $x$ and $x^r$ are the position vector of the nodes in current configuration and reference configuration, respectively. The Cauchy stress tensor $\sigma$ associates the superficial force vector $t$ (or stress vector per unit of area in the current configuration) with the unit normal vector at the surface in current configuration $n$ defined by

$t=\sigma n.$ (112)

Using the Nanson theory and considering

$J=\det F$ (113)

and

$t^r=P{n}^r,$ (114)

is possible to demonstrate the following relation

$P=J\sigma F^{-T}.$ (115)

The Kirchhoff-Trefftz stress tensor is defined by

$\Sigma =P{F}^T=FS{F}^T=J\sigma .$ (116)

Note that the second Piola-Kirchhoff stress tensor $S$ is symmetric and doesn’t have a physical interpretation, the first Piola-Kirchhoff stress $P$ is nonsymmetrical, the Kirchhoff-Trefftz stress tensor $\Sigma$ is symmetrical, the Cauchy stress tensor $\sigma$ is symmetrical and they are energetically conjugate by this relation $S:\dot{E}=P:\dot{F}$. The displacements gradient $L$ is defined by

$L=\nabla u=F-I,$ (117)

where $u$ is the displacements vector of the nodes and $I$ is the identity matrix. The quadratic stretches tensor or right Cauchy-Green strain $C$ is defined by

$C=F^{T}F.$ (118)

The Green-Lagrange strain tensor is calculated by

$E=\frac{1}{2}\left(C-I\right).$ (119)

The stretch tensor $U$ is calculated by

$U^2=C.$ (120)

The Green-Lagrange strain tensor can be written also as

$E=\frac{1}{2}\left(U^2-I\right).$ (121)

The following strains belong to the Hill family:

${\epsilon }_m=\{\begin{array}{c}\frac{1}{m}\left({\lambda }^m-1\right),ifm\ne 0\\ \ln \left(\lambda \right),ifm=0\end{array},$ (122)

where $m$ is a integer value and $\lambda$ is the stretch. The equation (121) corresponds to (122) when the parameter $m$ is equal to 2.

The rotation tensor $R$ can be calculated by the following polar decompositions

$F=RU=VR,$ (123)

where $V$ is the left stretch tensor. The rotation tensor is orthogonal, i.e., $\det R=+1$.

For hyperelastic materials, the following relations that correlates stresses with strains are valid:

$P=\frac{\partial \psi }{\partial F},$ (124)

$S=\frac{\partial \psi }{\partial E},$ (125)

where $\psi \left(F\right)$ and $\psi \left(E\right)$ 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

$\mathbb{D}=\frac{{\partial }^2\psi }{\partial F^2}=\frac{\partial P}{\partial F},$ (126)

$\mathbb{D}=\frac{{\partial }^2\psi }{\partial E^2}=\frac{\partial S}{\partial E}.$ (127)

APPENDIX C: OVERHAUL OF ELASTIC-PLASTICITY MATERIAL

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. The Figure 25b represents the elastic-plastic model with hardening.

Figure 25
Uniaxial test of elastic plastic material: a) elastic-perfectly-plastic model, b) elastic-plastic model with hardening.

When the initial limit tension ${\sigma }_{y0}$ is reached, the limit tension value ${\sigma }_y$ 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

$F=F^e{F}^n,$ (128)

where $F$ is total deformation gradient, $F^e$ is the elastic part of the deformation gradient from the natural to the current configuration and $F^n=F^p$ corresponds to the plastic part from the initial reference configuration to the natural one, see Figure 26.

Figure 26
Geometrical interpretation of the tensors for elastic-plastic constitutive equations.

In most cases, one sets $F^p\left(0\right)=I$. The elastic and the plastic parts have the following polar decompositions

$F^e=R^e{U}^e=V^e{R}^e,$ (129)

$F^p=R^p{U}^p=V^p{R}^p,$ (130)

where $U^e$ is the right elastic stretching tensor, $V^e$ is the left elastic stretching tensor, $R^e$ is the elastic rotation tensor, $U^p$ is the right plastic stretching tensor, $V^p$ is the left plastic stretching tensor and $R^p$ is the plastic rotation tensor.

As the natural configuration is not fixed, this means that ${\dot{F}}^p\ne 0$ 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 $F_t$ and given by

$F_t=F{F}^n{}^{-1},$ (131)

see again Figure 26 for a geometrical representation. This gradient has the following polar decomposition

$F_t=R_t{U}_t.$ (132)

At each instant, the following identities hold

$F_t=F^e,$ (133)

$R_t=R^e,$ (134)

$U_t=U^e.$ (135)

Considering these relations $\frac{\partial v}{\partial x}=\frac{\partial v}{\partial x^r}\frac{\partial x^r}{\partial x}=\dot{F}{F}^{-1}$ and $\frac{\partial v}{\partial x^n}=\frac{\partial v}{\partial x^r}\frac{\partial x^r}{\partial x^n}=\dot{F}{F}^n{}^{-1}$, where $\frac{\partial v}{\partial x^r}=\dot{F}$ and $\frac{\partial v}{\partial x^n}={\dot{F}}_t$, the velocity gradient from the instantaneously fixed natural configuration is

${\dot{F}}_t=\dot{F}{F}^n{}^{-1}.$ (136)

Changing $F^p$ and ${\dot{F}}^p$ to the same reference configuration, one concludes that

$F_t^p=R_t^p=U_t^p=I,$ (137)

${\dot{F}}_t^p={\dot{F}}^p{F}^p{}^{-1}.$ (138)

Hence, one has

${\dot{F}}_t={\dot{F}}^e+F^e{\dot{F}}_t^p{F}^e{}^{-1}.$ (139)

Since the equation (137) is valid, one has

$E_t^p=U_t^p.$ (140)

For the same reason one may instantaneously write

${\dot{E}}_t^p={\dot{U}}_t^p.$ (141)

The elastic Green-Lagrange strain tensor, as shown in equation (121), is defined by

$E^e=\frac{1}{2}\left(U^e{}^2-I\right).$ (142)

Considering $c_i^e$ the principal axes of the right Cauchy-Green strain tensors $C$, once $U^e$ has the spectral decomposition $U^e={\lambda }_i^e{c}_i^e\otimes c_i^e$, where ${\lambda }_i^e=\sqrt{c_i^e}$, and using the relations present in (121) one may write $E^e$ as the spectral decomposition below

$E^e={\epsilon }_m^e{c}_i^e\otimes c_i^e,$ (143)

where

${\epsilon }_m^e\left({\lambda }_i^e\right)=\{\begin{array}{c}\frac{1}{m}\left({\lambda }_i^e{}^m-1\right),ifm\ne 0\\ \ln \left({\lambda }_i^e\right),ifm=0\end{array}.$ (144)

The total strain from the fixed natural configuration is defined by

$E_t=E_t\left(U_t\right).$ (145)

In view of (133), at current time, the following identity holds

$E_t=E^e.$ (146)

Isotropic plasticity

In plasticity the elastic domain encompasses all admissible stress states and is described by the yield function

$F\left(S^e,\alpha \right)\le \mathrm{0,}$ (147)

where $F:\mathcal{S}\times {ℝ}^{n_h}\to ℝ$ and $\alpha$ is the $\left(n_h\times 1\right)$ vector of the hardening parameters. $F$ does not depend on hardening parameters in the particular case of ideal plasticity. The surface in $\mathcal{S}$ described by $F\left(S^e,\alpha \right)=0$ is called yield surface.

In isotropic plasticity, first is assumed that elasticity is isotropic and described by

$S^e={\mathbb{D}}^e{E}^e,$ (148)

${\dot{S}}^e={\mathbb{D}}^e{\dot{E}}^e,$ (149)

and the yield function (147) is written as displayed below

$F\left(S^e,\alpha \right)=\overline{\sigma }\left(S^e\right)-{\sigma }_y\left(\alpha \right)\le \mathrm{0,}$ (150)

in which $\overline{\sigma }$ is an isotropic function of its argument and equal to the uniaxial stress in a simple tension test. This means the following

$S^e=\sigma e\otimes e\Rightarrow \overline{\sigma }=\sigma ,$ (151)

where $e$ is the unitary vector in the axial direction of the test and $\sigma$ is the corresponding uniaxial stress. $\overline{\sigma }$ is then called equivalent stress. On the other hand ${\sigma }_y$ is the yield stress measured in a simple tension test and $\alpha$ is the isotropic hardening parameter. The plastic dissipation per unit volume in the natural configuration is given by the internal dissipation, as displayed below

${\mathcal{D}}^n={\mathcal{P}}_{\mathrm{int}}^n-\dot{\psi }=S^e:\left({\dot{E}}_t-{\dot{E}}^e\right),$ (152)

where $\psi \left(E^e\right)$ is the strain energy function per unit volume in the natural configuration of an elastic material and the operator is defined as $A:B=tr\left(A^{T}B\right)$, where $tr(\cdot )$ indicates the trace operator. In order to formulate the constitutive rate equations of plasticity, is stated that the Principle of Maximum Plastic Dissipation

$\underset{S^e\in \mathcal{S}}{\max }\left({\mathcal{D}}^n\right),subjecttoF\left(S^e\right)\le \mathrm{0,}$ (153)

is a necessary and sufficient condition for the determination of $S^e$. Mathematical Programming asserts that (153) can be transformed to the following unconstrained problem

$\frac{\partial }{\partial S^e}\left[-S^e:\left({\dot{E}}_t-{\dot{E}}^e\right)+\lambda F\right]=\mathrm{0,}$ (154)

where $\lambda$ is the Lagrange multiplier. Thus, one has

${\dot{E}}_t-{\dot{E}}^e=\lambda N,$ (155)

where

$N=\frac{\partial F}{\partial S^e}=\frac{\partial \overline{\sigma }}{\partial S^e}.$ (156)

$N$ is the normal to the yield surface $F=0$ in the stress space. Since $F$ is an isotropic function of $S^e$ and the elasticity is isotropic, $N$, $S^e$ and $U^e$ are coaxial. Thus, (155) implies that ${\dot{E}}_t-{\dot{E}}^e$, $N$, $S^e$ and $U^e$ are coaxial. It is easy to demonstrate the following result

${\dot{E}}_t={\dot{E}}^e+{\dot{E}}_t^p,$ (157)

which is similar to the additive decomposition of the total strain in small strain plasticity. From

${\dot{S}}_t^p=-{\mathbb{D}}^e\left({\dot{E}}_t-{\dot{E}}^e\right)$ (158)

and

${\dot{S}}^e={\mathbb{D}}^e\left({\dot{E}}_t-\left({\dot{E}}_t-{\dot{E}}^e\right)\right),$ (159)

one arrives at

${\dot{S}}_t^p=-{\mathbb{D}}^e{\dot{E}}_t^p$ (160)

and

${\dot{S}}^e={\mathbb{D}}^e\left({\dot{E}}_t-{\dot{E}}_t^p\right).$ (161)

From (157), (160) and (161), one necessarily has

${\dot{E}}_t^p=\lambda N.$ (162)

Mathematical Programming asserts also that $F$ must be convex and the following Kuhn-Tucker optimality conditions must hold $\lambda \ge 0$, $F\le 0$ and $\lambda F=0$. So, from (162) follows

$\Vert {\dot{E}}_t^p\Vert =\lambda \Vert N\Vert .$ (163)

When $\lambda =0$, i.e., when $F<0$, or when $F=0$ and $\dot{F}<0$, there is no plastic strain rate and ${\dot{E}}_t^p=0$. Hence, one has ${\dot{E}}_t={\dot{E}}^e$. From ${\dot{S}}^e={\mathbb{D}}^e{\dot{E}}^e$ follows ${\dot{S}}_t={\mathbb{D}}^e{\dot{E}}_t$. When $\lambda >0$, i.e., when $F=0$ and $\dot{F}=0$, one has from (150), by differentiation, the following equation

$\dot{F}=N:{\dot{S}}^e-h\dot{\alpha }=\mathrm{0,}$ (164)

where

$h=\frac{d{\sigma }_y}{d\alpha }.$ (165)

Defined $\alpha$ such that

${\mathcal{D}}^n=\left(N:S^e\right)\dot{\alpha },$ (166)

then, from (155) follows

$\lambda =\dot{\alpha }.$ (167)

Introducing (161) in (164), one has

$\dot{F}=N:{\mathbb{D}}^e\left({\dot{E}}_t-{\dot{E}}_t^p\right)-h\dot{\alpha }=0.$ (168)

With the aid of (162) the equation (168) can be rewritten by

$\dot{F}=N:{\mathbb{D}}^e{\dot{E}}_t-N:{\mathbb{D}}^e\lambda N-h\dot{\alpha }=0.$ (169)

Introducing (167) into (169), obtain

$\dot{F}=N:{\mathbb{D}}^e{\dot{E}}_t-\left(N:{\mathbb{D}}^{e}N+h\right)\dot{\alpha }=0.$ (170)

The above equation leads to the following result

$\dot{\alpha }=\frac{1}{N:{\mathbb{D}}^{e}N+h}N:{\mathbb{D}}^e{\dot{E}}_t.$ (171)

Note that for $\dot{\alpha }>0$ and introducing the equations (162) and (167) into (161) can be rewritten by

${\dot{S}}_t={\mathbb{D}}^e\left({\dot{E}}_t-\dot{\alpha }N\right).$ (172)

Thus, one may write

${\dot{S}}_t=\{\begin{array}{cc}{\mathbb{D}}^e{\dot{E}}_t,& ifF<0orifF=0and\dot{F}<\mathrm{0,}\\ {\mathbb{D}}^{ep}{\dot{E}}_t,& ifF=0and\dot{F}=\mathrm{0,}\end{array}$ (173)

where

${\mathbb{D}}^{ep}={\mathbb{D}}^e-\frac{1}{N:{\mathbb{D}}^{e}N+h}\left({\mathbb{D}}^{e}N\right)\otimes \left({\mathbb{D}}^{e}N\right).$ (174)

${\mathbb{D}}^{ep}$ is called the tensor of elastic-plastic tangent moduli for the pair $\left\{S_t,E_t\right\}$. If one rephrases (153) as

$\underset{S^e\in \mathcal{S}}{\min }\left(\frac{\partial {\psi }^{*}}{\partial S^e}:{\dot{S}}_t^p\right),subjecttoF\left(S^e\right)\le 0.$ (175)

From (175) results the following unconstrained minimization problem

$\underset{S^e\in \mathcal{S}}{\min }\left(\frac{\partial {\psi }^{*}}{\partial S^e}:{\dot{S}}_t^p+\lambda F\left(S^e,\alpha \right)\right).$ (176)

Mathematical Programming then asserts that

$\frac{\partial }{\partial S^e}\left(\frac{\partial {\psi }^{*}}{\partial S^e}:{\dot{S}}_t^p+\lambda F\right)=\mathrm{0,}$ (177)

which leads to

$\frac{{\partial }^2{\psi }^{*}}{\partial S^e{}^2}:{\dot{S}}_t^p+\lambda N={\mathbb{D}}^e{}^{-1}{\dot{S}}_t^p+\lambda N=0.$ (178)

Introducing ${\dot{S}}_t^p=-{\mathbb{D}}^e\left({\dot{E}}_t-{\dot{E}}^e\right)$ 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

$J^p=\det U^p=\mathrm{1,}$ (179)

which implies

$Sph\left(\ln U^p\right)=Sph\left(\ln {\dot{U}}^p\right)=Sph\left({\dot{E}}_t^p\right)=0$ (180)

and

${\Sigma }_i={\Sigma }_i^e.$ (181)

According to (155), if $N$ is deviatoric, then ${\dot{E}}_t^p$ is deviatoric too. This implies (179). In the case of isochoric plasticity, one has ${\mathcal{P}}_{\mathrm{int}}^r={\mathcal{P}}_{\mathrm{int}}^n$ and ${\mathcal{D}}^r={\mathcal{D}}^n$.

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 $\left\{t_0,t_1,t_2,\cdots ,t_i,t_{i+1},\cdots \right\}$.

Consider the material state at instant $t_i$, at which $F_i$, $U_i^p$ and ${\alpha }_i$ are supposed to be known. Hence, from (128) one has

$F_i^e=F_i{U}_i^p{}^{-1}.$ (182)

If the polar decomposition

$F_i^e=R_i^e{U}_i^e$ (183)

is performed, the stress tensor at this instant can be computed through

$S_i^e=\frac{\partial \psi }{\partial E^e}\left(E_i^e\right),$ (184)

where

$E_i^e=\frac{1}{2}\left(U_i^e{}^2-I\right).$ (185)

The above equation lies to

${\dot{E}}_i^e=U_i^e{\dot{U}}_i^e.$ (186)

At instant $t_{i+1}$ only the deformation gradient $F_{i+1}$ is assumed to be known. The following trial elastic state is defined, assuming that the plastic part of the deformation has not changed,

$F_t^e=F_{i+1}{U}_i^p{}^{-1}.$ (187)

The polar decomposition

$F_t^e=R_t^e{U}_t^e,$ (188)

is then performed and the strain tensor $E_t^e$ is computed with the aid of

$E_t^e=\frac{1}{2}\left(U_t^e{}^2-I\right).$ (189)

The stress tensor in this state is given by

$S_t=\frac{\partial \psi }{\partial E^e}\left(E_t^e\right).$ (190)

If $F_t=F\left(S_t\right)<0$, then the trial state is admissible, $\Delta \alpha =0$, and

$S_{i+1}^e=S_t,{\alpha }_{i+1}={\alpha }_{i}and{U}_{i+1}^p=U_i^p.$ (191)

If $F_t\ge 0$, there is an incremental plastic stretch $\Delta U^p$ such that

$F_{i+1}^e=F_t^e\Delta U^p{}^{-1}and{F}_{i+1}^p=\Delta U^p{U}_i^p.$ (192)

The plastic stretch tensor at the end of the increment, denoted by $U_{i+1}^p$, is obtained from $F_{i+1}^p$ through

$U_{i+1}^p={\left(F_{i+1}^p{}^T{F}_{i+1}^p\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}.$ (193)

The incremental stress integration algorithm is formulated by the following constrained minimization problem, which must be solved for $S_{i+1}^e$,

$\min {\psi }^{*}\left(\Delta S^p\right)subjecttoF\left(S_{i+1}^e,{\alpha }_{i+1}\right)=\mathrm{0,}$ (194)

where

$\Delta S^p=S_{i+1}^e-S_t.$ (195)

The equation (195) is based on (175) and is equivalent to the following unconstrained minimization problem

$\min \left[{\psi }^{*}\left(\Delta S^p\right)+\Delta \alpha F\left(S_{i+1}^e\right)\right],$ (196)

where $\Delta \alpha$ is the Lagrange multiplier. (196) leads to the following system

$\{\begin{array}{c}\frac{\partial {\psi }^{*}}{\partial S^e}\left(S_{i+1}^e-S_t\right)+\Delta \alpha N_{i+1}=\mathrm{0,}\\ F\left(S_{i+1}^e,{\alpha }_i+\Delta \alpha \right)=0.\end{array}$ (197)

In (197), the following normal has been introduced

$N_{i+1}=\frac{\partial F}{\partial S^e}\left(S_{i+1}^e\right).$ (198)

The solutions of (197) are $S_{i+1}^e$ and $\Delta \alpha$. If ${\psi }^{*}$ and $F$ are convex, (194) has a unique solution and the Newton Method applied to (197) is globally convergent. From (194) results the following Kuhn-Tucker optimality conditions

$\Delta \alpha \ge \mathrm{0,}F\le 0and\Delta \alpha F=0.$ (199)

When ${\psi }^{*}$ is given by

$\psi \left(E^e\right)=\frac{1}{2}{E}^e:{\mathbb{D}}^e{E}^{e}and{\psi }^{*}\left(S^e\right)=\frac{1}{2}{S}^e:{\mathbb{D}}^e{}^{-1}{S}^e,$ (200)

delivers

$\{\begin{array}{c}{\mathbb{D}}^{e-1}\left(S_{i+1}^e-S_t\right)+\Delta \alpha N_{i+1}=\mathrm{0,}\\ F\left(S_{i+1}^e,{\alpha }_i+\Delta \alpha \right)=0.\end{array}$ (201)

Equation (201) is the incremental counterpart of (178). Hence, one may write

$\Delta E_t^p=-{\mathbb{D}}^e{}^{-1}\left(S_{i+1}^e-S_t\right)=\Delta \alpha N_{i+1},$ (202)

where (202) is the incremental counterpart of (178). After solving (197), one has $S_{i+1}^e$ and $\Delta \alpha$. Afterwards, the following quantities can be computed

$\begin{array}{l}\Delta E_t^p=\Delta \alpha N_{i+1},\\ \Delta U^p=U_i^p{}^{-1}\Delta E_t^p,\\ F_{i+1}^p=\Delta U^p{U}_i^p,\\ U_{i+1}^p={\left(F_{i+1}^p{}^T{F}_{i+1}^p\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.},\\ {\alpha }_{i+1}={\alpha }_i+\text{Δ}\alpha .\end{array}$ (203)

Note that, when $N_{i+1}$ is deviatoric, $\Delta E_t^p$ is also deviatoric and $\Delta U^p$ is isochoric. It is important to remark additionally that $N_{i+1}$, $S_{i+1}^e$, $S_t$, $\Delta U^p$ and $\Delta E_t^p$ are coaxial. Equation (201) yields also

$S_{i+1}^e=S_t-{\mathbb{D}}^e\Delta E_t^p=S_t-\Delta \alpha {\mathbb{D}}^e{N}_{i+1}.$ (204)

The analytical solution of (204) is given by

$S_{i+1}^e=Sph\left(S_t^e\right)+\frac{{\sigma }_{y_{i+1}}}{{\overline{\sigma }}_t}Dev\left(S_t^e\right),$ (205)

where

${\sigma }_{y_{i+1}}={\sigma }_y\left({\alpha }_{i+1}\right)={\sigma }_y\left({\alpha }_i+\Delta \alpha \right)and{F}_t={\overline{\sigma }}_t-{\sigma }_y\left({\alpha }_i\right).$ (206)

This solution is known as the radial return method. Thus, from (204), one has

$S_t-\Delta \alpha {\mathbb{D}}^e{N}_{i+1}=Sph\left(S_t^e\right)+\frac{{\sigma }_{y_{i+1}}}{{\overline{\sigma }}_t}Dev\left(S_t^e\right).$ (207)

Considering that

$S_t=Sph\left(S_t^e\right)+Dev\left(S_t^e\right),$ (208)

equation (207) can be rewritten as

$Sph\left(S_t^e\right)+Dev\left(S_t^e\right)-\Delta \alpha {\mathbb{D}}^e{N}_{i+1}=Sph\left(S_t^e\right)+\frac{{\sigma }_{y_{i+1}}}{{\overline{\sigma }}_t}Dev\left(S_t^e\right).$ (209)

The value of $\Delta \alpha$ is obtained from the above equation, just solving the following relation

$\Delta \alpha {\mathbb{D}}^e{N}_{i+1}=\left(1-\frac{{\sigma }_{y_{i+1}}}{{\overline{\sigma }}_t}\right)Dev\left(S_t^e\right)=\left(\frac{F_t}{{\overline{\sigma }}_t}\right)Dev\left(S_t^e\right).$ (210)

The so-called tensor of algorithmic tangent moduli is defined by

${\mathbb{D}}^{a\lg }=\frac{\partial S_{i+1}^e}{\partial E_t}.$ (211)

Remember that (211) is referenced to the natural configuration at the instant $t_i$ 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.

Figure 27
Interpretation of the update procedure for elastic-plastic constitutive equations.

For computational purposes the radial return algorithm on principal axes for the von Mises plasticity with isotropic hardening is summarized below.

Algorithm 3
Plasticity algorithmic.

I) Trial step:
a) $F_t^e=F_{i+1}{U}_i^p{}^{-1}$, $C_t^e={\left(F_t^e\right)}^T{F}_t^e$; b) $U_t^e={\left(C_t^e\right)}^{1/2}$, $E_t=\frac{1}{2}\left(U_t^e{}^2-I\right)$; c) $R_t^e=F_t^e{U}_t^{e-1}$, $S_t^e={\mathbb{D}}^e{E}_t$; d) $F_t={\overline{\sigma }}_t-{\sigma }_y\left({\alpha }_i\right)$
II) Radial return algorithm:
If $\left(F_t<0\right)$ then Elastic step: 1. ${\alpha }_{i+1}={\alpha }_i$; 2. $U_{i+1}^p=U_i^p$ ; 3. $S_{i+1}=S_t^e$; 4. ${\mathbb{D}}_{i+1}={\mathbb{D}}^e$; Else if $\left(F_t\ge 0\right)$ then Elastic-plastic step: 1. $\Delta \alpha$, ${\alpha }_{i+1}={\alpha }_i+\Delta \alpha$; ${\sigma }_y\left({\alpha }_{i+1}\right)={\sigma }_y\left({\alpha }_i\right)+h\Delta \alpha$; 2. $\Delta E^p=\Delta \alpha N_{i+1}$, $\Delta U^p=U_i^p{}^{-1}\Delta E^p$, $F_{i+1}^p=\Delta U^p{U}_i^p$, $U_{i+1}^p={\left[{\left(F_{i+1}^p\right)}^T{F}_{i+1}^p\right]}^{\frac{1}{2}}$ 3. $S_{i+1}=S_t^e-{\mathbb{D}}^e\Delta E^p$; 4. ${\mathbb{D}}_{i+1}={\mathbb{D}}^{ep}$; End if.

Edited by

Editor:

Marcílio Alves.

Publication Dates

• Publication in this collection
  24 June 2020
• Date of issue
  2020

History

• Received
  20 Feb 2020
• Reviewed
  25 Apr 2020
• Accepted
  28 Apr 2020
• Published
  04 May 2020