Contrary to the second-order Phase field model (PFM) of fracture, fourth-order PFM provides a more precise representation of the crack surface by incorporating higher-order derivatives (curvature) of the phase-field order parameter in the so-called crack density functional. As a result, in a finite element setting, the weak form of the phase-field governing differential equation requires \(C^1\) continuity in the basis function. \(C^0\) Sibson interpolants or Natural element interpolants are obtained by the ratio of area traced by the second-order Voronoi cell over the first-order Voronoi cells, which is based on the natural neighbor of a nodal point set. \(C^1\) Sibson interpolants are obtained by degree elevating the evaluated \(C^0\) interpolants in the Bernstein-Bezier patch of a cubic simplex. For better computational efficiency while accounting only for the tensile part for driving fracture, a hybrid PFM is adopted. In this work, the numerical implementation of higher-order PFM with \(C^1\) Sibson interpolants along with some benchmark examples are presented to showcase the performance of this method for simulating fracture in brittle materials.
