Unraveling the myth of closure corrections Sharpening the definition of opening and closure stresses with an energy approach

The substantiation of fatigue crack closure corrections is disputed, based on the closure stress definition. The Δ K eff equation lacks a physical explanation. An inconsistency is observed between the opening stress S op phen as used by this equation and the physical opening stress S op phys . This S op phys is related to S op phen through an energy equivalent area approach. Furthermore, an elastic spring model is used as a physical approach to crack closure effects. An FEA approach generates S op phys values, which are reworked into S op phen . This physical model agrees well with existing closure corrections, and is able to provide a physical explanation for their necessity.


Introduction
Fatigue crack growth in metals can be described using linear elastic fracture mechanics (LEFM). One of the key concepts of LEFM is the range of stress intensity factors ΔK indicating the severity of the stress distribution around the crack tip, as function of the applied far field stress range ΔS and half crack length a.
In 1961 Paris et al. [1] observed that the crack growth rate da/dN is a function of both ΔK and the stress ratio R. The R dependency can be accounted for by replacing ΔS with an effective stress range ΔS eff = S max -S op , which results in the effective stress intensity factor: The stress level S op is considered the stress level corresponding to the first moment in the loading cycle where the crack tip is fully opened. The existence of such a stress value is widely reported, among others in Refs. [3][4][5][6][7][8][9].
The use of the effective stress range results in crack growth curves collapsing onto each other. Elber [10] was one of the first to relate this reduced stress range to the observed phenomenon of crack closure. The current state of the art models [11][12][13] use the stress intensity factor ΔK eff as similitude parameter, because different fatigue cases can then be compared using the unique relationship between ΔK eff and da/dN, independent of R. The similitude parameter ΔK eff is often linked to the crack growth rate using a power law, such as the Paris equation, Ref. [14], which is a purely empirical correlation, for which the physical explanation is unknown, Ref. [2].
In Ref. [10] it is reported that during the loading phase of a fatigue cycle an initial nonlinear relationship is observed between the crack opening displacement and the applied stress. This initial nonlinearity is ascribed to crack closure effects, imposed by plasticity. During unloading the plastically deformed area around the crack tip can close before S min is reached, and similarly during loading the crack tip starts to open at a stress level above S min . This crack tip plasticity plays an important role in crack closure, yet it does not correspond well with the LEFM theory on which the ΔK eff parameter is based. A main LEFM assumption is that plasticity is concentrated in an infinitesimally small area at the crack tip. The plasticity related contradiction between LEFM and the ΔK eff parameter therefore suggests that ΔK eff and related closure corrections are not complete. Moreover, ΔK eff is often assumed to be the driving force for crack growth, while this statement has little physical basis: it is not a force, but it rather is a representation of stress effects at the crack tip.
There is no reason why the ΔK approach using ΔS would give a correct result, as the method is not physically correct. The improved approach of ΔK eff using ΔS eff with S op still suffers from the same oversight, it tries to correct a method that was not physically correct to begin with. In this paper a discrepancy in opening stress values and closure corrections is explained, and a physical method is presented to transform true opening stresses to S op values used in the ΔK eff closure correction approach. As a first step, the incompleteness of closure corrections based on this ΔK eff approach are further elaborated on below.
Eq. (1) is a function of S op to account for closure effects. Closure corrections such as Elber [10], Schijve [15], De Koning [16], and Newman [17], describe closure effects using the nondimensional parameter S op /S max . The latter three predict crack closure for all stress ratios R, even though there is a range of R values where the crack tip never closes during the full load cycle, as suggested by Refs. [18][19][20][21]. Three other closure corrections, Refs. [22][23][24] report similar closure behavior but for steel instead of aluminum, suggesting that the existence of non-closure near R = 1 is not material dependent. Furthermore, all seven mentioned closure corrections start to differ significantly for negative R values, showing that there is no true consensus on the crack closure behavior at compressive S min .
Environmental effects on crack surface roughness are not a significant cause of the spread in closure corrections, as shown by Refs. [26][27][28]8] for different metals and different environments.
A large literature base, Refs. [29][30][31][32][33][34][35][36][37], suggests that known closure corrections are not sufficient to collapse da/dN versus ΔK curves for various R onto a single da/dN versus ΔK eff curve. Refs. [30,38] discuss a theoretical model to account for asperity effects behind the crack tip, showing that closure is only significant when the crack is fully closed, and that the measured asperity effects are small. Vasudevan et al. [30] note that despite several decades of literature, there seems to be no accurate method of observing crack closure.
According to Refs. [39,17] S op /S max varies with ΔK max /K o , and Ref. [40] observed the same with a strip-yield model. This suggests S op / S max = f(R,S max ,S yield ), implying that, according to Refs. [10,[15][16][17], all crack closure corrections as f(R, S max ) cannot accurately or uniquely describe closure for every R and metal (or isotropic material), and that S yield or σ 0 needs to be taken into account.
If S op /S max = f(R, S max , S yield ), it appears logical that there is also a finite width effect. The finite width raises the net-section stress in the crack plane, increasing the stress around the crack tip even more as the crack grows. Therefore it could be that S op /S max = f(R, S max , S yield , a/W).
The incompleteness of common closure corrections to describe closure for all R and for different materials raises a question about the definition of opening stress S op in the ΔK eff equation: is it actually S op ? This line of thought is further developed below.
In this paper a distinction is made between the phenomenologically observed stress S op phen , assumed to be the opening stress mentioned in literature and used in crack closure corrections through Eq. (1), and the true physical opening stress S op phys . The choice of the opening stress value S op phen in the ΔS stress range is addressed, and is linked to the true opening stress S op phys using a theoretical model based on multiple linear elastic springs. A FEA simulation is used to generate realistic S op phys values, and together with an energy based model the corresponding S op phen for the ΔS range is obtained. The resulting S op phys derived closure correction resembles existing closure corrections, but is based on a physics approach. It is shown that this new crack closure correction is a function of R, S max , S yield , and a/W.

How experiments support closure corrections
The crack closure corrections, proposed in [10,[15][16][17] are of similar trend but of increasing complexity: see Table 1  nondimensional fitting parameters, to incorporate plane strain or plane stress conditions, or to include a dependency on the flow stress σ 0 . And all these closure corrections are based on curve fitting of a limited amount of measurements. This has several implications for the accuracy of the corrections, relating to applicable R range, measurement techniques, and curve fitting. The measurements on which closure corrections are based are inevitably made over a limited R range. The corrections are then easily extrapolated beyond their original R range, a danger that Schijve [15] warns of. Elber's Eq. (1) is a clear example, since the original R range is known. Extrapolating this equation for negative R ratios gives unrealistic results, and it is generally assumed that it holds correct only for R⩾ − 0.1.
The empirical nature and subsequent curve fitting practice of some closure corrections is illustrated by Overbeeke et al. [24]: a discontinuity is present at R = − 0.5 while there is no physical reason for such a discontinuity to exist.
Furthermore, measuring S op is not straightforward since no measurement device can measure directly at the crack tip. For example, Elber [10] placed a clip gauge system 2 mm behind the crack tip. Such strain gauge measurements or COD measurements give an indication of   [41], propose a crack opening ratio parameter with an elaborate experimental setup to assess fatigue crack closure. Notwithstanding the validity of the method, it remains an indirect, phenomenologically derived measurement of crack opening or closure. The closure corrections shown in Table 1 do not provide scatter bands or error bars on the original measurements. One possible reason is the indirect construction of S op /S max values: these values can be constructed from da/dN versus ΔK eff crack growth curves by having them coincide through introducing a shift in ΔK eff . This shift indicates the ΔS eff needed and hence defines the opening stress S op , or more precise: To give a demonstration of this ambiguity; Schijve [4] used the same data, Ref. [5], as De Koning [16], yet arrived at a slightly different crack closure correction. Using the same data set while arriving at two different equations suggests that there is no consensus on the method how to extract the curve from the data. For this particular data set Schijve gives values of the coefficient of variation CV. The ΔK eff shift as described above is used to generate S op /S max values for six different R values − 1⩽R⩽0.54, which result in CV values ranging from 4.5% to 8.6%. When data point R = − 1 is omitted, the CV range improves towards 4.2% to 6.5%. This illustrates the significant differences at low R values of closure corrections, as shown in Fig. 1. The associated curve fittings are thereby affected too. These variations are partially due to measurement scatter, but are also partially caused by the inadequacy of the closure correction to match the observed ΔK eff shift.
The issues mentioned above affect closure correction measurements, relating to ΔS eff . ΔK eff , using ΔS eff , is an improvement over ΔK. This improvement is not complete either, as the overall validity of the ΔK eff equation is questioned in literature. Kujawski [42] observes that ΔK tends to under-predict and ΔK eff tends to over-predict crack closure effects. Scaling parameters need to be applied to both ΔK and ΔK eff to match observations more closely. Castro et al. [43] mention experimental results which cannot be fully explained by using either ΔK or ΔK eff as similitude parameter. Furthermore, while the incompleteness to describe crack closure of ΔK is noticed, it is not explained why ΔK eff is incomplete too.
The closure correction measurement issues presented here suggest that ΔK eff does not fully or not correctly account for crack closure. There appears to be an inconsistency in the definition of S op : the phenomenologically derived opening stress is likely not the physical opening stress. S op phen can be understood as a virtual stress representing the applied work used for crack growth, in accordance with the ΔK eff approach. It is hypothesized that S op phen is not necessarily equal to S op phys , and this line of thought is further tested as explained in the next sections.

Why S op is phenomenological
Consider Eq. (1); it is often implicitly assumed to hold for all R. Note that for R close to unity, the crack might not close at all, because during the load cycle, the stress never drops below the physical crack opening stress; in other words S min > S op phys . In such a case, in principle there should not be any closure correction. Literature confirms that closure is not measured at high R values; R = 0.7 is the maximum value for which closure has been reported [17,29].
Furthermore, plane strain conditions (mid-thickness) see significantly less plasticity induced closure compared to plane stress conditions (near or at the surface). Observing crack closure at a specimen surface therefore tends to overestimate the amount and influence of closure on the crack growth rate. The combination of geometry and S max results in differences in the transition from plane strain to plane stress conditions during crack growth, which affects the observation of closure in different tests. It follows that S op phys /S max = f(R, S max ): even though closure corrections are normalized by S max , there is still an S max dependency related to the internal stress conditions. This finding is in line with the aforementioned doubts that Kujawski [42] states about the accuracy of ΔK eff .
The choice of S max (for example as ratio of flow stress σ 0 ) at a given R also affects S op phen /S max , as the amount of crack tip plasticity changes. This effect is mentioned since the 1980s by Newman [29,45] and McClung et al. [46,47], but a physical explanation is not given in literature. Fig. 3 is reproduced from Ref. [29], and shows the dependency of S op phen on S max . It clearly demonstrates the significant changes of S op phen / S max versus R for plane stress.
Furthermore, McClung [39] notes that the correlation of ratio of S max over the flow stress σ 0 and S op phen /S max works for CCT specimens only, as the correlation of this ratio and S op phen /S max for other geometries and/or loading conditions is not successful. A new correlation is obtained by introducing a normalized stress intensity parameter ΔK max /K 0 , which appears to work well for small-scale yielding. This example serves to illustrate that improvements or corrections to Eq. (1) are sought using other parameters, rather than looking at the discrepancy between S op phys and S op phen . This results in corrections to the phenomenological description of ΔK, and does not necessarily constitute a physically correct approach.
It is shown that closure corrections, even with fitting parameters or alterations to the ΔK eff method, are not complete in describing crack closure. This paper suggests therefore another approach: an explanation and a solution for the incompleteness can be found in the opening stress itself: S op phen ∕ = S op phys . They are related, but not equal. The true background of S op phen is explained in more detail further on.

An energy equivalent area approach to S op phen
Alderliesten [48] provides a first step to an energy related explanation of crack closure corrections. An energy equivalent area analogy is presented, showing that ΔS eff = S max − S op phen is correlated to the actual cyclic energy applied between S min and S max , through an equivalent area in the stress-strain curve. It is best explained graphically, using the bilinear force-displacement curves in Fig. 4. For a given S max , S min , and an elastic material, the area under the curve between ε min and ε max relates to the cyclic energy ΔU. A rectangular area is spanned by a strain range 0⩽ε⩽ε max and a stress range S op phen ≤ S ≤ S max such that the area is equal to ΔU, and is called the equivalent area. Alderliesten assumed for this example that the crack opens or closes at S = 0, seen as a change of slope at this stress level. In reality the crack opening happens at a positive stress value S op phys , Ref. [10], and the change of slope is actually a gradual nonlinear change over a certain small stress range. Crack tip plasticity effects will smoothen the transition between a closed and an open crack.
For modeling purposes, the stress-strain loading curve is simplified as a bilinear elastic curve with an instant change of slope at a positive nonzero S op phys value. Both simplifications; the bilinear approximation, and only considering the loading curve, do not alter the applicability or trend of the energy equivalent area analogy. Both cases in Fig. 4 are rather similar, but show how the true crack opening stress alters the stress-strain slope and subsequent decrease of ε max .
The difference in absolute U values in Fig. 4 is irrelevant since these are two different cases. The equivalent area approach holds for each individual case separately, and links the ΔU area to the rectangular equivalent area within that particular case. The increase of S op phys thus results in a larger ΔU and equivalent area, which is only possible with a smaller S op phen . This improved analogy obtains S op phen /S max values which closely follow known closure corrections, although the equivalent area has no direct physical meaning. It explains how S op phen can be derived from the actual cyclic energy ΔU, but shows that this value is not equal to S op phys .
The analogy holds for all R. For sufficiently large R close to 1, where S min > S op phys , the energy equivalent area approach shows correctly that the crack stays open during the full cycle: S min > S op phen . The closure corrections however still predict S op phen > S min .
Closure corrections can be described as corrections for closure effects with respect to the ideal bilinear elastic case where closure happens at S op phys = 0, shown in Fig. 5. The stress-strain curve for the ideal bilinear elastic case is similar to case (a) in Fig. 4. Alderliesten [48] demonstrated that the ideal bilinear elastic case can be derived using the energy equivalent area approach. The area under the curve (such as shown in Fig. 4) is a function of R 2 (relates to the stress values), and a function of S op phys (which influences the corresponding strain). Existing closure corrections are effectively correcting this ideal bilinear elastic case, S op phys = 0, for cases where S op phys ∕ = 0, however not by using S op phys but using the virtual S op phen value. The corrections become more pronounced for decreasing R as plasticity and reverse plasticity effects increase. Fig. 4 is a function of R 2 (relates to the stress values), and a function of S op phys (which influences the corresponding strain). For S < 0, a compressive energy is present, but this does not affect crack growth as the crack is closed during this phase of the load cycle. The ideal bilinear elastic case therefore remains constant for R < 0, while the actual total cyclic energy increases as the sum of the tensile and compressive energy. The ideal curve can be described analytically: Correcting the ideal bilinear case for S op phys ∕ = 0 results in S op phen values similar to existing closure corrections. The approach of Alderliesten [48] therefore also holds for S op phys ∕ = 0, providing a link between S op phys and S op phen .   (1) , Newman (2) , Schijve (3) , De Koning (4) , Iwasaki (5) , Kurihara et al. (6) , Overbeeke et al. (7) , Correia et al. (8) , and the ideal bilinear elastic case as described by Alderliesten with closure at S op phys = 0.

A spring analogy to crack closure using S op phys
The physics based energy approach put forward by Alderliesten is extended below with a model of the crack closure effect. Fig. 6 shows the storable energy in a plate with a crack for two different R values: R < 0 and R > 0 with in both cases S op phys > S min . It is a similar schematic stress-strain curve as Fig. 4, but compressive energy is taken into account for negative R. It is assumed that the crack is either fully open or fully closed: there is a distinct region with linear elastic stiffness E 0 (crack closed), and a region with linear elastic stiffness Consider a uniaxially loaded fatigue plate specimen. It is either loaded in tension or compression, or unloaded. With a crack present, it is convenient to model the plate as two springs of different stiffness working in parallel, with an equal displacement constraint. Fig. 7 provides a schematic view of this model. If the plate has developed a crack, the spring system is no longer linear elastic for all S. Below the opening stress S op phys , the total stiffness still equals E 0 . At the opening stress S op phys the stiffness changes to just E 1 (note that the line with slope E 1 does not start at the origin). The change in stiffness is expressed as The cyclic energy, the change in energy during loading, is equal to the area under the force-displacement curve f(ε) of the plate, and indicated with ΔU: More specifically, the cyclic energy ΔU is the sum of the three regions or a part thereof, depending on the value of S min and S op phys : ΔU is the total cyclic energy stored in the plate. The change in stiffness during loading results in increased cyclic energy. The energy equivalent area approach then obtains a lower S op phen compared to the linear elastic case. Contrary to the ideal bilinear elastic curve outlined earlier, it is assumed that the compressive component of the cyclic energy is also involved even though the crack is closed during this part of the cycle: elastic stresses and reverse plasticity occur around the crack tip, influencing crack opening in the next loading phase. Furthermore, this physics based approach does not need fitting parameters as used in several closure corrections such as Newman [17,16]. The use of FEA to obtain S op phys values already includes S yield and a/W effects. This also reduces the need for fatigue tests to gather phenomenological fitting data.

Finite element analysis to obtain S op phys
In order to find true opening stress S op phys , and given the difficulties of determining it experimentally, a finite element simulation approach was chosen. Literature contains many finite element analysis studies investigating crack closure effects. The majority of these are 2D simulations, Refs. [49][50][51][52][53][54][55]39,47,56], using either plane stress or plane strain conditions. Newman [57] notes that since the mid-1980s relatively few 3D FEA studies have been undertaken, Refs. [58][59][60][61]. Kotousov et al., Ref. [62], note that progress in 3D FEA is still well behind that of 2D FEA. In nearly all FEA studies, the crack is instantaneously extended at maximum load, which corresponds well with measurement data. FEA of small crack growth, Ref. [63], might raise questions regarding the mesh size in comparison to the plastic zone size, but apart from that, literature shows that closure can be modeled well with FEA. The FEA analysis discussed below uses instantaneous crack extension at maximum load, and simulates a developed crack and plastic zone, in accordance with FEA studies in literature.
A finite element analysis was performed with SIMULIA Abaqus software on a CCT plane stress plate under constant amplitude (CA) loading. Symmetry conditions apply, therefore only one quarter of the plate was simulated. An infinitely stiff beam was used to model contact along the crack path. Fig. 8 shows the quarter plate and a detail of the mesh around the crack tip. A 2D mesh of quadrilateral element type CPS4R was generated, with enhanced hourglassing control for error reduction. An FEA mesh study showed that a two times finer mesh resulted in similar convergence rates, and similar S op phys values (< 3%). The material was an elastic-plastic model of Al 2024-T3, Ref. [64], of which the stress versus engineering strain curve is given in Fig. 9. This is a monotonic elastic-plastic curve. A cyclic curve may be used as well without changing the validity of the method. Crack opening was defined as contact removal between the node pair directly behind the crack tip. The starter crack length equals a = 10 mm, on a total plate width of 160 mm. FEA simulations were made at different S max values over the range − 1⩽R⩽1 in steps of 0.05. For each R, the model was run for 24 full cycles, with each half cycle divided in 200 equidistant partial loading steps. At maximum load of each cycle, a node on the crack center line was released to simulate crack propagation of one element distance; 0.05 mm. Note that the simulation does not try to mimick realistic increase of the crack growth rate over a significant part of the crack life,  but rather it tries to obtain the opening stress at a crack length which is (nearly) constant. This strenghtens the assumption of constant da/dN (equal to one element length), and constant a/W. The initial crack length is 10.00 mm, the end crack length is 11.20 mm, giving a final crack length over width ratio of a/W = 0.07. Fig. 10 shows details of the loading sequence and the node release.
The instant during the cycle where crack opening occurs is crossreferenced with the applied loading curve, to obtain the corresponding opening stress value S op phys . The FEA results showed that for every R value, within a few cycles the S op phys value had already converged to a steady state value, and the S op phys value of the last cycle was taken as the steady state value for that particular R value. Fig. 11 shows such convergence of the opening stress for R = 0.5: within a few cycles, the crack tip plastic zone is well developed and reaches a steady state. This plastic zone development is the result of starting the FEA at a crack length without any plasticity history from prior cycles.
Katcher [65], Zhang [66], and Newman [49] mention that releasing a node at maximum load is a realistic approach to crack growth simulation, and together with Schijve [15] they also suggest that steady state opening and closing stresses are generally observed to be quite similar which vindicates this approach to find the opening stress to assess closure effects.
The results of the FEA analysis is shown in Fig. 12, for different values of S max . Note that these are S op phys /S max values. Four distinct R regions are observed in each data set. In region 1 and 2 the values of S op phys /S max are linear functions of R, but the slopes are different. It is hypothesized that the more shallow slope of region 1 is a consequence of the force-displacement curve of the plate being influenced more by the closure effect at negative R. The compressive energy increases the reverse plasticity volume, thereby slightly raising the S op phys value for the subsequent cycles. compared to the trend of region 2. Region 3 shows a transition to a crack that never fully closes. The FEA simulation shows that crack opening happens increasingly early with increasing R, until S op phys = S min . This upper limit matches the apparent maximum R limit where closure can be experimentally observed (such as R = 0.7 in Newman's data, Ref. [17]). In region 4 (R ⪆ 0.7) the crack stays open during the full cycle, and as a result the simulation assumes S op phys = S min , meaning that the crack opens immediately at the beginning of the cycle.
Three FEA data sets are shown in Fig. 12. Even though S op /S max values are normalized with respect to S max , there is still an obvious S max dependency. In regions 1 and 2 it holds that for a given R, an increase in S max results in a decrease of S op phys /S max . This trend becomes stronger for smaller R, and is the equivalent of closure corrections moving away (down) from the ideal bilinear elastic curve. The increased S min values for a given R result in more compressive stress, which needs to be overcome before the crack opens. This provides a physical explanation Fig. 9. Stress versus engineering strain curve of Al 2024 T3 material used for the opening stress FEA. Curve based on data from [64]. Dots indicate the points used in the Abaqus material model, which linearly interpolates inbetween them.   for the observations of Newman [29], as shown in Fig. 3, but now explained using FEA. Furthermore, the boundary between region 1 and 2 is observed to stay around R = 0 for all cases as compressive stress only becomes significant for R < 0. The slope of region 1 slightly increases at larger S max , lowering the aforementioned S op phys /S max values. Increasing S max also increases the R value at which the crack stays permanently open (boundary between region 3 and 4).
It should not be attempted to directly compare the FEA opening stress results with other closure corrections from literature, as S op phys ∕ = S op phen . Assuming that the FEA results are representative of real fatigue tests, ΔU values for each point can be obtained by determining the area under the actual force-displacement graph for each corresponding S max and R value. As explained earlier, the force displacement curve is modeled as bilinear, with a change of slope at S = S op phys . The energy equivalent area analogy then reworks the S op phys based ΔU values into S op phen /S max values, which are shown in Fig. 13. Note how these reworked FEA results closely follow the known closure correction trends, and the obtained values are mainly within the spread of various published closure corrections. This implies that the S op phen values found by the known crack closure corrections are a good measure of the equivalent area, but do not constitute the physical opening stress S op phys .
Another observation from the FEA is also seen in literature for different metals and alloys: the R dependency of the effective stress intensity factor ratio. This ratio is defined as: Maljaars et al. [68] have compiled several equations of U SIF for steel and aluminum from literature, Refs. [23,69,24,22,70,3,71], including the respective valid R ranges. Fig. 14 shows (a), a comparison of the effective SIF ratio U SIF for steel reproduced from Ref. [68], and (b), the calculated effective SIF ratio for aluminum 2024-T3 from the FEA data. The trend for the curves is similar, even though the absolute values differ, likely due to the different material parameters considered, such as S yield and E. The flat section of U SIF = 1 near R = 1 is correctly predicted by FEA, showing R values for which the crack is permanently open (S min > S op phys ).

Closure corrections and the finite width correction
Fatigue specimens have a finite width, contrary to most theoretical models, including crack closure models. The finite width causes the crack growth to increase faster compared to the theoretical ideal infinite plate case. As the crack grows in a finite width plate, the remaining cross-section decreases, causing the net-section stress over the crack plane to increase. This effect is usually accounted for in ΔK eff , Eq. (1), by a finite width correction β. There are various analytical expressions for β, often as a function of a and W; the Feddersen equation Ref. [72] is well known. Similar to the closure corrections it holds that these equations can be excellent approximations (of experimental data) but they do not describe the underlying physics. The finite width correction β follows from the change in specimen compliance due to reduced stiffness resulting from the decreasing cross-section. Zhao et al. [73,74] and Alderliesten [75] describe this effect from a physics standpoint in the non-isotropic composite material GLARE, but it holds for isotropic materials as well since it is a consequence of the geometry, not the material. Crack closure is a local phenomenon around the crack tip, and is related to the local stress state and amount of crack tip plasticity. In a finite width specimen, this local stress increases related to the global specimen compliance, for reasons explained above. For a given R and S max , crack closure is thus affected by the ratio a/W. During crack growth a/W and β increase, leading to increased stress around the crack tip, which results in the value of S op phen /S max decreasing. This movement of the closure correction is schematically shown in Fig. 15. It appears that closure corrections from literature implicitly assume S op phen ∕ = f(a), likely because the effect is small for positive R values. The R range of the original closure correction by Elber [10] (Fig. 1) contains mostly positive R values, and presumably the crack length dependency therefore went unnoticed. The spread in closure corrections from literature mentioned earlier can be partly explained by different a/W ratios in the measurements. Since these a/W ratios are not reported, it is difficult to assess the extent of the finite width effect on known closure corrections from literature. To properly observe crack closure by experimental methods such as COD, it is beneficial to have a sufficiently large crack length a and a corresponding large change in specimen compliance upon loading (large da/dN). Because the finite width effect acts on local phenomena such as crack closure, it is therefore not the cause of the difference between S op phen /S max and S op phys /S max which are based on global parameters. It does mean that both S op phen /S max and S op phys /S max are f(R, S max ,S yield ,a/W). The FEA approach of finding S op phys already contains the a/W dependency, removing the need for phenomenological fitting and correction parameters in the proposed model.

Variable amplitude crack growth prediction with energy based closure correction for ΔK eff
Many models for variable amplitude (VA) fatigue prediction make use of the ΔK eff Eq. (1): crack closure models and strip yield models such as ONERA [76], CORPUS [77,78], PREFFAS [79], NASGRO [12], and Refs. [80][81][82][83], based on Dugdale's original work regarding strip yield models, Ref. [84]. In this section it is shown that the energy based closure correction may be used to obtain more accurate values for S op phen from S op phys , on a cycle-to-cycle basis, compared to strip yield models. A typical VA model would predict the S op phen value for each load cycle, taking into account a certain amount of the spectrum history. Schijve [15] states that strip yield models are superior to other types of (global stress based) crack closure models, since the crack geometry around the crack tip is directly modeled. This statement corresponds with the idea that a strip yield model should be based on S op phys , Fig. 12, rather than on S op phen , Fig. 13. However, Matias et al. [85] evaluated the strip yield model of NASGRO [12] for various aircraft loading spectra (VA), and reported that it correlates reasonably only for negative R ratios. Note that for a given S max , positive R values affect the tensile energy terms ΔU a and ΔU b of Eq. (5) which are dependent on S op phys , while negative R values affect only the compressive energy term ΔU c where the crack is closed. The energy approach outlined earlier can therefore be a more suitable candidate for fatigue modeling at any R: it works over the full R range unlike several closure corrections and strip yield methods, and it does not need scaling or fitting parameters as the energy based physics approach is used to obtain S op phys and S op phen . It is thereby also applicable for VA spectrum modeling. We will now discuss the qualitative example of an overload in a CA spectrum introducing crack growth retardation.
Consider a cracked plate loaded with a CA fatigue spectrum at a small positive R with a single overload. A schematic view of the force--displacement curves of the plate at various cycles during the spectrum is given in Fig. 16. The cyclic energy ΔU during the few CA cycles just before the overload can be considered essentially constant. When the single overload cycle occurs, the increase in S max increases ΔU. It also raises the amount of plasticity and therefore S op phys = f(S min , S max ) significantly. While for subsequent cycles the S max value has returned to the CA spectrum level, the crack needs to grow through the enlarged plastic zone generated by the overload, which slowly lowers S op phys from the elevated level back to its original level. Our method correctly predicts this behavior from the physical S op phys instead of the phenomenological S op phen .
Refs. [86,15] report that the opening stress S op phen indeed changes after an overload. Fig. 17 is reproduced from Ref. [86]. From the energy equivalent area approach it follows that S op phys will have a similar behavior. It proves that our theoretical model of the effect of an overload on a CA spectrum is in agreement with what is reported in literature, such as results from COD measurements.
In a truly random VA spectrum, every half cycle sees a new (S min , S max ) pair. S op phys and ΔU change virtually every half cycle, but a steady state is never reached during each half cycle. This means that the S op phys and ΔU are continuously chasing a non-existing CA equilibrium which changes each half cycle. The history of many consecutive cycles (ideally all previous cycles), needs to be taken into account to properly model the behavior of S op phys , ΔU, plasticity, and da/dN, in order to understand and predict VA fatigue crack growth. Amsterdam [87] describes a method that uses a maximum reference stress different from S max to account for VA spectra and pivot points. Pivot points connect multiple power law exponents at different crack length ranges of a Paris type crack growth curve. The altering of the maximum stress results in better power law curve fitting of the crack growth rate, but implicitly leaves S op phys unaffected. Applying the method proposed here, i.e. calculating S op phys for each half cycle, would give a similar outcome to the Amsterdam approach, as it would affect the ΔS eff used in the ΔK eff equation. Our method therefore gives a physical explanation for the ΔK eff value used in the Amsterdam approach.
The shortcomings of the closure corrections and the inconsistency between plasticity effects and the LEFM ΔK eff approach suggest that an energy based approach (S op phys , the elastic spring model, and the energy equivalent area approach) is an improvement over existing methods to properly include closure and plasticity effects. It results in more accurate descriptions of S max , S op phys , and a/W effects on the crack growth behavior, which can improve modeling accuracy for CA and VA fatigue. VA fatigue modeling will likely still require a cycle-by-cycle prediction method to include the (full) loading history effects on S op phys .
This new physics based approach is more realistic than existing closure models, as it is based on the true opening stress. It may be the way forward to increase prediction and modeling accuracy for both CA and VA fatigue crack growth cases.

Conclusions
Ample fatigue crack growth closure corrections exist which result in opening stress S op values used in the ΔK eff approach. Literature shows that these corrections still do not fully account for closure effects, and  Fig. 15. Illustration of the effect of finite width on crack closure. The finite width of a specimen increases the crack growth rate due to increased netsection stress on the remaining crack plane cross-section. As a result, crack tip plasticity increases too, lowering S op phen /S max . For a given S max , the crack closure correction moves down as the crack progresses. that they have a large spread. This led to the hypothesis that there are two distinct S op values: S op phen used in the ΔK eff approach, and S op phys which is the true crack opening stress. The subsequent investigation has shown that this is true, and resulted in the following conclusions: 1. Many known closure corrections, Refs. [10,15,17,16] result in different corrections for the same phenomenon, for unknown reasons. Our research suggests that existing closure corrections do not properly take into account the physical opening stress S op phys . 2. The same known closure corrections, Refs. [10,15,17,16], contain empirically determined correction factors, and different corrections result in different S op /S max values based on S op phen for the same load cycle. Our results suggest that these observations may be explained by a failure to correctly account for the effect of S max and a/W on the physical opening stress S op phys . Thus the corresponding S op phen needs to be corrected to make up for this. By applying the equivalent energy approach, Ref. [48], a value of S op phen can be determined based on the correct S op phys , as found via FEA. The S op phen values found in this way matches the S op phen found via previously known correction methods, but without relying on empirical correction factors. Thus our new method is more physically realistic, and potentially requires less experimental calibration than existing methods. 3. FEA shows that S op phys follows four distinct regions over the full R range. From low to high R, these are: tensile-compressive loading with closure, tensile-tensile loading with closure, transition to an always open crack, and an always open crack. In the last case, it can be assumed that S op phys = S min . Using the energy equivalent area analogy to obtain S op phen values for the full R range, the FEA results show close agreement with known correction curves, especially for R⩾0. This confirms that S op phen in S op phen /S max is not the true opening stress. 4. FEA confirmed that it at least holds that S op phys = f(R, S max ) and that the derived S op phen /S max = f(R, S max ) too.
During crack growth, the local net-section stress in the crack plane Fig. 16. Schematic view of ΔU (top) during a CA spectrum (bottom) with a single overload. Four specific cycles in the spectrum are shown, from left to right: a CA spectrum before an overload, a single overload, the first CA cycle after overload, a CA cycle many cycles after the overload. Fig. 17. CA spectrum with a single overload, associated crack growth retardation, and associated S op phen measured with a COD technique. Figure reproduced from Schijve [86].
increases due to the finite width effect. This does not affect S max , but it does affect the local phenomenon which is crack closure. During crack growth at constant R and S max , the S op phen /S max or S op phys /S max values slowly decrease because of the decreasing cross-section area. This dependency on a/W is not mentioned in literature, but may partly explain the spread of closure corrections. Furthermore, several closure corrections include a fitting parameter to match measurement data, which tends to be dependent on R. In literature it is noted that S op phen also depends on the ratio between S max and σ 0 , where σ 0 is related to S yield .
These observations combined suggest that S op phen = f(R, S max , S yield , a/W) and S op phys = f(R, S max , S yield , a/W).
Inclusion of the energy approach into the ΔK eff equation may improve the accuracy over known models. S op phys might be found through FEA or strip yield models. For VA fatigue a cycle-by-cycle integration of the loading history might still be necessary to keep track of changes in S op phys and subsequent S op phen .

Declaration of Competing 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.