A Damage-Based Framework for Flexible Cohesive Softening Laws in Abaqus without User Subroutines

1  Introduction

Cohesive zone models (CZMs) have become a standard tool for the numerical analysis of debonding, delamination, and fracture initiation and propagation in a wide range of quasi-brittle materials and structures. Unlike approaches directly rooted in classical fracture mechanics, CZMs are particularly well suited to problems in which a finite process zone develops ahead of the crack tip, allowing tractions to degrade progressively from a peak value to complete separation. This makes them especially attractive for problems involving material and geometric nonlinearities, three-dimensional configurations, and crack-initiation phenomena, which cannot be naturally addressed within traditional fracture-mechanics frameworks [1,2]. Since the seminal works of Dugdale [3] and Barenblatt [4], a large variety of cohesive formulations have been proposed, all relying on the introduction of an interface across which a displacement discontinuity is permitted. Such interfaces may be prescribed a priori based on physical considerations or introduced adaptively using partition-of-unity-based techniques [5–8]. At each interface point, a traction–separation law governs the relationship between the displacement jump and the corresponding tractions, with normal and shear components defining the possible fracture modes. An essential characteristic of any CZM is therefore its one-dimensional traction–separation response, which arises either under pure-mode loading or in mode-independent formulations where different modes do not exhibit explicit coupling.

In pure-mode and mode-independent settings, the traction–separation law is commonly characterized by an initial linear elastic response with stiffness Kp, followed by the attainment of a peak traction σc and a subsequent softening regime leading to complete decohesion. The total energy dissipated during this process, given by the area under the traction–separation curve, is equated to the fracture energy Gc. While the elastic stiffness is typically treated as a penalty parameter whose value can vary within a wide admissible range without significantly affecting the solution, the peak traction is a material-dependent quantity related to the size of the fracture process zone and the intrinsic strength of the interface. Its value is often difficult to measure experimentally and is usually inferred indirectly. From a numerical standpoint, large values of σc may also require increasingly fine discretizations in the vicinity of the process zone, leading to a substantial increase in computational cost [9,10]. As a consequence, the shape of the softening branch of the traction–separation law plays a critical role in controlling both the physical response and the numerical behavior of cohesive simulations.

With the widespread adoption of commercial finite element solvers such as Abaqus [11], the simulation of complex interfacial fracture problems has become increasingly accessible within native modeling environments. However, despite their general versatility, such solvers typically offer only a limited set of built-in cohesive laws. In Abaqus, cohesive elements are natively restricted to bilinear and linear–exponential traction–separation laws, which are often insufficient to capture the wide range of interfacial behaviors observed across different material systems. This limitation is particularly restrictive for bioinspired composites and architected interfaces, where softening behavior is known to deviate significantly from simple linear or exponential forms and is frequently modeled using custom finite element codes or user subroutines [12–15]. Although user subroutines provide considerable modeling freedom, their development requires advanced expertise, is platform dependent, and can significantly hinder rapid model development, testing, and parametric exploration [16–18].

Abaqus does, however, allow the damage evolution to be specified directly in tabular form as a function of the effective separation, a feature primarily intended for cases in which experimental damage data are available. In practice, obtaining reliable experimental traction–separation data is notoriously challenging and often infeasible [10,19]. In this work, we exploit this existing capability in a fundamentally different way by constructing the damage evolution law analytically from a prescribed traction–separation softening curve. By defining the damage variable directly from a desired softening response, any admissible analytical traction–separation law, or one obtained by fitting experimental data, can be implemented within the native Abaqus framework without the need for user subroutines. The proposed approach is formulated for pure-mode and mode-independent cohesive behavior, assuming similar interfacial properties—initial stiffness Kp, cohesive strength σc, and fracture energy Gc—in the normal and shear directions.

To demonstrate the generality of the proposed formulation, seven representative traction–separation laws are considered, including bilinear [20], linear–exponential [21], linear–cubic [22], bilinear–exponential [23], trapezoidal [24,25], a modified Park–Paulino–Roesler (PPR) law [26], and a general softening form capable of reproducing a broad range of cohesive responses. These laws encompass most cohesive behaviors commonly reported in the literature and allow a systematic assessment of the influence of softening shape under controlled conditions [19,27]. The proposed methodology is validated through a series of numerical examples, including mode-I and mode-II patch tests, a double-cantilever-beam benchmark, and the fracture analysis of a bioinspired nacre-like composite layer. Together, these examples demonstrate that complex and physically meaningful softening behaviors can be modeled accurately and efficiently in Abaqus without user subroutines, enabling rapid exploration and design of interfacial fracture models within a standard commercial finite element environment.

2  Method

2.1 Traction–Separation Relationship in Abaqus

The traction–separation (T–S) formulation implemented in Abaqus cohesive elements assumes an initially linear-elastic response followed by damage initiation and subsequent softening. At each integration point, Abaqus defines nominal tractions as force components divided by the original area, and nominal strains as separations divided by an original constitutive thickness t0 (with the transverse shear components evaluated using the standard Abaqus averaging procedure). Denoting the nominal traction vector by T=(Tn,Ts,Tt)T and the separation vector by Δ=(Δn,Δs,Δt)T, the nominal strains are

εn=Δnt0,εs=Δst0,εt=Δtt0(1)

The default value is t0=1, so that nominal strains become equal to relative separations (i.e., relative displacements of the cohesive faces). The constitutive thickness governing the T–S response is generally distinct from the geometric thickness of the interface; here we adopt t0=1 throughout.

In the undamaged regime, Abaqus relates tractions and separations through an elastic constitutive matrix,

T=KΔ,K=(knnknskntknsksskstkntkstktt)(2)

To isolate the effect of softening behavior in the following sections, we consider an uncoupled penalty formulation prior to damage, with knn=kss=ktt=Kp and all coupling terms set to zero (kns=knt=kst=0). Under mixed normal and shear separations, Abaqus introduces an effective separation measure [20],

Δm=⟨Δn⟩2+Δs2+Δt2(3)

where ⟨⋅⟩ denotes the Macaulay bracket, ensuring that compressive normal separation does not contribute to the driving measure. For Δm≤δmc (here, δmc denotes the effective separation at the onset of damage initiation), the response is elastic and the traction components reduce to

Tne=knnΔn,Tse=kssΔs,Tte=kttΔt(4)

Damage initiates when a prescribed failure criterion is satisfied. In this work, we adopt the quadratic nominal-stress criterion [28], with compressive normal tractions excluded from initiation,

(⟨Tn⟩τn)2+(Tsτs)2+(Ttτt)2=1(5)

where τn,τs,τt denote the critical tractions (interfacial strengths) in the normal and shear directions. Beyond initiation, for δmc<Δm<δmf (here, δmf denotes the effective separation at complete decohesion), the interface softens through a scalar damage variable D∈[0,1] that evolves monotonically from D=0 at the onset of damage to D=1 at complete failure. Abaqus degrades the elastic tractions according to

Tnd={(1−D)Tne,Tne≥0Tne,Tne<0Tsd=(1−D)TseTtd=(1−D)Tte(6)

where Tie are given by Eq. (4). For a pure-mode loading path with i∈{n,s,t}, the resulting one-dimensional T–S response can be written in the standard piecewise form,

Ti(Δi)={kiiΔi,0≤Δi≤δic(1−D)kiiΔi,δic<Δi<δif0,Δi≥δif(7)

with δic and δif denoting the mode-i separations at damage initiation and complete failure, respectively.

2.2 Controlling Softening Shape through Damage Variable

The post-peak softening branch of a traction–separation law (TSL) governs the evolution of the fracture process zone and can strongly influence both the predicted response and the numerical requirements of cohesive simulations. In Abaqus, however, the native cohesive formulation provides only linear and exponential softening options; other constitutive forms typically require user-defined elements implemented through subroutines. Abaqus does allow the damage evolution to be prescribed directly in tabular form as a function of the effective separation, a feature primarily intended for cases in which experimental T–S data are available. Because such data are difficult to obtain reliably, we instead construct the damage evolution analytically from a prescribed softening curve, enabling rapid implementation and systematic comparison of cohesive softening behaviors within the native Abaqus framework.

The proposed construction is formulated for a mode-independent cohesive response with similar interfacial properties in the normal and shear directions. Specifically, the elastic response is linear and uncoupled with a common penalty stiffness knn=kss=ktt=Kp and vanishing coupling terms; the interfacial strengths are equal, τn=τs=τt=σc; and the fracture energies are equal, Gnc=Gsc=Gtc=Gc, with damage driven by the effective separation and complete failure attained when the total work of separation equals Gc. Under these assumptions, the cohesive response can be represented by a single effective TSL in terms of the effective separation Δm defined in Eq. (3). The corresponding effective traction Tm is then written as

Tm(Δm)={KpΔm,0≤Δm≤δmc(1−D)KpΔm,δmc<Δm<δmf0,Δm≥δmf(8)

where δmc and δmf denote the effective separations at damage initiation and complete failure, respectively. Let Tmd(Δm) denote a prescribed post-initiation (damaged) softening curve over δmc<Δm<δmf. Combining Eq. (8) with the Abaqus damage degradation form yields a direct mapping from the desired softening law to the damage variable,

D(Δm)={0,0≤Δm≤δmc1−Tmd(Δm)KpΔm,δmc<Δm<δmf1,Δm≥δmf(9)

It is noted that the prescribed softening law must be admissible under this mapping, in the sense that the resulting damage evolution remains bounded and physically meaningful, with 0≤D(Δm)≤1 and monotonic non-decreasing behavior over the softening regime. Equivalently, in δmc<Δm<δmf the damaged traction must satisfy 0≤Tmd(Δm)≤KpΔm, which prevents negative damage and ensures that traction degradation is consistent with the Abaqus constitutive structure.

The fracture energy, defined as the area under the effective T–S curve, naturally decomposes into the elastic contribution prior to damage initiation and the dissipated contribution during softening. Defining Ge, Gd, and Gc as the pre-initiation, post-initiation, and total fracture energies, respectively, one obtains

Ge=∫0δmcTme(Δ)dΔ=12Kp(δmc)2Gd=∫δmcδmfTmd(Δ)dΔGc=Ge+Gd=12Kp(δmc)2+∫δmcδmfTmd(Δ)dΔ(10)

For a prescribed (Kp,δmc,Gc) and a chosen softening curve Tmd(Δm), the complete-failure separation δmf follows from Eq. (10).

Building on the above formulations, specific softening behaviors can be prescribed by selecting appropriate analytical expressions for the damaged traction Tmd(Δm). As a first step, we consider the linear and exponential softening laws that are natively available in Abaqus, which serve both as reference cases and as a direct verification of the proposed damage-based construction. For linear softening, the damaged traction is defined as

Tmd(Δm)={σcδmf−Δmδmf−δmc,δmc<Δm<δmf0,Δm≥δmf(11)

which corresponds to a constant degradation rate from the peak traction σc to complete failure. The resulting evolution of both traction and damage variable as functions of the effective separation is shown in Fig. 1a. An exponential softening response can be defined in an analogous manner to

Figure 1: Evolution of traction and damage variable as functions of the effective separation for reference softening laws. (a) Linear softening law. (b) Exponential softening law for several values of the shape parameter αr, illustrating its influence on the post-peak degradation rate.

Tmd(Δm)={σc[1−1−exp(−αr(Δm−δmc)/(δmf−δmc))1−exp⁡(−αr)],δmc<Δm<δmf0,Δm≥δmf(12)

where αr is a dimensionless shape parameter controlling the rate of traction decay. Smaller values of αr lead to a more gradual reduction in traction, whereas larger values produce a sharper post-peak drop. Fig. 1b illustrates the corresponding traction–separation and damage evolution curves for several representative choices of αr.

For both cases, the damage variable D(Δm) follows directly Eq. (9) once Kp, σc, and Gc are specified. Because Abaqus natively supports linear and exponential traction–separation laws, these two softening forms provide a convenient and rigorous validation of the proposed formulation. Mode-I and mode-II patch tests, presented in Section 3.1, demonstrate that cohesive responses obtained using the analytically generated damage evolution are indistinguishable from those produced by the default Abaqus implementations when identical material parameters are employed.

Having established equivalence with the linear and exponential softening laws available in Abaqus, the same construction can be applied to a wide range of cohesive laws reported in the literature, enabling their direct implementation through analytically generated damage evolution. As a representative example, cubic softening has been proposed to model damage evolution in composite structures subjected to multiaxial loading, where a smooth cubic decay captures gradual stiffness degradation after peak traction [22]. In this case, the damaged traction can be written as

Tmd(Δm)={σc[1−3r2+2r3],δmc<Δm<δmf0,Δm≥δmf(13)

where r=(Δm−δmc)/(δmf−δmc) and produces a smooth transition to zero traction with vanishing slope at complete failure. The corresponding traction and damage evolution are shown in Fig. 2a.

Figure 2: Representative non-native traction–separation laws implemented through analytically generated damage evolution. (a) Cubic softening. (b) Linear–exponential softening with a transition at δmb. (c) Trapezoidal TSL with bilinear softening with a finite traction plateau. (d) Modified PPR softening for several values of the shape parameter α, illustrating the transition between convex and concave post-peak responses.

To account for fiber-bridging effects during mode-I delamination of composite laminates, particularly under environmental degradation, a linear–exponential softening law has also been proposed [23]. In this formulation, an initial linear degradation from σc to a bridging stress σb is followed by an exponential decay, giving

Tmd(Δm)={σc−(σc−σb)Δm−δmcδmb−δmc,δmc<Δm≤δmbσbexp[−γΔm−δmbδmf−δmb],δmb<Δm<δmf0,Δm≥δmf(14)

where δmb denotes the transition separation and γ controls the exponential decay rate. Fig. 2b illustrates this response for δmb=12δmf, σb=711σc, and γ=5.

TSLs with trapezoidal form have been widely employed to model crack initiation and growth resistance in elastic–plastic materials, where a finite plateau in traction represents stable crack growth prior to final separation [24,25]. The corresponding damaged traction is defined as

Tmd(Δm)={σc,δmc<Δm≤δmsσcδmf−Δmδmf−δms,δms<Δm<δmf0,Δm≥δmf(15)

where δms defines the extent of the traction plateau. The resulting traction and damage evolution are shown in Fig. 2c.

Finally, to capture a broad spectrum of cohesive responses within a unified framework, we consider a modified form of the PPR potential-based model [26], specialized here to a mode-independent setting. The damaged traction is expressed as

Tmd(Δm)={σcA[m(1−Δmδmf)α(mα+Δmδmf)m−1−α(1−Δmδmf)α−1(mα+Δmδmf)m],δmc<Δm<δmf0,Δm≥δmf(16)

where α controls the convexity or concavity of the softening response and the normalization constant A is given by

A=m(1−λ)α(mα+λ)m−1−α(1−λ)α−1(mα+λ)mm=α(α−1)λ21−αλ2,λ=δmcδmf(17)

Fig. 2d illustrates the influence of α on the traction–separation response and the corresponding damage evolution.

While the preceding TSLs cover a range of cohesive behaviors reported in the literature, their functional forms are typically tailored to specific mechanisms or applications. To provide a compact yet flexible representation that can span a broad family of softening responses within a single definition, we propose the following general mode-independent softening form:

Tmd(Δm)={σc,δmc<Δm≤δmsσc(δmf−Δmδmf−δms)p,δms<Δm<δmf0,Δm≥δmf(18)

where δms defines the extent of a traction plateau and p controls the curvature of the post-plateau decay. For δms=δmc the plateau vanishes and Eq. (18) reduces to a one-parameter family of direct softening laws, including the linear (bilinear) law when p=1 and an elastic-perfectly plastic response when p=0. For δmc<δms<δmf, the initial plateau represents a stable separation phase (commonly known as shear flow) prior to softening, and the decay exponent p governs the subsequent degradation rate: p=1 produces a trapezoidal law, p>1 yields a more abrupt (convex) post-peak drop, and 0<p<1 produces a more gradual (concave) decay. Fig. 3 illustrates representative traction–separation and damage evolution curves for this general form and highlights its ability to reproduce several of the reference shapes considered above through appropriate parameter choices. The influence of these shape controls on the macroscopic response and fracture patterns is demonstrated in Section 3.3 for a bioinspired nacre-like composite layer.

Figure 3: General softening law in Eq. (18) and the corresponding damage evolution obtained from Eq. (9). The parameter δms controls the length of the traction plateau or shear flow, while the exponent p governs the curvature of the post-plateau decay, enabling linear (p=1), trapezoidal (p=1 with δms>δmc), concave (0<p<1), convex (p>1), and perfectly plastic (p=0) limiting behaviors.

The same framework also accommodates softening laws defined directly from experimental data. If a T–S response can be obtained or fitted as a discrete function Tmd(Δm) over δmc≤Δm≤δmf, the corresponding damage evolution immediately follows Eq. (9). The resulting law can then be implemented in Abaqus by specifying a tabular damage evolution in terms of the relative effective displacement (Δm−δmc), with the two-column data format [(Di,Δm,i−δmc),…].

3  Examples

To assess the proposed damage-based construction of cohesive softening laws in Section 2.2, three numerical examples are considered: (i) mode-I and mode-II patch tests used to verify the implementation against Abaqus-native cohesive laws under loading and unloading/reloading, (ii) a double cantilever beam (DCB) benchmark compared against published reference results [29], and (iii) a tensile fracture simulation of a bioinspired nacre-like composite layer illustrating the sensitivity of the macroscopic response and crack patterns to the chosen traction–separation softening behavior.

3.1 Patch Tests: Mode-I and Mode-II

We first verify the formulation using two patch tests designed to generate a uniform stress state in the bulk while localizing separation across a single cohesive interface. The geometries and boundary conditions are shown in Fig. 4. In both cases, a 1.0mm×1.0mm square plate is modeled with E=200GPa and ν=0.3, and a single two-dimensional cohesive element (COH2D4) is inserted either along the mid-section for mode-I (Fig. 4a) or along the diagonal for mode-II (Fig. 4b). The cohesive parameters are taken as Kp=4.0×103, σc=8.0MPa, and Gc=0.1kJ/m2. Since Abaqus provides linear and exponential softening natively, we use these two cases as references and compare: (i) the analytically generated damage evolution corresponding to linear softening in Eq. (11) and (ii) the exponential case in Eq. (12) with αr=2, against the default Abaqus cohesive implementation under identical parameters. Unloading and reloading segments are included in the prescribed displacement histories to confirm that the alternative construction reproduces the native linear unloading/reloading response.

Figure 4: Patch-test configurations and boundary conditions: (a) mode-I opening with a cohesive interface along the mid-section and (b) mode-II shearing with a cohesive interface along the diagonal.

For the mode-I test, the bulk is discretized with two four-node plane strain elements (CPE4) separated by a 2D cohesive element (COH2D4) placed along Y=0.5mm. With the bottom edge fixed, the top edge is displaced upward to 0.008mm, then reversed to −0.001mm, and finally re-extended to drive complete failure. The resulting traction–separation responses for linear and exponential softening are shown in Fig. 5a,b. As loading proceeds, the cohesive traction increases linearly to the peak σc and subsequently follows the prescribed softening law; during unloading and reloading, the response follows the native linear path. The curves obtained using the analytically generated damage evolution are indistinguishable from those produced by the Abaqus default cohesive law for both softening types.

Figure 5: Verification of the damage-based construction against the Abaqus-native cohesive implementation. Mode-I patch test: (a) linear and (b) exponential softening with αr=2. Mode-II patch test: (c) linear and (d) exponential softening with αr=2. Arrows indicate unloading/reloading paths; the damage-based and Abaqus-native responses are indistinguishable in all cases.

For the mode-II test, the bulk is discretized with two three-node plane strain elements (CPE3) and the cohesive element is placed along the diagonal. With the left and bottom edges fixed, displacements are imposed on the top and right edges to create a shear-dominated separation history, including reversal to enforce unloading and reloading prior to complete failure. The corresponding traction–separation curves are shown in Fig. 5c,d, where the peak traction reaches σc before softening and the unloading/reloading segments follow the expected linear paths. In both mode-I and mode-II patch tests, the damage-based construction reproduces the Abaqus-native bilinear and linear–exponential cohesive responses exactly, confirming the correctness of the implementation and providing a validated foundation for the non-native softening laws introduced in Section 2.2.

3.2 DCB Analysis

The influence of cohesive softening shape on a structural-scale response is examined through a mode-I double cantilever beam (DCB) delamination benchmark. Following the reference study in [29], the specimen is a unidirectional T300/1076 graphite/epoxy laminate of length 150mm and width 25mm, consisting of two 1.5mm-thick arms joined by a 0.01mm-thick cohesive layer, with an initial delamination length a0=30.5mm. The orthotropic elastic properties of the arms and the cohesive parameters used for all simulations are reported in Table 1.

The finite element model comprises two layers of four-node plane strain elements (CPE4) connected by four-node cohesive elements (COH2D4). A global mesh size of 0.25mm is adopted with three elements through each arm thickness direction, consistent with common DCB modeling recommendations [29,30]. For the given material parameters and initial crack length a0=30.5mm, the referenced benchmark studies report that a global mesh size of 0.5mm is sufficient to achieve mesh convergence; the present model therefore employs a refined mesh of 0.25mm to further minimize discretization effects. The same discretization is maintained for all softening laws to ensure consistent cohesive-zone resolution and enable fair comparison of crack-growth behavior without mesh-dependent bias. As shown in Fig. 6a, the lower-arm corner node is pinned and a vertical displacement up to δ=8mm is applied at the upper-arm corner node; the resulting deformed configuration and crack propagation are illustrated in Fig. 6b. To isolate the role of the post-peak softening shape, four softening laws are considered while keeping the cohesive properties identical (Table 1): (i) linear softening (Eq. (11)), (ii) exponential softening with αr=2 (Eq. (12)), (iii) the proposed general form with a zero plateau (δms=δmc) and p=2 (Eq. (18)), and (iv) the modified PPR form with α=2 (Eq. (16)). Fig. 6c compares the predicted load–displacement curves against the benchmark data of [29]. The benchmark reports a critical load Pcrit=62.0N at a critical displacement δcrit=1.52mm. All four softening laws reproduce the global response closely, yet the modified PPR case provides the best overall agreement, with R2=0.984, predicting Pcrit=62.42N at δcrit=1.89mm. This comparison highlights the practical value of the proposed damage-based construction: cohesive laws such as the modified PPR response, which are not available as native options in Abaqus and are commonly introduced via user subroutines [26,31], can be implemented and assessed directly within the standard Abaqus workflow while preserving consistency in (Kp,σc,Gc).

Figure 6: DCB benchmark setup and sensitivity to cohesive softening shape. (a) Boundary conditions and loading configuration. (b) Deformed configuration illustrating delamination growth. (c) Load–displacement responses for different traction–separation softening laws compared with benchmark data from [29].

3.3 Nacre-Like Composite Analysis

The mechanical response and failure patterns of bioinspired organic–inorganic composites such as nacre are governed primarily by the properties of the organic interfaces rather than by the stiffness of the mineral tablets themselves [12,16]. CZM has therefore become a natural framework for studying deformation and fracture in such systems [13,17]. However, the wide diversity of interfacial behaviors observed in bioinspired materials cannot be adequately represented by bilinear or linear-exponential TSLs alone. As a result, previous studies have often relied on custom finite element codes or user subroutines to implement tailored cohesive laws [14,15,26,31], which limits flexibility and hinders rapid model development. Here, we demonstrate how the proposed general softening law in Eq. (18) enables systematic exploration of diverse interfacial behaviors within the native Abaqus environment.

A single layer of a nacre-like composite is considered in the form of a rectangular domain of 8400nm×8400nm. The mineral phase is generated as a two-dimensional polygonal array using centroidal Voronoi tessellation (CVT) technique, producing a statistically uniform yet disordered tablet arrangement representative of natural nacre; details of the generation procedure are provided in [32]. The two-dimensional microstructure is then extruded to a thickness of 300nm to form a three-dimensional mineral–organic layer, and zero-thickness cohesive elements are inserted systematically along all tablet boundaries to represent the organic matrix. The tablets are discretized with approximately 2300 linear hexahedral elements (C3D8), while about 2200 three-dimensional cohesive elements (COH3D8) are used along the interfaces. The element size along the interfaces is selected following established cohesive-zone modeling guidelines [33], where the characteristic cohesive-zone length may be estimated as lcz≈EGc/σc2. For the adopted material parameters, the mesh is chosen such that the cohesive element size satisfies h≲lcz/5, ensuring adequate resolution of the fracture process zone. We choose the minimal bound to avoid unwanted errors due to mesh resolution, and maintain the same discretization for all softening laws to enable consistent comparison of their structural effects without introducing mesh-dependent bias. The mineral tablets are modeled as a stiff elastic phase with Young’s modulus E=100GPa, Poisson’s ratio ν=0.28, and density ρ=3.0g/cm3, consistent with calcium carbonate and related biominerals. The cohesive interfaces are assigned with Kp=3.694×105, σc=50.0MPa, and Gc=0.5kJ/m2 [34]. Quasi-static uniaxial tension is applied by prescribing vertical displacement on the top surface (Y=8400nm) while fixing the bottom surface (Y=0), up to a nominal tensile strain of 0.6%, where all other surfaces remain traction free. The Abaqus/Explicit solver is employed to robustly capture progressive interface softening, contact interactions, and crack evolution.

To examine the influence of softening shape, four representative combinations of the plateau parameter δms and the decay exponent p are considered within the general law in Eq. (18): linear softening without a plateau (δms=δmc, p=1.0), elastic–perfectly plastic behavior (δms=δmc, p=0.0), shear flow followed by a concave decay (δms=0.005, p=0.5), and shear flow followed by a convex decay (δms=0.005, p=2.0). Although the geometry and material parameters are identical in all cases, Fig. 7a–d shows that both the stress distribution and the resulting crack paths differ markedly depending on the chosen softening behavior, as different post-peak traction reductions alter the balance between stress redistribution and local energy dissipation along competing interfaces, thereby favoring either distributed interfacial separation across multiple tablet boundaries or the early localization of fracture along a dominant crack path. These differences become even more pronounced at the macroscopic level: the bulk stress–strain curves in Fig. 7e closely mirror the underlying interfacial softening characteristics, demonstrating that the global response of the composite is strongly controlled by the post-peak traction–separation behavior through its direct influence on load transfer, progressive decohesion, and the accumulation of inelastic separation across the tablet interfaces. This sensitivity is consistent with prior observations in bioinspired composite modeling [12,14,15] and clearly illustrates the necessity of flexible softening descriptions. Collectively, these results confirm that the proposed general softening law enables systematic, physically meaningful exploration of interfacial behavior in complex composite systems without recourse to user subroutines.

Figure 7: Fracture response of a nacre-like composite layer governed by the general softening law in Eq. (18). Final stress distribution (σ22) and crack patterns at 0.6% applied tensile strain for (a) linear softening without a plateau (δms=δmc, p=1.0), (b) elastic–perfectly plastic behavior (δms=δmc, p=0.0), (c) shear flow followed by a concave decay (δms=0.005, p=0.5), and (d) shear flow followed by a convex decay (δms=0.005, p=2.0). (e) corresponding macroscopic stress–strain responses, highlighting the strong dependence of the bulk behavior on the interfacial softening shape (Fig. 3 illustrates the corresponding softening shapes).

4  Limitations and Future Directions

In Sections 2 and 3, we present the proposed framework strictly formulated within the native cohesive constitutive structure of Abaqus, where a scalar damage variable degrades the nominal traction response through tabular damage evolution defined with respect to an effective separation measure [11]. Accordingly, the range of implementable traction–separation behaviors is constrained primarily by the underlying constitutive formulation of cohesive elements in Abaqus rather than by the proposed mapping itself. In particular, a mode-independent cohesive response with identical normal and shear parameters is adopted under an effective separation definition. Though this assumption provides limited scope for applying the proposed softening formulations, it reflects common practice in cohesive-zone modeling when experimental validation is unavailable to justify different interfacial behavior and properties in normal and shear directions and reliable mixed-mode calibration data are limited to construct a compact and numerically robust representation of constrained coupling functions [9,20,35,36]. It is nonetheless recognized that many interfaces exhibit pronounced mode dependence and evolving mode mixity, and that mixed-mode cohesive predictions are sensitive to the selected interaction law and experimental calibration [30,33,37,38]. Within Abaqus, general mixed-mode coupling typically requires additional constitutive definitions—such as mode-mix-dependent evolution or distinct normal/shear degradation laws—beyond a single prescribed softening curve [11]. Constructing a fully consistent damage surface compatible with displacement-type tabular evolution is therefore nontrivial and not supported by an equally direct solver-native workflow. The proposed method should thus be interpreted as a rigorous and practical approach for implementing a broad family of cohesive softening shapes within the scalar-damage/effective-separation framework of Abaqus, while fully general mixed-mode coupling remains a natural extension that may require additional constitutive assumptions or user element/subroutine implementations.

Additionally, the present study is limited to rate-independent behavior to isolate the influence of softening-shape flexibility and preserve compatibility with standard quasi-static cohesive settings. Rate effects can significantly influence interfacial strength and fracture energy in many material systems [39], and incorporating general rate dependence typically introduces additional internal variables that are commonly implemented via user subroutines. Future extensions may therefore consider rate-dependent or visco-cohesive formulations while maintaining the present mapping philosophy wherever compatible solver input structures exist.

Finally, in the present study, we adopt the cohesive parameters from the literature for verification and benchmarking, while isolating the effect of softening-shape flexibility and demonstrating that non-native laws can be realized without UEL/VUEL implementations. For predictive applications, systematic calibration of (Kp,σc,Gc) and additional shape parameters is essential and is increasingly addressed through inverse finite-element and optimization frameworks [40–42]. Because the proposed approach regenerates tabular damage evolution directly from a prescribed traction–separation law, it is naturally compatible with such inverse pipelines, enabling iterative parameter identification without altering the underlying element formulation. These considerations outline the intended applicability of the solver-native framework while identifying clear directions for future developments in mixed-mode coupling, rate dependence, and systematic parameter calibration.

5  Conclusion

This work presents a damage-based construction that enables flexible and physically consistent control of cohesive softening behavior within the native Abaqus framework, without the use of user subroutines. By defining the damage evolution analytically from a prescribed post-initiation traction–separation response, the approach allows a broad class of admissible softening laws to be implemented for pure-mode and mode-independent cohesive behavior, yielding similar effective stiffness, strength, and fracture energy values in the normal and shear directions. Mode-I and mode-II patch tests confirm that the proposed definition reproduces the Abaqus-native bilinear and linear–exponential cohesive responses exactly, including unloading and reloading behavior. A double cantilever beam benchmark further illustrates how variations in softening shape influence the structural response and demonstrates that non-native laws, such as modified PPR softening, can be assessed accurately within the standard Abaqus workflow. Finally, simulations of a bioinspired nacre-like composite layer show that the proposed general softening form provides compact yet powerful control over interfacial behavior, directly governing macroscopic stress–strain response and fracture patterns in complex microstructures. These results are obtained within the scalar-damage and effective-separation constitutive definitions of cohesive elements in Abaqus and are therefore limited to rate-independent, mode-independent formulations with identical normal and shear parameters, as assumed in this study. Within this scope, the proposed framework establishes an accessible and robust pathway for implementing, comparing, and exploring cohesive softening behaviors in Abaqus, providing a practical solver-native workflow for parametric and sensitivity studies under the stated assumptions, without the need for user subroutines.

Acknowledgement: Not applicable.

Funding Statement: A grant from the University of Texas at San Antonio, Office of the Vice President for Research, funded this research.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Md Jalal Uddin Rumi, Xiaowei Zeng; methodology, Md Jalal Uddin Rumi; software, Md Jalal Uddin Rumi; validation, Md Jalal Uddin Rumi; formal analysis, Md Jalal Uddin Rumi; investigation, Md Jalal Uddin Rumi; resources, Md Jalal Uddin Rumi, Xiaowei Zeng; data curation, Md Jalal Uddin Rumi; writing—original draft preparation, Md Jalal Uddin Rumi; writing—review and editing, Md Jalal Uddin Rumi, Xiaowei Zeng; visualization, Md Jalal Uddin Rumi; supervision, Xiaowei Zeng; project administration, Xiaowei Zeng; funding acquisition, Xiaowei Zeng. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: All data generated or analyzed during this study are included in this article. For the detailed implementation of Voronoi modeling of the nacre-like composite layer described in Section 3.3, please see the following repository: https://github.com/Rumi381/PolyGen.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Wells GN, Sluys L. A new method for modelling cohesive cracks using finite elements. Int J Numer Methods Eng. 2001;50(12):2667–82. doi:10.1002/nme.143. [Google Scholar] [CrossRef]

2. Wells G, De Borst R, Sluys L. A consistent geometrically non-linear approach for delamination. Int J Numer Methods Eng. 2002;54(9):1333–55. doi:10.1002/nme.462. [Google Scholar] [CrossRef]

3. Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids. 1960;8(2):100–4. doi:10.1016/0022-5096(60)90013-2. [Google Scholar] [CrossRef]

4. Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech. 1962;7(6):55–129. doi:10.1016/s0065-2156(08)70121-2. [Google Scholar] [CrossRef]

5. Möes N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech. 2002;69(7):813–33. doi:10.1016/s0013-7944(01)00128-x. [Google Scholar] [CrossRef]

6. Remmers JJC, Borst Rd, Needleman A. A cohesive segments method for the simulation of crack growth. Comput Mech. 2003;31(1):69–77. doi:10.1007/s00466-002-0394-z. [Google Scholar] [CrossRef]

7. Melenk J, Babuška I. Approximation with harmonic and generalized harmonic polynomials in the partition of unity method. Comput Assist Methods Eng Sci. 2023;4(3–4):607–32. [Google Scholar]

8. Vellwock AE, Libonati F. XFEM for composites, biological, and bioinspired materials: a review. Materials. 2024;17(3):745. doi:10.3390/ma17030745. [Google Scholar] [PubMed] [CrossRef]

9. Alfano G, Crisfield M. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. Int J Numer Methods Eng. 2001;50(7):1701–36. doi:10.1002/nme.93. [Google Scholar] [CrossRef]

10. Kim HG, Howe R, Wiebe R, Spottswood SM, O’Hara PJ, Salviato M. Experimental characterization of cohesive laws for mode-II interlaminar fracture in geometrically scaled composites using through-thickness deformation analysis. Eng Fract Mech. 2025;326:111361. doi:10.2139/ssrn.5227657. [Google Scholar] [CrossRef]

11. Smith M. ABAQUS/standard user’s manual, version 6.9. Providence, RI, USA: Dassault Systèmes Simulia Corp; 2009. [Google Scholar]

12. Yang X, Rumi MJU, Zeng X. Computational investigation of the mechanical response of a bioinspired nacre-like nanocomposite under three-point bending. J Compos Sci. 2024;8(5):173. doi:10.3390/jcs8050173. [Google Scholar] [CrossRef]

13. Lin L, Samuel J, Zeng X, Wang X. Contribution of extrafibrillar matrix to the mechanical behavior of bone using a novel cohesive finite element model. J Mech Behav Biomed Mater. 2017;65:224–35. doi:10.1016/j.jmbbm.2016.08.027. [Google Scholar] [PubMed] [CrossRef]

14. Lin L, Wang X, Zeng X. Computational modeling of interfacial behaviors in nanocomposite materials. Int J Solids Struct. 2017;115:43–52. doi:10.1016/j.ijsolstr.2017.02.029. [Google Scholar] [PubMed] [CrossRef]

15. Lin L, Wang X, Zeng X. An improved interfacial bonding model for material interface modeling. Eng Fract Mech. 2017;169:276–91. doi:10.1016/j.engfracmech.2016.10.015. [Google Scholar] [PubMed] [CrossRef]

16. Yang X, Maghsoudi-Ganjeh M, Zeng X. Computational investigation of the mechanical behavior of a Bone-inspired Nanocomposite material. J Compos Sci. 2023;7(8):341. doi:10.3390/jcs7080341. [Google Scholar] [CrossRef]

17. Maghsoudi-Ganjeh M, Lin L, Wang X, Zeng X. Computational investigation of ultrastructural behavior of bone using a cohesive finite element approach. Biomech Model Mechanobiol. 2019;18(2):463–78. doi: 10.1007/s10237-018-1096-6. [Google Scholar] [PubMed] [CrossRef]

18. Sharma H, Singh A. Numerical implementation of a modified cohesive zone model for HCF behavior of adhesively bonded composite laminates under mixed mode loading. Int J Fatigue. 2024;181:108128. doi: 10.1016/j.ijfatigue.2023.108128. [Google Scholar] [CrossRef]

19. Wei C, Zhang J, Liechti KM, Wu C. Data driven modeling of interfacial traction-separation relations using a thermodynamically consistent neural network. Comput Methods Appl Mech Eng. 2023;404:115826. doi: 10.1016/j.cma.2022.115826. [Google Scholar] [CrossRef]

20. Camanho PP, Dávila CG. Mixed-mode decohesion finite elements for the simulation of delamination in composite materials. Washington, DC, USA: NASA; 2002. Report No.: NASA/TM-2002-211737. [Google Scholar]

21. Peerlings RH, de Borst R, Brekelmans W, Geers MG. Gradient-enhanced damage modelling of concrete fracture. Mech Cohesive Frict Mater Int J Exp Model Comput Mater Struct. 1998;3(4):323–42. [Google Scholar]

22. Donadon MV, De Almeida SFM, Arbelo MA, de Faria AR. A three-dimensional ply failure model for composite structures. Int J Aerosp Eng. 2009;2009(1):486063. [Google Scholar]

23. Wong KJ, Chong WWF, Goh KL. Fiber bridging mechanism in moisture-induced mode I delamination in carbon/epoxy composites: finite element analysis and experimental investigation. Polym Compos. 2023;44(2):1392–407. doi:10.1002/pc.27179. [Google Scholar] [CrossRef]

24. Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. J Mech Phys Solids. 1992;40(6):1377–97. doi:10.1016/0022-5096(92)90020-3. [Google Scholar] [CrossRef]

25. Tvergaard V, Hutchinson JW. The influence of plasticity on mixed mode interface toughness. J Mech Phys Solids. 1993;41(6):1119–35. doi:10.1016/0022-5096(93)90057-m. [Google Scholar] [CrossRef]

26. Park K, Paulino GH, Roesler JR. A unified potential-based cohesive model of mixed-mode fracture. J Mech Phys Solids. 2009;57(6):891–908. doi:10.1016/j.jmps.2008.10.003. [Google Scholar] [CrossRef]

27. Tao C, Zhang C, Ji H, Qiu J. A Paris-law-informed neural fatigue cohesive model and its application to open-hole composite laminates. Int J Solids Struct. 2023;267:112158. [Google Scholar]

28. Ye L. Role of matrix resin in delamination onset and growth in composite laminates. Compos Sci Technol. 1988;33(4):257–77. doi:10.1016/0266-3538(88)90043-7. [Google Scholar] [CrossRef]

29. Krueger R. A summary of benchmark examples to assess the performance of quasi-static delamination propagation prediction capabilities in finite element codes. J Compos Mater. 2015;49(26):3297–316. doi:10.1177/0021998314561812. [Google Scholar] [CrossRef]

30. Turon A, Davila CG, Camanho PP, Costa J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng Fract Mech. 2007;74(10):1665–82. doi:10.1016/j.engfracmech.2006.08.025. [Google Scholar] [CrossRef]

31. Park K, Paulino GH. Computational implementation of the PPR potential-based cohesive model in ABAQUS: educational perspective. Eng Fract Mech. 2012;93(3):239–62. doi:10.1016/j.engfracmech.2012.02.007. [Google Scholar] [CrossRef]

32. Rumi MJU, Han HC, Ye J, Feldman M, Gruslova A, Nolen D, et al. Polygen: an efficient framework for polycrystals generation and cohesive zone modeling in arbitrary domains. Eng Comput. 2025;41(5):3651–71. doi:10.1007/s00366-025-02180-6. [Google Scholar] [CrossRef]

33. Song K, Dávila CG, Rose CA. Guidelines and parameter selection for the simulation of progressive delamination. In: 2008 ABAQUS User’s Conference. Providence, RI, USA: Dassault Systèmes Simulia Corp; 2008. p. 1–15. [Google Scholar]

34. Rumi MJU, Zeng X. Computational investigation of the interplay between interfacial and bulk material properties of bioinspired nacre-like composites. JOM. 2025;77(12):9226–40. doi:10.1007/s11837-025-07754-9. [Google Scholar] [CrossRef]

35. Park K, Paulino GH. Cohesive zone models: a critical review of traction-separation relationships across fracture surfaces. Appl Mech Rev. 2011;64(6):060802. doi:10.1115/1.4023110. [Google Scholar] [CrossRef]

36. Wciślik W, Pała T. Selected aspects of cohesive zone modeling in fracture mechanics. Metals. 2021;11(2):302. doi:10.3390/met11020302. [Google Scholar] [CrossRef]

37. Benzeggagh M, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos Sci Technol. 1996;56(4):439–49. doi:10.1016/0266-3538(96)00005-x. [Google Scholar] [CrossRef]

38. Nairn JA, Aimene YE. A re-evaluation of mixed-mode cohesive zone modeling based on strength concepts instead of traction laws. Eng Fract Mech. 2021;248(4):107704. doi:10.1016/j.engfracmech.2021.107704. [Google Scholar] [CrossRef]

39. Zhang ZJ, Paulino GH. Cohesive zone modeling of dynamic failure in homogeneous and functionally graded materials. Int J Plast. 2005;21(6):1195–254. [Google Scholar]

40. Shi T, Pang M, Wang Y, Zhang Y. Inverse parameter identification framework for cohesive zone models based on multi-island genetic algorithm. Eng Fract Mech. 2024;300:110005. doi:10.1016/j.engfracmech.2024.110005. [Google Scholar] [CrossRef]

41. Zhang W, Tang M, Mu H, Yang X, Zeng X, Tuo R, et al. Inverse design of nonlinear mechanics of bio-inspired materials through interface engineering and Bayesian optimization. Extrem Mech Lett. 2025;78(2):102359. doi:10.1016/j.eml.2025.102359. [Google Scholar] [CrossRef]

42. Gao Y, Xiao L, Li Y, Duan K, He Y, Li L. A deep learning-assisted inverse design framework for nacre-like composites with excellent specific damping performance. Compos Part A Appl Sci Manuf. 2025;198(37):109050. doi:10.1016/j.compositesa.2025.109050. [Google Scholar] [CrossRef]