Improved Staggered Algorithm for Phase-Field Brittle Fracture with the Local Arc-Length Method

Highlights

An improved staggered algorithm that uses the local arc-length method is formulated and proposed.

The energy and damage tolerance convergence criterion are employed in the numerical implementation.

The proposed algorithm results are not dependent on the selected sizes of the loading increments.

1  Introduction

Brittle fracture is one of the most feared failure modes in solids and structures. The numerical simulation method of brittle fracture is significant in engineering applications. Numerical simulation models for crack growth problems are generally divided into discrete crack models and diffuse crack models. The discrete crack model is a geometric description method based on mesh extension, which extends on mesh and mesh boundaries, including the element deletion method [1,2] and the interface element method [3,4]. The element deletion method is a relatively simple fracture modeling method. It only needs to take the stress of the element zero or remove the mass of the fracture element from the overall mass matrix when a particular physical quantity associated with the element reaches the fracture criterion. The element deletion method cannot resolve the crack bifurcation problem [2,5]. The interface element method inserts interface elements at the interfaces between elements, and there is cohesion between nodes of interface elements. When the fracture criterion is reached, the interface element fails and is broken. The interface element method only allows cracks to propagate at the element boundary, with apparent mesh dependence [5] and numerical instability. The diffuse crack model, including the extended finite element method (XFEM) [6–8] and the phase-field fracture method [9,10] is a nongeometric description method that can be separated from the mesh and expanded within the mesh. The XFEM allows cracks to expand inside the element by adding degrees of freedom to the nodes. The XFEM can be used to solve some simple fracture bifurcation problems in combination with existing bifurcation criteria [11–13]. However, it is still challenging to simulate arbitrary crack propagation in three-dimensional problems. The difference between the XFEM and phase-field fracture method lies in the crack surface tracking technology deployed by the complex crack systems simulation method. XFEM is a discontinuous interface method requiring special crack tracking techniques to catch crack surfaces, such as level sets [14] and fast marching [15]. These methods are particularly complicated when addressing three-dimensional and multicrack problems.

The research field of phase-field fracture models can be divided into the field of physics and the field of mechanics. The physical basis of the phase-field fracture models are based on the Ginzburg-Laudau theory, as shown by the model proposed by Hakim et al. [16]. The phase-field fracture model in mechanics is based on the variational fracture approach developed by Francfort et al. [17] and is an extension of the traditional Griffith’s brittle fracture. Bourdin et al. [18,19] realized the rate-independent variational principle of diffuse fracture numerically for the first time. The diffuse phase-field boundary is employed to describe sharp crack boundaries. The phase-field fracture model is described by introducing order parameters and using continuous functions and controls the evolution of order parameters through the phase-field governing equations. Therefore, the phase-field fracture model obtains the path and location of the crack through automatic evolution of the phase-field variables without tracking the geometry of the crack [20,21]. Crack propagation is controlled by the principle of free energy minimization for the phase-field fracture model [22,23].

In terms of algorithm implementation, the phase-field crack model can be solved with either a monolithic algorithm [24,25] or a staggered algorithm [26,27]. In a monolithic algorithm, the displacement field and phase-field are solved simultaneously [28,29]. However, the two fields are calculated separately in the staggered algorithm. Due to the strong nonlinearity of the partial differential equations, however, the monolithic algorithm requires a large number of iteration steps. The elastic strain energy decomposition method based on strain spectrum decomposition enhances the nonlinearity of the equations, and a monolithic algorithm often results in divergent consequences [30–33]. A staggered method was developed by Miehe et al. [34] with a single staggered iterative process. It is widely used due to its simple implementation and strong robustness [35,36]. To alleviate the problem of the common single iteration staggered algorithm, Karlo et al. [37] presented a comprehensive discussion on the use of a stopping criterion within the broadly used staggered algorithm that imposes the Newton-Rapson incremental loading procedure.

Nevertheless, the Newton-Rapson incremental loading procedure has difficulty capturing the complete failure process and the entire load-displacement curve because of the rapid failure of the structure due to the large loading step. To obtain accurate results, however, a small loading increment is necessary for the staggered algorithm. Several loading algorithms can track the whole load-displacement curve well, such as the global arc-length method [38,39], local arc-length method [40], energy dissipation control method [29,41], fracture surface control method [10,42], and J integral control method [43]. Wu [44] proposed an iterative alternate minimization algorithm enhanced with path-following strategies which is based on the fracture surface control and the indirect displacement control. For the construction of complex engineering symmetric curves and surfaces in 2D and 3D phase-field modeling, an efficient method can be considered. Muhammad et al. [45] proposed a novel approach of hybrid trigonometric Bézier curve and applied it to engineering and computer aided design/computer aided manufacturing [46]. The efficient modeling techniques and processes are helpful for computational models.

In this work, the phase-field method based on the rate-independent variational principle of diffuse fracture is considered. The formulations of the staggered algorithm with the local arc-length method are derived and proposed to avoid the sensitivity to the local failure of the structure and the softening of materials. The energy and damage iterative tolerance convergence criterion based on the residual norm of the displacement and phase-field is proposed to save computing resources.

The method of separating the node and degree of freedom of the element is employed to implement the improved staggered algorithm in the ABAQUS package, and they were alternating between the phase-field degree of freedom and displacement degree of freedom. An energy and damage iterative tolerance convergence criterion is embedded in each incremental loading of the local arc-length method. A one element benchmark test compared with the analytical solution was conducted to validate the feasibility and accuracy of the developed method. Several numerical experiments are conducted to simulate crack initiation, intersection, coalescence, branching, and propagation.

The paper is organized as follows. The representation of the phase-field approximation of diffuse crack topology is outlined in Section 2.1. The phase-field formulation for brittle fracture and classical staggered algorithm with Newton-Raphson method is derived in detail in Sections 2.2 and 2.3. The improved staggered solution for phase-filed with the local arc-length method and the implementation in ABAQUS/UEL, which includes the whole solution process, are discussed in detail in Section 3. Several static and dynamic benchmark tests and numerical examples are simulated and compared with other existing results in Section 4. Finally, conclusions are presented in Section 5.

2  Brief Review of the Phase-Field Model for Brittle Fracture

2.1 Representation of the Phase-Field Approximation of Diffuse Crack Topology

First, a simple model is considered to introduce the variational principle of diffuse crack topology. A cracked 1D bar with cross section Γ and length L∈[–∞, +∞] is depicted in Fig. 1a. Assuming that there is a fully opened crack at x = 0, Γ can represent the crack surface that is completely broken. In this work, the auxiliary variables ϕ(x) = [0,1] were employed to describe this sharp topology. When ϕ = 0, broken materials are present, whereas ϕ = 1, unbroken materials are present. The following exponential function [9] is used to approximate the nonsmooth (diffuse) phase-field:

ϕ(x)=exp(−|x|/l0)(1)

where l0 is a regularization parameter describing the actual width of a diffuse crack. The axial crack in Fig. 1b can be transferred into Fig. 1c employing the phase-field approximation of the diffuse crack. Regularization is determined by l0. According to Fig. 1b, considering the Dirichlet-type boundary conditions ϕ(0) = 1 and ϕ(±∞) = 0. The solution of the homogenous PDE can be written as:

ϕ(x)−l02ϕ′′(x)=0(2)

The one-dimensional description of the diffuse crack topology can be directly extended to the multidimensional case. According to the derivation process, the surface density function γ(ϕ,▽ϕ) of cracks in unit volume can be expressed as:

γ(ϕ,∇ϕ)=12l0ϕ2+l02|∇ϕ|2(3)

Figure 1: Sharp and diffuse crack topology (a) Cracked 1D bar (b) Sharp crack (c) Diffuse crack

Assume that we are given a sharp crack surface topology in Ω, as shown in Fig. 2a. The regularized crack phase-field in the domain Ω can be obtained according to the minimization principle below (as shown in Fig. 2b).

Figure 2: (a) Discontinuous discrete crack surface (b) Diffuse crack surface expressed by the phase-field variable

Francfort et al. [17] proposed a variational fracture model based on Griffith’s energy theory. The total potential energy of elastic solids consists of two parts: elastic energy Eb and crack surface energy Es:

E(u,Γ)=Eb+Es=∫Ωψe(ε(u))dΩ+∫ΓGcdΓ(4)

where ψe(ε(u)) is the elastic strain energy density function, ε(u) is the strain tensor, u∈Rd (d∈{1,2,3}) and Gc is the critical energy release rate. The variational model considers that at time t∈[0, T], the initiation, propagation, and branching of cracks at any location should minimize the total potential energy and meet some irreversible condition; namely, once a crack has formed, it cannot be restored. The regularized surface energy can be written as:

∫ΓGcdΓ=∫ΩGcγ(ϕ,∇ϕ)dΩ(5)

Assuming that linear elasticity exists in the unbroken solid, the elastic energy density can be calculated as:

ψ0(ε)=12εT⋅C0⋅ε(6)

where C0 is the material linear elastic stiffness matrix. Elastic energy degradation due to damage follows the function:

g(ϕ)=(1−ϕ)2+k(7)

the elastic energy density can be rewritten as:

ψe(ε)=g(ϕ)⋅12εT⋅C0⋅ε(8)

According to the above fracture variation criterion, the surface energy of cracks results from elastic strain energy, that is, elastic strain energy drives the evolution of the phase-field variables. If the elastic strain energy is not decomposed, it will cause false branching of the crack. Given this problem, this work, based on the method proposed by Miehe et al. [23], observes the tension and compression decomposition of the elastic strain energy occurring so that the tension component of the elastic strain energy drives the evolution of the phase-field. The spectral decomposition of the strain tensor and elastic energy density is carried out via:

ψe±=12λ⟨tr(ε)⟩±2+μtr(ε±2)(9)

where λ and μ are Lame constants. To correctly calculate the reduction effect of the phase-field variables on the stiffness of the material, assuming that the phase-field variables only affect the stretched part of the elastic strain energy density, the elastic strain energy density can be expressed as:

ψe(ε)=[(1−ϕ)2+k]ψe+(ε)+ψe−(ε)(10)

where the k ∈ [0, 1) is a model parameter used to avoid numerical singularity when ϕ = 1.

Considering the crack surface energy expression in Eq. (3) and the elastic strain energy density expression in Eq. (6), and considering the kinetic energy Ek as:

Ek=∫Ω12ρu˙⋅u˙dV(11)

We can substitute part of the energy into the variational principle expression, to obtain:

∫Ω(ρd2udt2−∇⋅σ−γ¯)⋅δu⋅dV+∫Ω[Gclc(ϕ−l02∇2ϕ)−2(1−ϕ)ψ0]⋅δϕ⋅dV+∫ΓΩ(σ⋅n−η¯)⋅δu⋅dΓΩ+∫Γd(Gcl0∇ϕ⋅n)⋅δd⋅dΓϕ=0(12)

where γ¯ and η¯ denote the body and face forces, respectively, and n is the unit outward normal vector. The minimization of the internal potential energy relative to the phase-field is now considered. In order to introduce the irreversible condition of damage, the historical field variable H need to be introduced:

H={ψ0(ϵ) if ψ0(ϵ)>HtHt otherwise (13)

Ht is the previously calculated energy obtained in history according to the displacement field ut at the time increment t.

By defining the change in external and internal potential energy and introducing the principle of virtual displacement, the governing equations of the model are obtained:

Displacement-field:  {∇⋅σ+γ¯=ρϕ2u/ϕt2in Ωu=u¯ on Γuσ⋅n=γ¯on Γγ(14)

Phase-field:{Gc(1l0d−l0∇2ϕ)=2(1−ϕ)H in Ω∇ϕ⋅n=0 on Γϕ(15)

In order to obtain the numerical solution of the partial differential Eqs. (13) and (14), the finite element method is used to deal with the weak form more conveniently.

∫Ω{σδε−b⋅δu}dV+∫∂ΩHH⋅δudA=0(16)

∫Ω{−2(1−ϕ)δϕψ0(ε)+Gc[1ℓϕδϕ+l0∇ϕ⋅∇δϕ]}dV=0(17)

2.2 Discretization

In 2D space, the displacement field and the phase-field can be discretized as:

u=∑i=1mNiuui,ϕ=∑i=1mNiϕi(18)

where the shape function matrix is expressed as:

Niu=[Ni00Ni](19)

then it is possible to express the gradients:

ε=∑i=1mBiuui(20)

∇ϕ=∑i=1mBiϕϕi(21)

The strain-displacement matrices are expressed as:

Biu=[Ni,x00Ni,yNi,yNi,x],Biϕ=[Ni,xNi,y](22)

According to the above expression, the residual of the discrete equation corresponding to the equilibrium condition for the displacement field and the crack phase-field can be expressed as:

riu=∫Ω[(1−ϕ)2+k](Biu)Tσ0dV−∫Ω(Niu)TbdV−∫∂ΩH(Niu)THdA(23)

riϕ=∫Ω{−2(1−ϕ)Niψ0(ϵ)+Gc[1l0Niϕ+l0(Biϕ)T∇ϕ]}dV(24)

To obtain the solutions of ru = 0 and rϕ = 0, and considering that their residuals are nonlinear, the incremental iteration scheme of the Newton-Raphson method is adopted.

{uϕ}t+Δt={uϕ}t−[KuuKuϕKϕuKϕϕ]t−1{RuRϕ}t(25)

where the static tangent matrices are calculated as:

Kijuu=∂riu∂uj=∫Ω[(1−ϕ)2+k](Biu)TC0BjudVKijuϕ=∂riu∂ϕj=∫Ω−2(1−ϕ)(Biu)Tσ0NjdVKijϕu=∂riϕ∂uj=∫Ω−2(1−ϕ)Niσ0TBjudVKijϕϕ=∂riϕ∂ϕj=∫Ω{[2ψ0(ε)+Gcl0]NiNj+Gcl0(Biϕ)T(Bjϕ)}dV(26)

2.3 Classical Staggered Algorithm with Newton-Raphson Method

In the case of unstable crack growth, the monolithic model is prone to convergence problems. The stress fields are rearranged due to the newly degenerated stiffness matrix when the crack propagates. Therefore, the solver cannot obtain a stable equilibrium solution. To obtain a stable model, the displacement field u and phase-field ϕ of the coupled system can be solved as sequential coupled staggered fields. The phase-field and displacement fields in the staggered scheme are coupled variable fields. Furthermore, the calculation of the current history field is the main idea for the algorithmic split of the coupled equations, and it enforces the irreversibility of the evolution of the crack phase-filed. Consequently, the history field satisfies the Kuhn-Tucker conditions,

ψ0−H≤0,H˙≥0,H˙(ψ0−H)=0(27)

for both loading and unloading cases. Thus, no penalty term is needed in contrast to the monolithic approach. The residual corresponding to the evolution of the crack phase field is then as follows:

riϕ=∫Ω{−2(1−ϕ)NiH+Gc[1l0Niϕ+l0(Biϕ)T∇ϕ]}dV(28)

Similarly, the residual corresponding to the displacement field is then as:

riu=∫Ω[(1−ϕ)2+k](Biu)Tσ0dV−∫Ω(Niu)TbdV−∫∂ΩH(Niu)THdA(29)

The following system of the equation can be solved iteratively by employing the Newton-Raphson method,

{uϕ}t+Δt={uϕ}t−[Kuu00Kϕϕ]t−1{rurϕ}t(30)

in which the tangent stiffness matrices are calculated as:

Kijuu=∂riu∂uj=∫Ω[(1−ϕ)2+k](Biu)TC0BjudVKijϕϕ=∂riϕ∂ϕj=∫Ω{[2H+Gcl0]NiNj+Gcl0(Biϕ)T(Bjϕ)}dV(31)

The new formula generates the above decoupling equations for sequential phase and displacement fields. Then, a staggered scheme of phase-field model is proposed. The displacement, phase and history fields ut, ϕt and Ht are known at the previous time increment t. Then, the prescribed load vectors b are updated at the current time increment t + Δt. The current phase-field ϕt + Δt and the current displacement field ut + Δt are now computed using the corresponding linear algebraic Euler equations.

Ktϕϕϕt+Δt=−rtϕKtuuut+Δt=−rtu(32)

3  Improved Staggered Solution Scheme with the Local Arc-Length Method

In this section, the improved staggered algorithm with the local arc-length method formulations and the corresponding numerical implementation are established. The local arc-length method [40] is applied to combine the staggered solution scheme with the energy and damage iterative tolerance convergence criterion.

3.1 Formulations of the Local Arc-Length Method

The arc-length method adjusts the external force load f or displacement load χ to the convergence trend by constructing a scaling parameter λ during the loading process. The size of the loading step can be expressed in simple linear form as follows:

f=λf^+f0,χ=λχ^+χ(33)

The f^ and χ^ are the scaled external force load and displacement load, respectively, and f0 and λ0 correspond to the initial values of f and χ. For the displacement load of the arc-length method that requires f^ = 0 and χ^ ≠ 0, the force load does the opposite. A constraint equation is added to meet the solution requirements because λ is an unknown parameter; for the local arc-length method,

Rλ=Δl−SΔχ=0(34)

where S is the weight function.

After linearization by the Newton-Raphson method, the residual Eqs. (25) and (26) are transformed into:

[KuuKuλKλuKλλ]{δχδλ}={RuRλ}(35)

where the expressions of tangent stiffness are as follows:

Kuu=∫ΩBT(∂σ∂ε)BdVKuλ=−∂Ru∂λ=−f^Kλu=−∂Rλ∂χ=SKλλ=−∂Rλ∂λ=0(36)

The unknown increment can be obtained by the following formula:

{δχδλ}=[KuuKuλKλuKλλ]{RuRλ}(37)

According to the Sherman–Morrison formula, Eq. (30) can transfer to:

{δχδλ}={χI0}−KλuχI+RλKλuχII+Rλλ{χII0}(38)

KuuχI=Ru,KuuχII=f^(39)

where χI and χII are obtained from the following equilibrium equations and correspond to the displacement of the node caused by residuals Ru and the external force f, respectively.

3.2 Staggered Algorithm with the Local Arc-Length Method and Its Numerical Implementation

Using the Newton-Raphson backward Euler scheme means that the numerical solution is unconditionally stable; that is, to achieve convergence means to obtain an equilibrium solution. Therefore, the solution is independent of the time increment. In the classic staggered scheme, the phase and the displacement fields are solved independently. This is unlike the monolithic case, where convergence is fulfilled for both coupled fields simultaneously, retaining unconditional stability. By solving for the phase and the displacement fields independently, the staggered scheme can track the equilibrium solution in a semi-implicit method, leading to a more robust numerical implementation. However, load increment sensitivity studies should be conducted. The residual staggered algorithm is improved to reduce the influence of the loading increment step on the solution process. The flowchart of the improved staggered solution scheme with the local arc-length method is shown in Fig. 3.

Figure 3: Process for the improved staggered algorithm with the arc-length method based on the energy and damage convergence criterion

Expression (37) is solved for the phase-field and the displacement field. The displacement field can be solved by the local arc-length loading method. The loading increment and iteration step should be distinguished in solving the displacement field. In some nonlinear analyses, the total load of each analysis step is decomposed into smaller incremental steps. The convergence results are obtained under a small incremental load to carry out the next step. The iteration step attempts to find a balanced solution in the loading increment step. At time increment t, Ktuu and Rtu are built from ϕt, and Ktϕϕ and Rtϕ are built from ut. When load increment step n is applied to load increment step n + 1, Rtu(ϕnk) can be calculated in the previous load increment step:

R(ϕnk)=Rtu(ϕnk−1)+Rtϕ(ϕnk)(40)

and rtu(ϕnk) can be obtained in the form:

Rtϕ(ϕnk)=Rtϕ(ϕnk−1)+Rtu(ϕnk−1)(41)

in the n + 1 load increment step, residuals are expressed as:

Rtu(ϕn+1k)=Rtu(ϕn+1k−1)+Rtϕ(ϕnk)(42)

Rtϕ(ϕn+1k)=Rtϕ(ϕn+1k−1)+Rtu(ϕnk)(43)

Then, following this process and repeating the iteration by assembling the residual stiffness matrix Ktuu,Ktϕϕ and combining Rtu and Rtϕ, ignores the effect of the load increment step directly. We obtain further verification from the numerical experiments section of this work.

By introducing the energy and damage convergence criterion, the improved staggered solution scheme is implemented to further improve the computational efficiency in solving phase-field brittle fracture problems. In the computation of the staggered solution, the solutions of the displacement field and phase-field both have their convergence criteria.

At the beginning of the calculation, the damage value is fixed to solve the displacement subproblem. The damage value can be solved by the local arc-length method. When the displacement field converges, the displacement field is taken as a known quantity, and the damage subproblem is solved. There are two boundary conditions for solving the damage subproblem. One is that the phase-field is bounded, which is ϕ ∈ [0,1], and the range of the unknown quantity needs to be constrained in the process of solving. The other is that to ensure the irreversible condition of damage evolution d˙ ≥ 0, the lower bound of the damage solution must be monotonically increasing. There are two kinds of convergence criteria adopted in this paper, one is based on the damage, and the other is based on energy:

|Ei+1−Ei|Ei≤TOL(44)

|ϕi+1−ϕi|≤TOL(45)

The process for the improved staggered algorithm based on the double convergence criterion is shown in Fig. 3. The Expression (44) is the energy convergence criterion, while Expression (45) is the damage convergence criterion. The improved staggered algorithm considers both the energy convergence criterion and the damage convergence criterion. Once one of the convergence criteria reaches the corresponding tolerance, the program will continue to calculate the next iteration step. TOL is the tolerance of convergence, generally 10–5. For materials with strong brittleness, TOL should be lower; otherwise, the peak point may be lost. In the actual solution process, when the structure is in the linear elastic state, there is only one staggered solution process for each loading step. However, when the structural response enters the nonlinear state, because of the strict convergence criterion of the staggered solution algorithm, multiple staggered solution iterations are often needed in one loading step.

The core of the traditional staggered solution scheme is to fix the phase-field ϕ and solve for the displacement u, fix the displacement u and solve for the phase-field ϕ, and repeat the cycle until convergence. In addition to the conventional displacement degree of freedom, the degree of freedom of the finite element node needs to be added to the phase-field, and it is realized from the perspective of a multifield finite element. The alternating process of the phase-field degree of freedom and displacement degree of freedom can be realized in the following three ways:

1.    Add the phase-field element and displacement element. The corresponding degrees of freedom are given in the staggered solution process.

2.    Add the phase-field element. Two submodels with the same mesh division should be established, one of which is the displacement element and the other being the phase-field element. In the staggered solution scheme, submodels are called for solutions.

3.    Separate the node and the degree of freedom of the element. Set the element type without the degree of freedom in the model file, and then give the element’s corresponding node the degree of freedom in the solution process.

Regarding program expansibility and simplicity, the third method is best. The phase-field brittle fracture model can be solved directly in the finite element framework. The commercial software ABAQUS is employed to solve the iteration process and realize the improved staggered algorithm for phase-field modeling.

The layered system of the finite element implemented in ABAQUS with UEL and UMAT is shown in Fig. 3. The process for solving the phase-field model in each layer is described in detail in Fig. 4.

Figure 4: Schematic illustration of the 3-layered element structure in ABAQUS with its subroutine

Since the phase-field equation and the balanced equation are solved separately, the implicit solver ABAQUS/Standard or the explicit solver ABAQUS/Explicit can be flexibly selected according to the model when solving the displacement field. The former is implemented through the UMAT subroutine interface provided by ABAQUS, and the latter is implemented through VUMAT.

As depicted in Fig. 4, the node and the integration points coincide in the phase-field model. The first layer of finite elements includes user elements implemented by the user subroutine UEL. The degree of freedom ϕi, is implemented instead of the ABAQUS rotation degree of freedom so that the phase-field residual norm and the residual displacement norm are entirely separated. The discretization matrices, stiffness matrix, and residual vector must be defined in the first layer. The second layer of finite element is for the displacement field.

The solution-dependent variables used to solve differential equations in two dimensions are listed in Table 1. The third layer is for the UMAT field element. Not only can it be used as a visualization tool, but it can also control the convergence of the iteration process.

4  Numerical Examples and Algorithm Validation

4.1 One Element Test

A finite element model with the UEL implementation should be developed to verify the feasibility of solving brittle fracture problems using the algorithm, as demonstrated through the most straightforward case in Fig. 5. A 2D plate is extended on top in direction Y, whereas bottom nodes are constrained in both directions. As the dimensions of the square plate are 1 × 1 mm, the characteristic of the element is 1 mm. The material and fracture properties are chosen: Young’s modulus E = 210 GPa, Poisson’s μ = 0.3, critical fracture energy density Gc = 2.7 × 10–3 kN/mm, and length scale parameter l0 = 0.1 mm.

Figure 5: The plane strain element with boundary condition

This case is solved using 5 different displacements Δu to reach u = 0.1 mm. The five-step displacements are Δu = 1 × 10–5, Δu = 5 × 10–5, Δu = 1 × 10–4, Δu = 5 × 10–4, and Δu = 1 × 10–3, corresponding to 10000 steps, 2000 steps, 1000 steps, 200 steps, and 100 steps. Fig. 6 shows the axial stress-strain response comparison of the numerical results of the traditional staggered algorithm and the improved staggered algorithm against the analytical solution.

Figure 6: Stress-strain curves obtained by (a) The traditional staggered algorithm and (b) The improved staggered algorithm in comparison with the analytical solution

According to Fig. 6a, when the axial strain is close to 0.01, the axial stress reaches the maximum value. With the increase of strain, the error of stress increases in the rising stage and decreases in the falling stage. At the peak of the curve, the error of stress is at its maximum. The traditional staggered algorithm is strictly dependent on the size of the increment step. The smaller the increment is, the smaller the error between the numerical result and the analytical solution. That is when the crack starts to propagate and due to the quick stiffness and stress reduction in some elements, stability problems appear in the traditional staggered algorithm. According to Fig. 6b, however, the numerical result of the improved staggered algorithm is independent of the increment step size. Regardless of the value of the increment step, the numerical solution can be in good agreement with the analytical solution. Similarly, according to Fig. 7, the numerical solution of the governing phase-field evolution through the axial strain matches well with the analytical solution. Hence, the improved staggered algorithm is validated through the one element test. The results show that the improved staggered algorithm for phase-field brittle fracture with the local arc-length method efficiently solves brittle fracture problems.

Figure 7: Comparison of the phase-field between the numerical and analytical results

4.2 Single Edge Notched Test

The single edge notched test is the most common test used to validate phase-field brittle fracture models. The geometry and the boundary conditions of the square specimen are depicted in Fig. 8. The side length of the square is 1 mm, and the edge notched is horizontal from the middle of the left side to the center of the square. α = 90° and the α = 0°, corresponds to the pure tension case and pure shear case, respectively. The material and fracture properties are chosen: Young’s modulus E = 210 GPa, Poisson’s μ = 0.3, critical fracture energy density Gc = 2.7 × 10–3 kN/mm, and length scale parameter l0 = 0.075 mm.

Figure 8: Geometry and boundary conditions of a single edge notched specimen

Fig. 9 shows the fracture pattern for unidirectional tension (α = 90°). Fig. 10 shows the fracture pattern for pure shear (α = 0°). The two fracture patterns in the improved staggered algorithm are in good agreement with the existing references [35]. For unidirectional tension (α = 90°), the crack begins in the center of the geometry and extends horizontally to the right side. For pure shear (α = 0°), a crack path can be seen starting with β = 61° from the deformation direction.

Figure 9: Crack patterns for the single notch edge tension test at different displacements

Figure 10: Crack patterns for the single notch edge shear test at different displacements

Fig. 11 shows the reaction force in tensile specimens, along with the results of Miehe et al. [9]. The maximum reaction force value is consistent, and only a small deviation can be observed during crack propagation. In the rising stage, the error of the reaction force is small. The large error in the descending stage mainly comes from the selection of mesh density.

Figure 11: Reaction force in the single notched edge tension test

In addition to the quasi-static crack propagation of the 2D model, the 3D model of the single edge notch test is solved using the same parameter. The geometry and boundary conditions for the 3D single edge notched tensile modeling are depicted in Fig. 12. The number of solution dependent state variables of two dimensions are different from the three dimensions due to the z-direction degree of freedom, which is defined as 26 SDVs. The crack starts at the end of the initial crack in the middle of the specimen and ends at the right boundary of the specimen, similar to the two-dimensional model. The crack patterns for the single notch edge tension test of the 3D model at different displacements are shown in Fig. 13. The case shows that the improved staggered algorithm for phase-field brittle fracture with the local arc-length method is efficient for solving brittle fracture problems in pure tension and shear situations in the 2D model and the 3D model.

Figure 12: Geometry and boundary conditions for the 3D single edge notched tensile test

Figure 13: Crack patterns for the single notch edge tension test of the 3D model at different displacements

4.3 Notched Bimaterial Tensile Test

This case study is a tensile sample of two materials. The upper material is harder, while the lower material is 10 times softer and has an initial notch. The sample is tensioned parallel to the separation line. The geometry and boundary conditions are shown in Fig. 14. The upper material properties are chosen: Young’s modulus E = 210 GPa, Poisson’s μ = 0.2, and critical fracture energy density Gc = 1 × 10–2 kN/mm. The lower material properties are chosen: Young’s modulus E = 3.77 GPa, Poisson’s μ = 0.2, and critical fracture energy density Gc = 1 × 10–3 kN/mm.

Figure 14: The geometry and boundary conditions for the bimaterial tensile test

The purpose of this benchmark test is to show that our implementation can model the crack branching phenomenon. Fig. 15 shows the crack patterns for the bimaterial tensile test at different displacements. After the crack begins, it normally propagates until it reaches the material transition zone. Then, because it takes less energy to branch over a short period of time than it does to continue into the upper hard material, the damage zone becomes larger, and the cracks separate along with the interface between the upper and lower material.

Figure 15: Crack patterns for the bimaterial tensile test at different displacements

Since the purpose of this benchmark test is to show crack branching, the center of the finite element mesh is dense. As the local amplification region is showm in Fig. 15, the cracks spreads along the direction of the crack, and the width of the crack decreases after separating along with the interface. As a result, cracks in harder materials are wider because the finite element size in this region is larger. This notched bimaterial tensile numerical experiment is conducted to simulate crack branching using the improved staggered algorithm for phase-field brittle fracture with the local arc-length method.

4.4 Symmetric and Asymmetric Double Notched Tensile Test

After the simple example, we will show a few more complex cases to demonstrate the utility of the new implementation. A symmetric and asymmetric double-notched specimen is used to investigate coalescence of cracks. Figs. 16 and 18 depict the exact geometry. The following material properties are used: Young’s modulus E = 210 GPa, Poisson’s μ = 0.3, critical fracture energy density Gc = 2.7 × 10–3 kN/mm, and length scale parameter l0 = 0.075 mm.

Figure 16: Geometry and boundary conditions of symmetric doubled notched specimen

Figure 17: Crack patterns for symmetric doubled notched specimen tension test at different displacements

Figure 18: Geometry and boundary conditions of asymmetric doubled notched specimen

The crack patterns for symmetric and asymmetric doubled notched specimen tension tests at different displacements are shown in Figs. 17 and 19. For asymmetric doubled notched specimen tension tests, the cracks expand simultaneously from the left and right sides. The upper crack extends downwards and the lower crack extends upwards and then meet in the center of the specimen. For symmetric doubled notched specimen tension tests, the cracks expand simultaneously from the left and right side, then meet in the middle and merge. The final fracture pattern is very consistent with the previous results. This symmetric and asymmetric double notched numerical experiment is conducted to simulate crack coalescence using the improved staggered algorithm for phase-field brittle fracture with the local arc-length method.

Figure 19: Crack patterns for the asymmetric doubled notched specimen tension test at different displacements

4.5 L-Shaped Specimen

The L-shaped specimen in [26] is taken as an example to simulate the evolution process of the crack path. The geometry and boundary conditions are shown in Fig. 20. The following material properties are used: Young’s modulus E = 28.8 GPa, Poisson’s μ = 0.18, critical fracture energy density Gc = 8.9 × 10–5 kN/mm, and length scale parameter l0 = 1.1185 mm. A finite element mesh consisting of 52634 elements was used to refine the expected crack propagation region. The improved staggered algorithm is used for this specimen to demonstrate its ability in more complex geometry.

Figure 20: Geometry and boundary conditions of the L-shaped specimen

Fig. 21 shows the phase-field and equivalent strain plots for the L-shaped specimen obtained from the phase-field fracture model at different displacements. The crack starts at the corner of the L-shaped plate and spreads slightly upward horizontally. The crack path is basically consistent with the path of equivalent strain plots. The results were compared with the experimental measurement results provided in [26]. The calculated results of the crack growth process are convergent and consistent with the published experimental results. This L-shaped specimen numerical experiment is conducted to simulate crack propagation in more complex modeling using the improved staggered algorithm for phase-field brittle fracture with the local arc-length method.

Figure 21: Phase-field and equivalent strain plots for the L-shaped specimen obtained from the phase-field fracture model at different displacements. (a) u = 0.20 mm, (b) u = 0.35 mm, (c) u = 0.70 mm

4.6 Dynamic Crack Branching

In addition to the quasi-static crack propagation examples, dynamic crack propagation is also handled by means of the improved staggered algorithm. The geometry and boundary conditions of the dynamic crack branching are shown in Fig. 22. An initial crack with a length of 50 mm was prefabricated in the middle of the plate, and the top and bottom sides of the plate were subjected to an impact load of 1 MPa. In order to reduce the errors of mesh size on the calculation results, the phase-field fracture modeling is composed of uniform quadrilateral elements with a size of 2.5 × 10–4. The material and fracture properties are chosen as: density ρ = 2450 kg/m3, Young’s modulus E = 32 GPa, Poisson’s μ = 0.18, critical fracture energy density Gc = 3 × 10–3 kN/mm, length scale parameter l0 = 0.01 mm, and Rayleigh wave velocity VR = 2125 m/s.

Figure 22: Geometry and boundary conditions for the dynamic crack problems

The crack patterns for the dynamic crack branching are shown in Fig. 23, and the velocity curve of the crack is shown in Fig. 24. The crack will branch when it suffers the impact load. From the curve, the maximum velocity of the crack tip reaches approximately 0.5 VR, which is basically consistent with the results of the test. It is noted that the crack tip velocity branches immediately after reaching the peak for the second time, and the branching occurs from 30 μs and 36 μs. The velocity of the crack tip decreases immediately after crack branching and starts to increase again after 40 μs. Borden et al. [28] first applied the phase-field method to model the case. Song et al. adopted other fracture numerical methods for the same model. In this work, we compare the curve of the velocity of the crack tip propagation with the results of Borden’s test, which is solved in terms of the monolithic algorithm. The numerical results of the algorithm are in good agreement, and the velocity of the crack tip for the staggered method is slightly lower than that of the monolithic method. The phase-field fracture model does not require additional branch criteria or numerical iteration of discontinuous displacement fields when simulating dynamic crack branching and propagation. Compared with the XFEM, phase-field fracture modeling is more efficient. This dynamic branching numerical experiment is conducted to simulate dynamic crack propagation using the improved staggered algorithm, which applies not only to the quasi-static model but also to the dynamic model.

Figure 23: Crack patterns for the dynamic crack problems

Figure 24: The crack tip velocity for the dynamic crack problems

5  Conclusions

In this study, an improved staggered scheme with the local arc-length method is proposed to study phase-field brittle fracture. The energy and damage iterative tolerance convergence criteria based on the residuals of displacement and phase-field is developed and implemented in ABAQUS. We adopt a one element benchmark numerical test and several numerical examples to verify the feasibility and accuracy of the phase-field model. From the presented work, the following conclusions can be drawn:

1.    The local arc-length method applied to the staggered solution scheme of the phase-field brittle fracture can track the crack propagation process and load-displacement curve well.

2.    The improved staggered scheme for brittle fracture is not sensitive to the size of the load increment. The accuracy and efficiency are greatly improved.

3.    The results of one element benchmark test and several complex numerical examples are in good agreement with the other numerical solutions.

4.    The improved staggered scheme for the brittle fracture model is an effective method to accurately simulate crack initiation, coalescence, intersection, and branching.

Funding Statement: Financial supports by the National Key R&D Program of China (No. 2018YFD1100401), the National Natural Science Foundation of China (No. 51578142), the Fundamental Research Funds for the Central Universities (No. LEM21A03), and Jiangsu Key Laboratory of Engineering Mechanics (Southeast University) are gratefully acknowledged.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

 1.  Fan, R., Fish, J. (2008). The Rs-method for material failure simulations. International Journal for Numerical Methods in Engineering, 73(11), 1607–1623. DOI 10.1002/(ISSN)1097-0207. [Google Scholar] [CrossRef]

 2.  Liu, Y., Filonova, V., Hu, N., Yuan, Z., Fish, J. et al. (2014). A regularized phenomenological multiscale damage model. International Journal for Numerical Methods in Engineering, 99(12), 867–887. DOI 10.1002/nme.4705. [Google Scholar] [CrossRef]

 3.  Xu, X. P., Needleman, A. (1994). Numerical simulations of fast crack-growth in brittle solids. Journal of the Mechanics and Physics of Solids, 42(9), 1397–1434. DOI 10.1016/0022-5096(94)90003-5. [Google Scholar] [CrossRef]

 4.  Camacho, G. T., Ortiz, M. (1996). Computational modelling of impact damage in brittle materials. International Journal of Solids and Structures, 33(20–22), 2899–2938. DOI 10.1016/0020-7683(95)00255-3. [Google Scholar] [CrossRef]

 5.  Song, J. H., Wang, H., Belytschko, T. (2008). A comparative study on finite element methods for dynamic fracture. Computational Mechanics, 42(2), 239–250. DOI 10.1007/s00466-007-0210-x. [Google Scholar] [CrossRef]

 6.  Wen, L., Tian, R. (2016). Improved XFEM: Accurate and robust dynamic crack growth simulation. Computer Methods in Applied Mechanics and Engineering, 308, 256–285. DOI 10.1016/j.cma.2016.05.013. [Google Scholar] [CrossRef]

 7.  Monteiro, A. B., Wolenski, A. R. V., Barros, F. B., Pitangueira, R. L. S., Penna, S. S. (2017). A computational framework for G/XFEM material nonlinear analysis. Advances in Engineering Software, 114(9), 380–393. DOI 10.1016/j.advengsoft.2017.08.002. [Google Scholar] [CrossRef]

 8.  Feng, S. Z., Li, W. (2018). An accurate and efficient algorithm for the simulation of fatigue crack growth based on XFEM and combined approximations. Applied Mathematical Modelling, 55, 600–615. DOI 10.1016/j.apm.2017.11.015. [Google Scholar] [CrossRef]

 9.  Miehe, C., Welschinger, F., Hofacker, M. (2010). Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field Fe implementations. International Journal for Numerical Methods in Engineering, 83(10), 1273–1311. DOI 10.1002/nme.2861. [Google Scholar] [CrossRef]

10. Singh, N., Verhoosel, C. V., de Borst, R., van Brummelen, E. H. (2016). A fracture-controlled path-following technique for phase-field modeling of brittle fracture. Finite Elements in Analysis and Design, 113, 14–29. DOI 10.1016/j.finel.2015.12.005. [Google Scholar] [CrossRef]

11. Belytschko, T., Chen, H., Xu, J. X., Zi, G. (2003). Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. International Journal for Numerical Methods in Engineering, 58(12), 1873–1905. DOI 10.1002/(ISSN)1097-0207. [Google Scholar] [CrossRef]

12. Menouillard, T., Belytschko, T. (2010). Dynamic fracture with meshfree enriched XFEM. Acta Mechanica, 213(1–2), 53–69. DOI 10.1007/s00707-009-0275-z. [Google Scholar] [CrossRef]

13. Xu, D., Liu, Z., Liu, X., Zeng, Q., Zhuang, Z. (2014). Modeling of dynamic crack branching by enhanced extended finite element method. Computational Mechanics, 54(2), 489–502. DOI 10.1007/s00466-014-1001-9. [Google Scholar] [CrossRef]

14. Gravouil, A., Moes, N., Belytschko, T. (2002). Non-planar 3D crack growth by the extended finite element and level sets–Part II: Level set update. International Journal for Numerical Methods in Engineering, 53(11), 2569–2586. DOI 10.1002/nme.430. [Google Scholar] [CrossRef]

15. Sukumar, N., Chopp, D. L., Bechet, E., Moes, N. (2008). Three-dimensional non-planar crack growth by a coupled extended finite element and fast marching method. International Journal for Numerical Methods in Engineering, 76(5), 727–748. DOI 10.1002/nme.2344. [Google Scholar] [CrossRef]

16. Hakim, V., Karma, A. (2009). Laws of crack motion and phase-field models of fracture. Journal of the Mechanics and Physics of Solids, 57(2), 342–368. DOI 10.1016/j.jmps.2008.10.012. [Google Scholar] [CrossRef]

17. Francfort, G. A., Marigo, J. J. (1998). Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids, 46(8), 1319–1342. DOI 10.1016/S0022-5096(98)00034-9. [Google Scholar] [CrossRef]

18. Bourdin, B., Francfort, G. A., Marigo, J. J. (2000). Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids, 48(4), 797–826. DOI 10.1016/S0022-5096(99)00028-9. [Google Scholar] [CrossRef]

19. Bourdin, B., Francfort, G. A., Marigo, J. J. (2008). The variational approach to fracture. Journal of Elasticity, 91(1–3), 5–148. DOI 10.1007/s10659-007-9107-3. [Google Scholar] [CrossRef]

20. Fang, J., Wu, C., Rabczuk, T., Wu, C., Ma, C. et al. (2019). Phase field fracture in elasto-plastic solids: Abaqus implementation and case studies. Theoretical and Applied Fracture Mechanics, 103, 102252. DOI 10.1016/j.tafmec.2019.102252. [Google Scholar] [CrossRef]

21. Zhou, S., Zhuang, X., Zhu, H., Rabczuk, T. (2018). Phase field modelling of crack propagation, branching and coalescence in rocks. Theoretical and Applied Fracture Mechanics, 96(1), 174–192. DOI 10.1016/j.tafmec.2018.04.011. [Google Scholar] [CrossRef]

22. Amor, H., Marigo, J. J., Maurini, C. (2009). Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments. Journal of the Mechanics and Physics of Solids, 57(8), 1209–1229. DOI 10.1016/j.jmps.2009.04.011. [Google Scholar] [CrossRef]

23. Miehe, C., Lambrecht, M. (2001). Algorithms for computation of stresses and elasticity moduli in terms of Seth-Hill's family of generalized strain tensors. Communications in Numerical Methods in Engineering, 17(5), 337–353. DOI 10.1002/cnm.404. [Google Scholar] [CrossRef]

24. Kuhn, C., Muller, R. (2010). A continuum phase field model for fracture. Engineering Fracture Mechanics, 77(18), 3625–3634. DOI 10.1016/j.engfracmech.2010.08.009. [Google Scholar] [CrossRef]

25. Amiri, F., Millan, D., Shen, Y., Rabczuk, T., Arroyo, M. (2014). Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 69(13), 102–109. DOI 10.1016/j.tafmec.2013.12.002. [Google Scholar] [CrossRef]

26. Borden, M. J., Hughes, T. J. R., Landis, C. M., Verhoosel, C. V. (2014). A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Computer Methods in Applied Mechanics and Engineering, 273(1), 100–118. DOI 10.1016/j.cma.2014.01.016. [Google Scholar] [CrossRef]

27. Borden, M. J., Verhoosel, C. V., Scott, M. A., Hughes, T. J. R., Landis, C. M. (2012). A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering, 217(8), 77–95. DOI 10.1016/j.cma.2012.01.008. [Google Scholar] [CrossRef]

28. Vignollet, J., May, S., de Borst, R., Verhoosel, C. V. (2014). Phase-field models for brittle and cohesive fracture. Meccanica, 49(11), 2587–2601. DOI 10.1007/s11012-013-9862-0. [Google Scholar] [CrossRef]

29. Verhoosel, C. V., Remmers, J. J. C., Gutierrez, M. A. (2009). A dissipation-based arc-length method for robust simulation of brittle and ductile failure. International Journal for Numerical Methods In Engineering, 77(9), 1290–1321. DOI 10.1002/nme.2447. [Google Scholar] [CrossRef]

30. Liu, G., Li, Q., Msekh, M. A., Zuo, Z. (2016). Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field model. Computational Materials Science, 121(4), 35–47. DOI 10.1016/j.commatsci.2016.04.009. [Google Scholar] [CrossRef]

31. Kristensen, P. K., Martinez-Paneda, E. (2020). Phase field fracture modelling using Quasi-Newton methods and a new adaptive step scheme. Theoretical and Applied Fracture Mechanics, 107(8), 102446. DOI 10.1016/j.tafmec.2019.102446. [Google Scholar] [CrossRef]

32. Wu, J. Y., Huang, Y. (2020). Comprehensive implementations of phase-field damage models in Abaqus. Theoretical and Applied Fracture Mechanics, 106(4), 102440. DOI 10.1016/j.tafmec.2019.102440. [Google Scholar] [CrossRef]

33. Molnar, G., Gravouil, A. (2017). 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture. Finite Elements in Analysis and Design, 130, 27–38. DOI 10.1016/j.finel.2017.03.002. [Google Scholar] [CrossRef]

34. Miehe, C., Hofacker, M., Welschinger, F. (2010). A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 199(45–48), 2765–2778. DOI 10.1016/j.cma.2010.04.011. [Google Scholar] [CrossRef]

35. Hesch, C., Weinberg, K. (2014). Thermodynamically consistent algorithms for a finite-deformation phase-field approach to fracture. International Journal for Numerical Methods in Engineering, 99(12), 906–924. DOI 10.1002/nme.4709. [Google Scholar] [CrossRef]

36. Hofacker, M., Miehe, C. (2012). Continuum phase field modeling of dynamic fracture: Variational principles and staggered Fe implementation. International Journal of Fracture, 178(1–2), 113–129. DOI 10.1007/s10704-012-9753-8. [Google Scholar] [CrossRef]

37. Seles, K., Lesicar, T., Tonkovic, Z., Soric, J. (2019). A residual control staggered solution scheme for the phase-field modeling of brittle fracture. Engineering Fracture Mechanics, 205, 370–386. DOI 10.1016/j.engfracmech.2018.09.027. [Google Scholar] [CrossRef]

38. Riks, E. (1979). Incremental approach to the solution of snapping and buckling problems. International Journal of Solids and Structures, 15(7), 529–551. DOI 10.1016/0020-7683(79)90081-7. [Google Scholar] [CrossRef]

39. Wriggers, P., Simo, J. C. (1990). A general procedure for the direct computation of turning and bifurcation points. International Journal for Numerical Methods in Engineering, 30(1), 155. DOI 10.1002/(ISSN)1097-0207. [Google Scholar] [CrossRef]

40. Deborst, R. (1987). Computation of postbifurcation and post-failure behavior of strain-softening solids. Computers & Structures, 25(2), 211–224. DOI 10.1016/0045-7949(87)90144-1. [Google Scholar] [CrossRef]

41. Gutierrez, M. A. (2004). Energy release control for numerical simulations of failure in quasi-brittle solids. Communications in Numerical Methods in Engineering, 20(1), 19–29. DOI 10.1002/cnm.649. [Google Scholar] [CrossRef]

42. May, S., Vignollet, J., de Borst, R. (2016). A new arc-length control method based on the rates of the internal and the dissipated energy. Engineering Computations, 33(1), 100–115. DOI 10.1108/EC-02-2015-0044. [Google Scholar] [CrossRef]

43. Barbieri, E., Ongaro, F., Pugno, N. M. (2017). A J-integral-based arc-length solver for brittle and ductile crack propagation in finite deformation-finite strain hyperelastic solids with an application to Graphene Kirigami. Computer Methods in Applied Mechanics and Engineering, 315(2), 713–743. DOI 10.1016/j.cma.2016.10.043. [Google Scholar] [CrossRef]

44. Wu, J. Y. (2018). Robust numerical implementation of non-standard phase-field damage models for failure in solids. Computer Methods in Applied Mechanics and Engineering, 340(14), 767–797. DOI 10.1016/j.cma.2018.06.007. [Google Scholar] [CrossRef]

45. Bibi, S., Abbas, M., Misro, M. Y., Hu, G. (2019). A novel approach of hybrid trigonometric Bezier curve to the modeling of symmetric revolutionary curves and symmetric rotation surfaces. IEEE Access, 7, 165779–165792. DOI 10.1109/ACCESS.2019.2953496. [Google Scholar] [CrossRef]

46. Usman, M., Abbas, M., Miura, K. T. (2020). Some engineering applications of new trigonometric cubic Bezier-like curves to free-form complex curve modeling. Journal of Advanced Mechanical Design Systems and Manufacturing, 14(4), JAMDSM0048. DOI 10.1299/jamdsm.2020jamdsm0048. [Google Scholar] [CrossRef]