Study on the cohesive edge crack in a square plate with the cohesive element method

The size of the fully developed process zone (FDPZ) is needed for the arrangement of displacement sensors in fracture experiments and choosing element size in numerical models using the cohesive element method (CEM). However, the FDPZ size is generally not known beforehand. Analytical solutions for the exact FDPZ size only exist for highly idealised bodies, e.g. semi-infinite plates. With respect to fracture testing, the CEM is also a potential tool to extrapolate laboratory test results to full-scale while considering the size effect. A numerical CEM-based model is built to compute the FDPZ size for an edge crack in a finite square plate of different lengths spanning several magnitudes. It is validated against existing analytical solutions. After successful validation, the FDPZ size of finite plates is calculated with the same numerical scheme. The (FDPZ) size for finite plates is influenced by the cracked plate size and physical crack length. Maximum cohesive zone sizes are given for rectangular and linear softening. Further, for this setup, the CEM-based numerical model captures the size effect and can be used to extrapolate small-scale test results to full-scale.

Dirac delta function λ π E U cr /2 0 , characteristic length λ L G c E / 2 0 , characteristic length to normalise the plate size U Green's function for the edge crack COD ν Poisson's ratio Arbitrary crack face pressure on 0 ≤ X ≤ A, Y = 0 ± In cohesive zone theory, the fully developed process zone (FDPZ) is the region in which material degradation is activated. The size of this zone is needed in experiments for measuring purposes, e.g. the arrangement of displacement sensors in ice fracture experiments (Dempsey et al. 1999a(Dempsey et al. , b, 2004Lu et al. 2015aLu et al. , 2019, and in numerical models based on the Cohesive Element Method (CEM), where a size estimate is required for appropriate discretisation. However, the FDPZ size is generally not known beforehand. There are formulas for rough estimates (Turon et al. 2007), but it is nontrivial to obtain its exact length. Analytical solutions only exist for highly idealised cracked bodies, e.g. an edge crack in a semi-infinite plate. A closely related problem is the scaling of laboratory fracture tests to full-scale considering the size effect, i.e. the nominal strength depending on structure size, which can theoretically be done with the CEM (Bažant and Yavari 2005).
To tackle these issues, a CEM-based model for an edge crack in finite square plates of different sizes is built, where several magnitudes lie between the smallest and the largest plate. FDPZ sizes and crack-opening loads are computed and the ability of the CEM to capture the size effect is investigated.
To be more specific, the CEM is popular in many fields concerning material or structural failure, e.g. (Miller et al. 1999;Xu and Needleman 1994;Wang et al. 2019a;Feng et al. 2016;Konuk et al. 2009a, b;Konuk and Yu 2010;Pandolfi et al. 2000), due to its straightforward incorporation in FEM models. It is implemented in commercial software, such as Abaqus and LS-Dyna, to simulate fracture problems e.g. crack propagation and fragmentation. Here, it is used for its ability to resolve the Fracture Process Zone (FPZ). Yet, for CEM-based models, it is not straightforward to choose a consistent set of parameters, e.g. the maximum cohesive traction, due to their mutual influence (Turon et al. 2007;Blal et al. 2013). The cohesive element size in particular is critical (Foulk 2010;Seagraves and Radovitzky 2010). Too large elements cannot accurately resolve the fracture process, whereas smaller elements increase computational effort. However, choosing an appropriate element size requires knowing the maximum size of the FPZ for a given scenario, the FDPZ. More precisely, the FDPZ is the zone when the half crack opening displacement at the traction-free crack tip reaches the critical separation of the cohesive model. Currently, the FDPZ size has to be computed and is rarely known. An exception is the work by Ha et al. (2015) who estimated its size for bending-and compact tension tests, for which however no analytical solutions exist.
Furthermore, a cohesive law is required. It describes the traction-separation relationship within the FPZ and is considered a material property, which needs to be identified with either lab-or field-measurements. For such measurements, an initial estimate about the size of the FDPZ is needed for sufficient sensor deployment around the fracture process zone. Figure 1 shows an example of a large-scale sea ice fracture experiment to measure the fracture properties, and the related displacement sensors arrangement ahead of the physical crack tip (see the physical crack size B in Fig. 1). This is also required for measurements of cohesive fracture properties for other quasi-brittle materials, e.g. concrete, rock, or ceramics (Bažant and Planas 1998). Moreover, because most of the test samples are of lim- A third and related issue is the scaling of small-scale fracture test results to full-scale while respecting the size effect, i.e. the dependence of the nominal strength of geometrically similar structures on size. For an indepth discussion of the size effect and different underlying theories see e.g. (Bažant 2005). Testing of full-scale geometries such as large ice floes or concrete structures is often not possible due to the required effort. Most tests are done with much smaller specimens. Hence computations, e.g. for the design of large structures, must rely on extrapolation of test results, see e.g. the discussion for concrete testing by Bažant (2002). The CEM should be able to capture the size effect and is a potential tool for this extrapolation (Lu et al. 2015b;Bažant and Yavari 2005;Elices et al. 2002). This was for instance confirmed by Morel and Dourado (2011) and Ha et al. (2018), though for beam tests and a limited range of specimen sizes.
Hence the first objective is creating and validating a CEM-based numerical model for finite plates. Analytical solutions based on Linear Elastic Fracture Mechanics (LEFM) are used to validate the force required to achieve crack propagation in a finite plate (Lu et al. 2015b). More importantly, both force and the FDPZ size for large plates are compared to the work by Wang et al. (2019b), who studied a cohesive edge crack in a semi-infinite plate. The authors derived the length of the FDPZ (R) and its evolution with the physical crack size (B) under a pair of concentrated forces (P) at the crack mouth. Their solutions can be seen as an analytical benchmark to the application of any Cohesive Zone Method based numerical method, e.g. (Lu et al. 2012;Turon et al. 2007;Park et al. 2012;de Borst 2003;Elices et al. 2002;Falk et al. 2001;Remmers et al. 2003;Kuutti et al. 2013;Paulino et al. 2008;Molinari et al. 2007;Unger et al. 2007;Zi and Belytschko 2003;Moës and Belytschko 2002).
After the validation, the second objective is to extend the results of Wang et al. (2019b) to that of a finite square plate as well as identify the unknown FDPZ size and its dependence on plate size and crack tip position. The third aim is to investigate the ability of the CEM model to capture the size effect with respect to the material strength and plate size.

Cohesive edge crack in a square plate
Consider an Edge-Cracked Square Plate (ECSP) (see Fig. 2) with finite width W = L, length L and thickness t. The cohesive law results in Eq. (1), which theoretically eliminates the stress singularity at the cohesive crack tip A.
Where K is the stress intensity factor and coh the normal tensile traction in the cohesive zone. H r (A, X ) is a weight function for the ECSP of related width to length ratios W/L and given by Dempsey and Mu (2014). Detailed information on the weight function is given in the Appendix.
Following the notation of Wang et al. (2019b), two cohesive stress profiles coh (A, B, X ) are considered, the rectangular softening in Eq. (2) and linear softening in Eq. (3).
0 is the local tensile strength, U is half of the crack opening displacement (COD) and U cr half of the critical separation at which the traction in the cohesive zone becomes zero, B is the physical crack length. When the cohesive size R = A − B is much smaller than the cracked body, i.e. R L, small scale yielding can be assumed and the cohesive crack model yields the same results (in terms of peak external force (X )| Peak and crack profiles U (A, B, X )) as those predicted by Linear Elastic Fracture Mechanics (LEFM). However, when the size L becomes comparable to the cohesive size R, the cohesive crack model deviates from the LEFM prediction.
For an edge cracked semi-infinite plate under a concentrated force P at the crack mouth, i.e. (X ) = P/tδ(X ) = P * δ(X ), Wang et al. (2019b) defined that the Fully Developed Process Zone (FDPZ) takes shape when U (A, B, B) = U cr ; and solved the FDPZ's evaluation with the physical crack length B.
Here, we extend the work to studying the finite ECSP by numerical simulations under the same P * δ(X ) at the crack mouth (Fig. 2). For generality, a characteristic length λ = π E U cr /(2 0 ) is introduced to normalise all the spatial terms, where E = E for plane stress and E = E 1−ν 2 for plane strain, ν is Poisson's ratio. The plate sizes are multiples of λ L , i.e. L = nλ L , where λ L = G c E / 2 0 and G c the energy release rate 1 . The definition of λ L is more general such that its numerical value does not change for different shapes of tractionseparation laws (TSL) for a given material (albeit with the same G c , E and 0 ). In this regard it differs from the previous λ = π E U cr /(2 0 ), which depends on TSL shape. Nonetheless they have the same meaning: they are both characteristic lengths depending solely on material properties. In this paper, for simplicity, we use λ L to normalise the size of the ECSP but express the FDPZ size in terms of λ; the results are interchangeable. See also corresponding section in the Appendix.
When the plate size L/λ L gets very large, the solution of the ECSP will converge to that predicted by Wang et al. (2019b). This is the benchmark of the current numerical simulation before new results concerning finite size ECSP are presented.

Cohesive element method
The cohesive zone model (or the fictitious crack model) was initially proposed by Hillerborg et al. (1976), who examined available fracture mechanics theories (the stress intensity factor, the energy balance approach, and the Dugdale and Barenblatt approaches) to describe crack initiation and propagation by means of the FEM. It models accumulated, localized damage as effective behaviour in the Fracture Process Zone (FPZ) ahead of the crack tip. In the process zone, also called cohesive zone, fracture is represented as a gradual process of separation of two virtual surfaces, resisted by tractions between these surfaces. It is a rather simple phenomenological model to generalize the Fracture Process Zone behaviour with different types of material bonding (Bažant 2008). It assumes a traction-separation relationship (cohesive law) to describe the line-like (2D) or surface-like (3D) cohesive zone (see Fig. 3). The benefit of this method is the elimination of the stress singularity at the crack tip (Eq. (1)) which makes its integration with existing numerical methods (e.g., the Finite Element Method, FEM) rather straightforward.
The cohesive zone theory, i.e. Eq. (1), can be integrated into the Finite Element Method through the so-called Cohesive Elements (CE). The constitutive behaviour of the CEs is described by the tractionseparation law (TSL), i.e., coh (x) versus U (x) in the FPZ (i.e. α − β). Various shapes of TSL exist for different purposes, e.g. ductile or brittle fracture. Eqs. (5) and (6) are two common traction and separation laws. In the numerical model these are slightly modified and used as intrinsic TSL, see Fig. 3. Intrinsic indicates a finite initial stiffness or slope of the TSL K coh . This is required to maintain compatibility and momentum transfer across elements because cohesive elements exist in the simulation model prior to fracture. By contrast, an extrinsic approach assumes an initially rigid response. This would require a dynamical insertion of cohesive elements during computation (Seagraves and Radovitzky 2010).
Here, to directly utilise the CEM framework in existing software and with little deviation from the theoretical requirement, we adopt the intrinsic CEMs and set the initial stiffness K coh to be as large as numerically admissible. Prior to reaching the peak traction 0 , the cohesive elements have a reversible elastic response.
For linear softening, passing the peak traction 0 initiates a damage process which decreases the element's stiffness, see left-hand side of Fig. 3. This model is based on (Dávila and Camanho 2001). For rectangular softening, the traction plateaus at the peak traction, see right-hand side of Fig. 3. This model is based on the work by Tvergaard and Hutchinson (1992) and its implementation as described in (Sandia National  ). The original model is a trapezoid with finite slopes on both sides. Here, there is no righthand finite slope which is a better approximation of the TSL from the analytical solution in (Wang et al. 2019b). The conversion from the variables of the original Tvergaard-Hutchinson model to the material properties used here is given in the Appendix.
After reaching the critical separation 2U cr the element is deleted. The area under the curve is the energy release rate G c . Cohesive models are usually defined in three directions, here only mode I is activated, i.e. separation normal to cohesive surfaces.
The size of the FDPZ is usually not known beforehand and needs to be solved in accordance to the equilibrium stated in Eq. (1). Practically, we insert cohesive elements along the centre line of the ECSP (see Fig. 5) and obtain its size by identifying the spatial locations of α and β at the critical moment of 2U (α, β, β) = 2U cr , i.e. when the leading element at the physical crack tip (at β) reaches the predefined critical separation 2U cr .  Table 1 Geometries of the simulated ECSP models. λ L = 1.111m, values are approximate

Numerical model
To achieve the two objectives of this paper, a square plate with different sizes is simulated. The cases are listed in Table 1. The numerical model maps the cohesive edge crack in a square plate (ECSP) with finiteand cohesive elements. The software LS-Dyna with an explicit time integration was used because future applications of similar models seek to simulate dynamic scenarios, e.g. a collision between an ice floe and a ship. Nevertheless, implicit solving should also be possible.
The plate consists of Belytschko-Lin-Tsay shell elements with unit thickness, one through-thickness integration point, and a Poisson's ratio of zero. The to-be-calculated FDPZ size R is scaled with the characteristic length λ L ; and the numerical accuracy of the calculated FDPZ is largely dependent on the cohesive element size. The cohesive element size should be small to capture the stress gradients along the crack, see Fig. 6. Through trial and error, a cohesive element size of smaller or equal to ≈ 0.007λ L was chosen, which balances accuracy and computational effort. This also indicates the spatial resolution of the results. The cohesive elements are inserted along the prescribed crack path and connect the edges of the shell elements. See (Livermore Software Technology Corporation 2019) for element formulations.
In the following, the length over which cohesive elements are inserted is termed D, see Fig. 5. A small initial crack of a length equal to 6-11 cohesive ele- ments is prescribed in the model to apply the crack mouth opening displacement (CMOD). The displacement is applied to one element less than initial crack length, i.e. over the length of 5-10 cohesive elements, see Table 1. This mimics the loading condition of (X ) = P/t δ(X ) = P * δ(X ) at the crack mouth. The splitting process is displacement controlled for numerical stability. The displacement is applied with an almost zero initial slope to avoid potential dynamic behaviour, see Fig. 4. Global mass weighted nodal damping is applied for the same reason but turned out to not influence the results. See (Livermore Software Technology Corporation 2019) for a description of the damping mechanism. The force required to achieve this displacement is the splitting load. One node in the middle of the far end of the plate is fixed.
An exemplary mesh is shown in Fig. 5. The first rows of elements around the crack path are rectangular shells. All other elements are triangular shells which allows for a rapid increase of element size away from the crack and keeps the total number of elements within a reasonable range. For plates of sizes L ≤ 30λ L cohesive elements are inserted throughout the whole plate, i.e. D = L. For larger plates D is capped at 30λ L . This is because the number of cohesive elements drives the total number of elements due to their small size and the small size of surrounding elements.
A simple elastic material model is used for the shell elements. Its parameters are E = 10 10 Pa, ν = 0, and a density of 900 kg · m −3 . User-defined material models are used for the cohesive elements based on the above-mentioned TSLs, see (Dávila and Camanho 2001;Tvergaard and Hutchinson 1992;Sandia National Laboratory 2003). The fracture energy G c and the initial stiffness K coh is the same for both TSLs. K coh is several magnitudes larger than the Young's modulus of the bulk material. Both user-defined material models were validated against material models from the LS-Dyna library, see Appendix. Values for parameters and properties are given in Table 2.

Simulation results processing
We capture the potential dynamics of crack propagation through explicit time integration. When a cohesive element reaches the peak traction or the maximum separation (points α and β in Fig. 6) an output with the time, position and its element id is generated. For every point in time where an element failed, the closest point in time where any other element reached peak stress is identified with a bisection algorithm. The X -coordinates of the element that reached peak traction and the element that reached maximum separation are used to compute α and β, respectively. Then the cohesive zone size is 2 R/λ = α − β. All results are normalised with λ or λ L , where λ is not the same for the different TSL, but λ L is. Non-normalised parameters are given in the Appendix Table. 2.
This yields our results of interest, i.e. R and P * versus the varying physical crack length B. Fig. 7 illustrates the cohesive stress profile coh 's distribution with varying β values (i.e. a running crack). The stress progresses as expected. Once the crack tip is close to D, the lack of additionally inserted cohesive elements along the "future" crack path distorts the stress distribution (β ≥ 8).

Validation with LEFM
First, the simulated splitting load is compared to an analytical solution for finite plates. When the size L/λ L of an ECSP gets large, the CEM-based simulation should yield the same results as those by LEFM. Based on LEFM theory, which basically means coh = 0 in Eq.
(4), this leads to the formulation for the normalised peak force in Eq. (7).
As an example, the normalised weight function h r (α, 0) for an edge cracked semi-infinite plate is given by Wang and Dempsey (2011) and repeated in Eq. (8).
In Fig. 8, a normalised splitting force P/t K √ L is plotted over a normalised crack tip position. Analytical solutions are given for a finite plate and it is natural to choose the physical crack tip position as the crack length. In contrast to LEFM, the fracture process in the CZM is not concentrated in one point. Theoretically, as the size of the plate gets large, the CZM model and LEFM model should yield the same results. The larger the FPZ size compared to the plate size, the more the solutions are expected to diverge from the analytical solution.
Simulated and analytical solutions are given in Fig. 8. The analytical solutions for the semi-infinite and finite plate diverge after B/L ∼ = 0.1. The larger the plate, the better the match between LEFM and CZM solutions. However, if the tail from the stress distribution (right hand side in Fig. 6) approaches D, the length along with cohesive elements are inserted, the results  Fig. 7. Hence for plates larger than 30λ L solutions are not independent of D anymore and therefore not shown in Fig. 8. Results for L ≥ 30λ L , and with D = L, are only used for computing the cohesive zone size until it reaches a plateau (Figs. 11 and 17). The plateau is reached before the stress distribution reaches D. On the lower end, the solution for the 1λ L plate does not match the LEFM solution. Further, close to B/L = 1 at the right-hand side of the plot, we can also see a small distortion due to the stress distribution reaching D.

Results
After the validation of the numerical model, the major results are presented in this section.

Size effects simulated by the CEM
According to the cohesive zone theory, there is a nonlinear size effect regarding the strength of the cracked body (Elices et al. 2002). The comparison made between the CEM and LEFM in the preceding section only makes sense when the sample size is large. For an ECSP of arbitrary size it is convenient to compare the CEM-based results with an available cohesive zone theory based analytical solution concerning the crack initiation force P/t K √ L. Fig. 9 Splitting force over plate length for a fixed ratio of B/L = 0.043 for linear softening. CZM curves are both normalised with λ L for comparability. Analytical solutions are from Eqs. (7) and (9) In cohesive zone theory, for the linear softening case (i.e., Eq. (6)), the solution to the peak splitting force can be reduced to an eigen-value problem in Eq. (9), where U(α, x, 0) is the normalised Green's function defined in Eq. (10). Equation (9) is derived in the Appendix.
The material strength is defined as P/t L in accordance with (Bažant 2005). The solutions of Eqs. (9) and (10) are compared with the CEM-based simulation results in Fig. 9 illustrating the size effect, i.e. the material strength's variation with the size of geometrically similar structures (i.e., L/λ L ). The results presented in Fig. 9 are only for the case with a physical crack length of B = 0.043L. However, the same trend can be found for other physical crack lengths as well, e.g. see reference (Lu et al. 2015b) for the case with B = 0.3L (see also Fig. 17 in that reference).
Generally, when the plate size L → ∞, the structural strength P/t L is scaled with √ L according to the LEFM scaling, thus the 2 : 1 slope for all curves at the lower-right end of Fig. 9. In this case, with B = 0.043L, we see the convergence from CZM to LEFM starts from around L = 20λ. For the case with , the convergence appears to take place from L = 12λ according to (Lu et al. 2015b).
When the size gets smaller, the predictions by CZMbased methods start to deviate from the LEFM-based results (black dashed line). In Fig. 9, three CZMbased methods are presented. These are (1) the CEMsimulated CZM with linear softening (blue solid line with triangular symbols), (2) the CEM-simulated CZM with rectangular softening (green dotted line with rectangular symbols), and (3) the CZM-based analytical solution with linear softening (dashed red line). All in all, the CEM-simulated linear softening result and its corresponding analytical solution coincide, signifying the correctness of the CEM-based implementation. Additionally, to the very left, we see that the rectangular softening based results predict a stronger material strength over those based on linear softening. This is an expected outcome as the rectangular softening characterises a stronger material compared to its linear softening counterpart.
For smaller ECSP sizes L ≤ 0.3λ L the behaviour is different. The external force already decreases when the leading cohesive element reaches critical separation, see Fig. 10. With increasing CMOD, the FPZ size increases until it spans the entire crack line (i.e. the D region in Fig. 5). The external force at the CM also increases, but at its maximum the leading cohesive element has not failed yet. Instead, when the critical separation is reached, the external force is already in a downward slope. From the point where maximum external force is reached (between 2 and 3 s), the separation accelerates, a "catastrophical failure". Achieving a stable crack growth in the simulation is challenging when the fractured sample is too small. Hence in Fig. 9 the maximum for P/Lt is taken instead of the value at B/L = 0.043 for the three leftmost points, i.e. the three smallest plates.

Cohesive zone size
Plots of the cohesive zone size versus the physical crack tip position are given in Figs. 11 and 17. The dashed black line is the analytical solution for semi-infinite plates, as derived in (Wang et al. 2019b). The horizontal dashed black line is the limiting case of a fully developed process zone as A → ∞. Other coloured curves are numerical model results. For plates with L ≥ 30λ L , the results are cut off when B is close to D for reasons given above (Sect. 3.2).
In the case of linear softening (Fig. 11) the curves for the larger plates with L ≥ 80λ L appear to approach the limiting case from the analytical solution for the range of B/λ simulated here. For the smaller plates, the cohesive zone size increases, reaches a maximum, and then decreases as the physical crack tip approaches the end of the plate. The curves for plates L ≥ 30λ L follow the analytical solution. The smaller plates deviate earlier from the analytical solution due to their size. The results for rectangular softening are very similar (see Appendix, Fig. 17). Maximum cohesive zone sizes (the plateau values, R max ) are given in Fig. 12 for both traction-separation relationships. The larger the plate, the closer the maximum cohesive size to the analytical limit case for semiinfinite plates. The R max values are larger for rectangular softening. This is in line with theory, as they represent an upper boundary compared to the linear softening case. The values for the smallest plates should be treated with caution as these model sizes are very small in the context of the applied methodology.
Two opposite cases, indicated in Fig. 12a, b, are schematically shown in Fig. 13. Figure 13a illustrates R max for a smaller plate, where the FDPZ almost covers the whole plate, whereas Fig. 13b shows what happens if R max approaches its limit case for large plates.

Splitting force
A different plot of the splitting load over crack tip position is given in Figs. 14 and 18 for linear and rectangular softening. These plots contain the same information as in the LEFM comparison plot, Fig. 8, but with a different normalisation, in line with the analytical solution. The dashed black line is the analytical solution for an edge crack semi-infinite plate (Wang et al. 2019b). The forces for the larger plates follow the analytical solution until ≈ 3 − 4β, whereas the forces for the smaller plates reach their maxima quickly and slowly decrease after. The forces for the small plates with L ≤ 1λ L don't plateau and decrease immediately. Overall, the forces are as expected for finite plates. The behaviour is similar for linear and rectangular softening.

Cohesive zone size for ice
Coming back to the motivation of this work, the first major contribution is the maximum cohesive zone size presented in Fig. 12. The results should be valid for all types of quasi-brittle materials, e.g., rock, cement concrete, ceramics and ice. Nonetheless, owing to the authors' background, this section uses ice as a prac-tical example of applying Fig. 12, for which material properties are given in Table 3. For lab ice, for linear softening, and using Eq. (38), we have λ = 0.0565m and λ L = 0.036m. From the plot in Fig. 12 we take point (L/λ L = 5.0; R max /λ = 0.358), so L = 5.0 · 0.036 = 0.18m. Then, R max = 0.358 · 0.0565 = 0.02m, or about 11% of the specimen length L. The same calculation for sea ice with λ = 0.377m and λ L = 0.24m gives R max = 0.135m.
The same applies to rectangular softening. Take the same point as before (L/λ L = 5.0; R max /λ = 0.4105). Then, L = 0.18m. Now, due to the different λ we have R max = 0.4105λ = 0.4105 · 0.0283 = 0.0116m, or about 6% of the specimen length L.
As λ is material-specific, estimating R max like this can help to design experiments similar to the model setup (Fig. 2). Exemplary tables for R max for laboratory and sea ice are given in the Appendix.

Cohesive element size
In any CEM-based simulation, the cohesive element size must be smaller than the FPZ size (Falk et al. 2001;Turon et al. 2007). How much smaller is subject of ongoing research and was not the focus of this work, but some remarks can be given.
As the cohesive zone size varies, any element size criteria based on its maximum (see exemplary calculations in above section) is an upper bound. That is, it is likely necessary to use smaller elements to capture the cases when B → 0 and B → L and the cohesive zone size R is small. Plots of the cohesive zone size, e.g. Fig. 11, as well as cohesive stress profiles along the crack such as Fig. 6 can help choosing the element size for simulating cracks in plates similar to our model (Fig. 2).

Splitting force and size effect
The splitting force versus the physical crack length is compared between the analytical solution and CEMbased numerical simulations. Two sets of analytical solutions are adopted: (1) the splitting force of a finite ECSP in Fig. 8; (2) the splitting force for a semi-infinite plate in Figs. 14 and 18.
In both scenarios, the analytical solutions are a manifestation of LEFM solutions. This is because firstly, the analytical solution in Fig. 8 is based on LEFM theory (see Eq. (7)); and secondly, in Figs. 14 and 18, the corresponding cohesive zone sizes are rather small in comparison to the plate size. Theoretically, we expect in the CEM-based simulations, as the size of the cracked plate increases, that the normalised splitting force (either P/t K √ L or P/tλ 0 ) converges to those analytical solutions. This is confirmed in the comparisons, showing our CEM-based numerical model works as intended.
Moreover, the CEM-based numerical model captures the size effect, at least in this setup, see Sect. 5.1 (Fig. 9 is strictly valid only for B/L = 0.043 and the given TSLs). It is another potential tool to extrapolate lab-scale measurements to field scale in addition to the "CZM + weight function" approach developed in (Lu et al. 2015b). This includes the parameters of the traction separation law, which are the same for all model sizes used here, similar to results by Ha et al. (2015). The discretization around the crack line (D in Fig. 5) is the same for sizes ≥ 1λ L . Below 1λ L element length must change to accustom the small geometries. The element sizes used here are too small for simulations of full-scale scenarios, also considering that only a single crack was simulated. Nevertheless it should be possible to also capture the size effect with a coarser mesh, especially if the aim is not to calculate accurate cohesive zone sizes, see also (Turon et al. 2007). Overall, it seems that the size effect is captured as long as the FPZ is sufficiently resolved by the element size. These results are encouraging, despite known problems of numerical CEM-based models (Seagraves and Radovitzky 2010; Rimoli and Rojas 2015;Lu et al. 2014).

Fully developed cohesive zone
The fully developed cohesive zone is defined as the size of β−α as in Fig. 6 when the leading cohesive element's separation reaches the critical separation 2U = 2U cr . From the size effect validation case in Fig. 9, for a large ECSP, we can extract the maximum splitting force reached so far before the leading cohesive element's failure right when its separation 2U = 2U cr (and the traction-free crack opening is at B/L = 0.043). This is before the peak splitting force is reached, since splitting force is increasing monotonically, even after the failure of the first cohesive element.
However, for the small cracked plate, e.g. L = 0.1λ L , the maximum splitting force (during the lifetime of the leading cohesive element) does not occur when the cohesive zone is fully developed (i.e., see Fig. 10, the peak force is reached before the failure of the leading cohesive element). Instead, when the cohesive zone is fully developed, the splitting force is already decreasing. This means that a crack would propagate in a small ECSP even with a decreasing splitting force, leading to an unstable crack propagation scenario. This indicates that it can be rather challenging to obtain a stable crack growth in an experiment when the fractured sample is too small (e.g., lab-scale experiments).

Static and running cracks
The cohesive zone size is extracted from a running physical crack (B) shown in Fig. 7 but it is compared to analytical solutions that are based on static analyses.
In the simulation, a slow application of prescribed displacement was applied to reduce dynamic effects. Yet this approach is limited since a slower displacement does not reduce the minimum time step but increases the number of time steps needed to achieve the same displacement, which in turn increases computation time. So although a rather slow displacement controlled loading scenario is simulated, a running crack will potentially bring in inertia effects (Seagraves and Radovitzky 2010).
However, given the satisfactory agreement of the cohesive zone size in the limiting scenario (i.e., favourable agreement in Figs. 11 and 17), we expect the presented results for other plate sizes are not far from their corresponding static solutions. Moreover, the kinetic energy was close to zero in all simulations.
So, for the time being, it is unknown what the exact difference between a static and a running crack's cohesive zone size is. Nevertheless, we expect the numerical results presented in Figs. 11 and 12 to be of practical interest for most lab experiments.

Numerical errors
The mapping from the analytical to the numerical model introduces sources of error. Despite the agreement between the numerical results and the analytical solutions in limiting scenarios, potential error sources are emphasized.
Generally, the analytical model is essentially a static 2D model, whereas the numerical model is solved as a dynamic quasi-2D (see Sect. above and 3.1) model with explicit time integration. The numerical results are also affected by discretization and element formulations. This was checked with a mesh convergence study and different element configurations, which did not change the outcome. If R max is compared to the cohesive element size, there are mostly 40 to 50 cohesive elements or more in the process zone, but not less than 26. This should be sufficient for accurately resolving the cohesive zone (see also Sect. 7.1), but more research is needed for general element size recommendations.
Also, the analytical solutions are based on a concentrated force at the crack mouth. This is not the case in the numerical model as a point force often leads to numerical errors, e.g. a distorted stress field near the corresponding node. Nonetheless, the area over which the force was applied is magnitudes smaller than the plate size, see Table 1.
Lastly, the numerical material models deviate slightly from the analytical behaviour due to the issue of initial stiffness, as described in Sect. 3. The plateau in the traction-separation law for the rectangular softening material model, see Fig. 3, leads to some oscillations in force and cohesive zone size R in the numerical results. Also, R sizes can be slightly bigger than for the semiinfinite analytical solution, which should not happen for the finite plate numerical model.

Conclusions
In a nutshell, the numerical CEM model works as expected. Several points can be derived from the results: -For this setup, the CEM-based numerical model captures the size effect. Like the CZM + weight function method (Lu et al. 2015b), it is a potential tool to extrapolate laboratory scale measurements to field scale. -For finite plate sizes of practical interest, Fig. 12 presents the maximum cohesive zone sizes R max in dependence of plate size. This is the major contribution of this paper and the results presented have an accuracy of about 0.0072λ L . λ L is a characteristic length and depends only on material properties. -In general, the fully developed process zone (FDPZ) size for finite plates is influenced by the cracked plate size and physical crack length. -With increasing physical crack length within a finite edge cracked plate, the FDPZ first increases, then plateaus to a maximum value R max , before it decreases as the physical crack reaches the far end boundary of the finite plate. This is illustrated in Figs. 11 and 17. -The plateaued FDPZ size R max increases with the plate size and with L → ∞ converges to 0.5λ and 0.465λ for the rectangular and linear softening traction and separation law, respectively. λ is a characteristic length and a material property, similar to λ L . -The presented CEM-based numerical model can be implemented to evaluate the cohesive zone size evolution for other cracked geometries.
Acknowledgements Open Access funding enabled and organized by Projekt DEAL. The authors would like to thank for the financial support by the following institutions: Leon Kellner was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-441262697. Wenjun Lu was funded by VISTA-a basic research programme in collaboration between The Norwegian Academy of Science and Letters, and Equinor.
Author contributions LK and WL conceived the study. LK created the numerical models, programmed the material model subroutines, designed and programmed the postprocessing and visualisation algorithms, ran the simulations, did the postprocessing, and contributed to the writing. WL provided the majority of the theoretical framework, provided the analytical solutions, and contributed to the writing. SE and KVH supervised research and contributed to the writing.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The data used to create the plots, i.e. the data resulting from the analyses, is available, see (Kellner and Lu 2020). Alternatively, please contact Leon Kellner, leon.kellner@tuhh.de, or Sören Ehlers, ehlers@tuhh.de.

Code availability
The postprocessing and visualisation codes are licensed under the GNU General Public License (gnu.org/ licenses/gpl-3.0.html), if you would like to obtain them, feel free to contact Leon Kellner, leon.kellner@tuhh.de.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Weight function method
In this paper, the weight function method, which is independent of external loading profiles, is adopted to calculate the Stress Intensity Factor (SIF). The weight function method developed by Bueckner (1970) and Rice (1972) characterizes the SIF for a cracked symmetric body with a crack length of A under a unit loading at a location X . For an arbitrary loading profile, its corresponding SIF K (A) can be calculated with Different cracked geometries have different weight functions H r (A, X ). Here, we adopt the weight function developed by Dempsey and Mu (2014) for the considered centrally cracked square plate. The weight function H r (A, X ) has the general form in the following equation in which the function needs to be separately established for specific cracked geometries, see e.g. (Wu and Carlsson 1991).
Note with A/L and l/λ, also A/L = α/l. For this paper, the explicit form of G i (expressed as a collection of other functions involving the characterization of crack surface displacement and crack surface area) can be found in the work by Dempsey and Mu (2014) and is not repeated.
One of the important properties of the weight function is that, if we introduce a scaling α = A/λ, β = B/λ and x = X/λ, we have h r (α, x) = √ λH r (A, X ). When we have a concentrated force acting at the crack mouth x = 0, the weight function reduces to h r (α, 0). For example, one simple form of the weight function for the cracked semi-infinite plate is expressed in Eq. (8). When it is under a pair of splitting forces acting at the crack mouth, we simply need to introduce x = 0 to obtain its corresponding weight function.
Derivation of Eq. (9) The derivation of Eq. (9) follows the procedures in the original works of Li and Bažant (1994) and Li and Liang (1993). Similar derivations can also be found in (Bažant and Planas 1998) and (Wang and Dempsey 2011). A detailed derivation of Eq. (9) is presented in Appendix A in (Lu et al. 2015b). For completeness, these derivations are repeated below.
According to the concept of cohesive zone theory, we have The general expression for the half COD can be written as Rearranging Eq. (15) leads to The Green function U(A, X, S) is defined as For a linear softening cohesive crack, we have Inserting the right-hand side of Eq. (18) into the lefthand side of Eq. (14), we have The two terms on the right-hand side of Eq. (19) (20).
According to the Leibniz integral rule, taking variations over A upon Eq. (22) leads to By definition, at peak value, δ P = 0 (Bažant and Planas 1998). Therefore, the first term on the right hand side of Eq. (23) vanishes. In addition, due to symmetry, U(A, X, A) = U(A, A, X ) = 0 in the third term. The fourth term also vanishes as δ B/δ A = 0. Because of the relationship in Eq. (24), which can be introduced to the second and sixth term in Eq. (23), these two terms cancel each other by virtue of Eq. (13).
Eventually, Eq. (23) is simplified into Eq. (25). (25) Replacing δ coh (X ) with a proportional cohesive stress profile (X ), we can obtain Eq. (26), which is a typical eigenvalue problem to solve. After identifying the eigenvector (X ), we can begin to calculate the peak splitting force P/tδ(X ).
Recalling that for linear softening the cohesive stress could be written as in the left-hand side version of Eq. (18). Multiplying both sides of this formula with the obtained eigenvector (X ) and integrating within the cohesive zone (i.e., B → A) lead to Eq. (27).
Eq. (29) is simplified into Eq. (31), which finally leads to Validation of user material models User material models must be used to obtain element information at the points of peak stress and failure (α and β). They were programmed as LS-Dyna subroutines (Erhart 2011). The user material models were validated against standard material models which are available from the LS-Dyna library. The linear softening model was validated against the cohesive mixed The Tvergaard-Hutchinson model characterises separation with a dimensionless gap vector . The separation for the initial and final peak tractions are 1 and 2 , the critical separation is fail . These are related to the material properties based on the area of a trapezoid without right-hand side slope, i.e. 2 = fail (Fig. 3). Both separations are normalised with u fail such that / fail = 1 for the critical separation.
The forces were compared for the same model with user defined and reference material model. Models with plate sizes of 5λ and 1λ were used. The comparison for the linear softening model is shown in Fig. 15. The curves match. The comparison of the splitting forces for rectangular softening and the Tvergaard-Hutchinson model is shown in Fig. 16. Again, due to the zigzag behaviour the curves were smoothed. The curves for the rectangular user-and reference material model also match. Relation between λ and λ L Two definitions are used for the characteristic length, one for normalising results λ, and a more generally defined one for describing lengths λ L . They are related as follows: For the linear case, since the fracture energy is G c = 0.5 0 2U cr = 0 U cr , and for the rectangular case, because the fracture energy is G c = 0 2U cr = 2 0 U cr , Additional plots The plots in Figs. 17 and 18 are analogous to Figs. 11 and 14 but for rectangular softening.

Fig. 17
Cohesive zone size for different plate sizes and rectangular softening, analytical solution including limiting case from Wang et al. (2019b). The cohesive size of the large plates (100λ L and 80λ L ) plateaus towards the limiting case of the semi-infinite plate's scenario. The range of β values is different compared to Fig. 11 due to λ being different for linear and rectangular softening. The plateau in the traction separation law for the rectangular softening material model (see Fig. 3) leads to some oscillations in force and cohesive zone size R in the numerical results. Hence curves are smoothed. Raw data is given in lighter colours  Table 2 gives non-normalised parameter values used for the elastic material model surrounding the crack as well as the cohesive models. K Ic is used for the 'LEFM + weight function' curve in Fig. 9 and calculated with √ G c E . Table 3 gives ice material properties  (2009) and for sea ice from Dempsey et al. (1999a). Table 4 gives exemplary cohesive zone sizes, see Sect. 6.