Three different methods were explored to model the skull fracture. First, the element is removed from the calculation when the failure criterion is satisfied. Element removal causes a sudden loss of mass, momentum, and energy. Therefore, numerical convergence could be slow and more expensive computation necessary to maintain desired accuracy. However, this method is efficient and widely used in the literature, for example, to model failure of tibia in Ref. [21]. Linear fractures in the skull were reproduced with element erosion in Ref. [12]. In the second method, the element is not deleted when the failure criteria are met. Upon meeting the failure criteria, the material model is changed to that of a fluid [22,23], i.e., the element loses its ability to carry tensile and shear stresses. In this method, mass and momentum of failed material is retained. However, the loss of shear strength could cause large deformation and ensuing reduction in time-step size. In the third method, the elements are converted into smoothed particle hydrodynamics (SPH) particles when the failure criterion is satisfied. The SPH particles carry mass, velocity, strength, and can interact with neighbor elements. Conversion to particles alleviates the mesh distortion issues with lesser penalty on efficiency and accuracy compared to the previous two methods. This algorithm is commonly employed to model brittle material fracture with reasonable results [24,25]. As a bench mark, we have also included results for the case, where material failure was not modeled by setting the failure strain to a large number. The four cases above are denoted as “erosion,” “fluid,” “SPH,” and “no failure,” respectively, in this paper.