Computing the stress intensity factor range for fatigue crack growth testing at 20 kHz

Inertia and damping influence the values of the stress intensity factors (SIFs) at high‐frequency loading and they must be included in computations. In the present study, different dynamic simulation procedures were carried out for two types of specimen geometries and the achieved SIF values were compared. Fast computation procedures such as harmonic modal analysis and direct steady‐state analysis were compared to the computationally expensive transient dynamic analysis. Two different methods for calculating the SIF, the J‐integral and the crack tip opening displacement (CTOD) methods, were applied and compared and the results showed a near perfect agreement in calculation of the mode I SIF. The Rayleigh damping model was introduced into the dynamic computation to investigate its effect and the results revealed a clear effect on the SIF at 20 kHz frequency. The fast direct steady‐state analysis showed good agreement to both harmonic modal and transient analysis with the different damping values used and is, after this study, the recommended procedure.


INTRODUCTION
In crack growth rate testing, the results are appropriately presented in a crack growth curve according to the Paris law, Equation (1).
For that, a relation between the stress intensity factor range (K) and the crack growth rate (da/dN) is required.While the crack growth rate is measured experimentally, the stress intensity factor (SIF) is computed analytically or numerically for a varying crack length.Crack growth testing using the ultrasonic fatigue testing system is carried out at 20 kHz load frequency, Figure 1.The 20 kHz load frequency enables testing for average crack growth rate in the range of 10 −8 − 10 −12 m∕cycle.The testing procedure is according to the ASTM E647 standard where the crack is let to propagate F I G U R E 1 Schematic figure displaying the main parts of the ultrasonic fatigue testing system. in small increments (0.1 − 0.5 mm) at varying ΔK-values.After the crack initiation, the ΔK-value is decreased each increment until crack arrest (i.e., reaching an increment where no crack growth occurs during 10 −7 cycles), and then increased each increment until final failure.The threshold value (ΔK th ) is estimated from the ΔK-values at the crack arrest and the Paris law curve is drawn with the data points from the ΔK-increasing part of the testing.Some experimental results are presented in Reference1 where three high strength automotive steels were tested at different load ratios.
A relationship between the controlling parameter of the testing system, that is, displacement at the top of the specimen, and the mode I SIF at the crack tip in the mid-section of the specimen is required.The 20 kHz fatigue crack growth (FCG) specimens have, during uniaxial oscillation, their maximum displacement at the ends of the specimens while at the mid-section the displacement is zero.The zero displacement at the mid-section of the specimens is highly convenient for the crack length measurements.
There are different models for calculating the SIF.In elastic-plastic fracture mechanics, mainly two models are recommended, the J-integral and the Crack Tip Opening Displacement (CTOD)-methods. 2,3The J-integral method is based on the accumulated strain energy while the CTOD method is based on, namely, the crack tip opening displacements close to the crack tip.
However, when calculating the SIF for high frequencies, for example, 20 kHz, static calculations are no longer sufficient.In static simulations, as used in References 4-7 for computation of SIF at 20 kHz frequency, the load is ramped up to a preselected level and held constant.No oscillation is applied and the dynamic conditions, that is, inertia forces and damping effects are not considered. 8Hence, dynamic analysis and simulations are highly recommended for more accurate results. 2,3Anyway, with both static and dynamic analysis, the SIF can be calculated using both the J-integral and the CTOD methods.
There are different computation procedures revealing the dynamic response of a system caused by harmonic excitation.Harmonic modal analysis, transient dynamic analysis and direct steady-state dynamic analysis are the different simulation procedures used and compared in this study.

Harmonic modal analysis
Harmonic modal analysis is, on the other hand, a very computationally efficient procedure.The analysis calculates the shape mode of the model at specific resonance frequencies in specific modes, that is, modal analysis, 3,8,9 and the computation procedure considers elastic material behavior and does not take the damping into account.However, this simulation method can potentially eliminate large computational time consumption and hence is investigated in this study.The J-integral is not possible to compute using the harmonic modal analysis, hence in this case the SIF is computed with the CTOD method.

Transient dynamic analysis
Discretisation of time yields a set of linear equations describing the response within each time increment.Such solution starts with an aperiodic oscillation and converges to a single frequency steady-state oscillation.The time, or number of cycles required to pass the non-periodic transients and reach steady oscillation depends on material damping, with increased sensitivity at low damping coefficient. 8The normalized SIF range variation during the first 20 simulated cycles is larger using smaller damping as demonstrated using three different values of damping selected in the range of the present interest, Figure 2.For these simulations, the time increment used was 10 −7 s giving 500 data points per cycle.
In a previous study, 10 a transient dynamic analysis has been applied where the J-integral and thereby the SIF at a crack in a fatigue crack growth (FCG) bar specimen has been computed.In References 11-13 the same simulation procedure (transient dynamic) has been used to conduct stress-strain analysis in 20 kHz fatigue strength hourglass specimens.Such analysis is computationally inefficient, especially when small values of material damping are used, and a large number of cycles are required to pass the transients and reach the steady state oscillation.In this study, this procedure is compared to more computationally efficient alternatives.

Direct steady-state dynamic analysis
Direct steady-state dynamic analysis is a linear perturbation procedure used to calculate the steady-state response of a system under a harmonic excitation at a prescribed range of frequencies.Conveniently, the frequency range is chosen F I G U R E 3 A discrete system with two mass points, two springs and two dashpots under the action of a periodic force F = f (t).
to cover the desired resonance frequency computed in a preceded harmonic modal analysis (modal frequency extraction analysis).The response is calculated in terms of physical degrees of freedom using the mass, damping and stiffness matrixes.The commercial FEM software Abaqus offers the direct steady-state dynamic analysis, which is appropriate when frequency dependent damping (i.e., Rayleigh damping) is present. 8he direct steady-state method is explained in Reference 14 where a simple model with two mass points m, two springs c, two dampers r and an applied periodic force F, Figure 3, has been analyzed.The steady-state solution of the same has also been computed by transient dynamic analysis and with a perfect match to the direct steady-state analysis.
The FEM software Abaqus yields the node displacements as a complex number where the imaginary part accounts for the phase shift due to damping.When using damping in direct steady-state analysis, as in the transient analysis, a phase shift exists in the oscillation of the model and increases as the distance to the region of the applied displacement increases.The node displacements are described by where |u| is the magnitude,  is the angular frequency,  is the phase angel and the real and imaginary part are In direct steady-state analysis the J-integral is not computed and the SIF is calculated using the CTOD method.Hence, the opening of the crack tip (i.e., distance between the upper and lower surfaces of the opened crack) is of interest.However, due to the present phase shift, the displacements of the upper u upper and lower u lower surfaces of the crack are not at their maximum at the same time.Hence, the maximum crack opening u opening is calculated as follows: where The direct steady-state analysis has been used in References 7,15-17 to compute the SIF for crack growth rate testing of thin sheet specimens of martensitic-bainitic hot rolled and complex-phase ferritic-martensitic steels at 20 kHz.However, the effect of damping was not taken into account.

Objectives
The objectives of the present study are to acquire more knowledge regarding the effects of specimen geometries (bar and sheet), methods of calculation (J-integral and CTOD) and simulation procedure on the SIF.

Experimental setup
The experimental setup consists of a generator transforming a 60 Hz voltage signal into a 20 kHz electrical sinusoidal signal, a piezoelectric transducer transforming the 20 kHz electric signal into a 20 kHz mechanical oscillation, a horn to magnify the displacement amplitude and the specimen, see Figure 3.The dimensions of both the horn and the specimen are carefully calculated for an individual uniaxial resonance frequency within the limits of the ultrasonic fatigue testing system, that is, 20 ± 0.5 kHz.Two specimen geometries, bar and sheet, according to, 3 were considered in this study, Figure 4.The mechanical characteristics used in the simulations are typical for the materials of the specimen and the horn: E d = 210 GPa,  = 0.3 and  = 7800 kg/m 3 for the steel specimen and E d = 110 GPa,  = 0.3 and  = 4420 kg/m 3 for the titanium horn.To implement the damping into the simulations, according to Equation ( 8), the Rayleigh mass-and stiffness-proportional damping coefficients  R and  R are required.
In this study the mass-proportional coefficient is kept constant at  R = 0.1 due to its low effect in metals at 20 kHz.While the stiffness-proportional coefficient is varied within the range of interest, as explained later.
When using complicated specimen geometry (i.e., varying cross-sections, notches, etc.) FEM software are conveniently used in particular when including the dynamic effects (e.g., inertia and damping) into the calculations.In this study, the FEM software Abaqus 2019 was used to compute the SIF for the two different specimen geometries using both the J-integral and CTOD methods.

Boundary conditions and mesh
The boundary conditions in the FEM models were used as follows.For the simulations a half symmetry 3D-model was used for both specimen geometries, and the horn was included, see Figure 5, with a tie constraint between the bottom

F I G U R E 5
The 1 / 2 symmetry 3D-models of the (A) bar, and (B) sheet specimens, and including the horn.
surface of the horn and the top surface of the specimen.As mentioned earlier, a previous study 10 has recommended method where the whole load train has been included (specimen, horn and oscillator) into the simulated 3D-model.However, for the sake of minimizing the number of degrees of freedom, the oscillator was removed from the present simulation procedure as it is a comparative study.
A hexagonal focused mesh with 8-node linear bricks (C3D8R) around the crack tip and 6-node linear triangular prism (C3D6) at the crack tip were used for the bar specimen.For the sheet specimen, a hexagonal focused mesh with 20-node quadratic bricks (C3D20R) around the crack tip and 15-node quadratic triangular prism (C3D15) at the crack tip were used.The rest of the model were meshed using 10-node quadratic tetrahedron elements (C3D10) for both specimen types in all simulations (Figures 5 and 6).Tie constraint couples the regions with different meshes.Using linear or quadratic elements at and around the crack tip shows negligible effect during these simulations.Two different sets of crack lengths were chosen for the bar and sheet specimens to be simulated for every variation in the simulation model.
The focused mesh around the crack tip contains eight contours.In order to achieve convergence in J-integral and crack surface displacement, the eight nodes on the crack lip, at the middle in the thickness direction and at the eight contours, Figure 6, were used to extract both the J-integral values and the displacement values used in the CTOD method.
In the transient analysis, a harmonic excitation was applied as a sinusoidal displacement (u 0 = 3.8 μm) to the top surface of the horn at a crack length specific frequency (experimentally measured or computed by harmonic modal analysis).The horn has a magnification factor of 2.5 increasing the sinusoidal displacement to approximately U 0 = 10 μm at the top of the specimen.
The simulations were run for 0.001 s enough for 20 cycles.Physical, elastic and damping properties were input to the analysis.Default parameters for Hilber-Hughes-Taylor time integration were used ( HHT = −0.05, HHT = 0.28, and  HHT = 0.55).This insures diminishing numerical damping.The SIF in all three modes (K I , K II and K III ) were automatically computed by Abaqus using the J-integral method.For the CTOD method however, the displacements in

F I G U R E 7
Experimentally measured frequencies and modal frequencies computed by the harmonic modal analysis plotted against a/w for bar and sheet specimen.
all three directions (u x , u y and u z ) at the eight mentioned nodes (red nodes in Figure 6) were plotted and used to calculate and extrapolate K I , K II and K III at the crack tip.
For the harmonic modal analysis, only physical (density) and elastic properties (elastic modulus and Poisson's ratio) were used since the damping is not considered in such analysis.The outcome of this analysis is the modal shape of the whole system (specimen, horn and oscillator) at the found resonance frequency.The modal shape reveals the normalized displacement of all nodes in the model.The displacements at the specimen top surface and the nodes at the crack lips are extracted and used to compute the SIF with the CTOD method.
The direct steady-state dynamic analysis requires, as the transient analysis does, physical, elastic and damping properties.The range of frequency was specified to cover the desired frequency, and a displacement amplitude was applied to the horn top surface, u 0 = 3.8 μm.The results of this analysis reveal, similarly to the harmonic modal analysis, the displacements of all nodes in the model.
The resonance frequency of the system will decrease as the crack length increases due to change in specimen compliance.For the direct steady-state analysis, the crack length specific harmonic modal frequency was used.The harmonic modal frequency was computed by harmonic modal analysis where specimen and horn were included into the FEM model.For the transient dynamic analysis, both the harmonic modal frequencies and experimentally measured frequencies at the individual crack lengths were used for the sake of comparison.The experimental frequency is simply the operational frequency displayed by the software during testing.For the bar specimen, the experimental frequencies were measured in a previous study 10 while the experimental frequencies for the sheet specimen were provided by collaborative researchers at Université Paris Ouest (UPO). 7The different sets of frequencies for both the bar and the sheet specimens are presented in Figure 7.

Computation of the SIF
The SIF for a varying crack length is achieved by interpolation of computed K for several specific crack lengths.In, 10 in relation to static simulation, a clear increase in the computed K was shown when using a dynamic simulation including the frequency and damping effect.For elastic-plastic materials with non-linear material deformation in a large region around the crack tip, the linear elastic fracture mechanics (LEFM) is no longer valid.However, the plastic deformation is taken into account in the elastic-plastic fracture mechanics (EPFM) model 2 where the crack tip conditions are described by the J-integral or the CTOD parameters.
The experimental testing at 20 kHz with the ultrasonic fatigue system is a displacement controlled testing procedure run according to Equation (9).The displacement of the top surface of the specimen U 0 is the controlled amplitude during testing. 3 where E d and  are the dynamic elastic Young's modulus and Poisson's ratio, respectively.a is the crack length and the function f(a/w) is the dimensionless shape function to be calibrated.From Equation ( 9) it is seen that the effects of U 0 , E d and  on the SIF are simple and clear.However, there are more material properties involved in this computation, for example, density and damping.These two properties used in the simulations calibrating the shape function f(a/w) and are hence embedded in it.
The Rayleigh damping model, Equation (10), is implemented in the FEM software Abaqus and was conveniently used in all simulations. 18

𝜉
where  is the angular frequency and  R and  R are the mass-and stiffness proportional damping coefficients, respectively.At high frequencies (e.g., 20 kHz), the effect of the mass proportional damping is insignificant, and the stiffness proportional damping becomes the dominant part of the Rayleigh damping, as seen in Equation (10).
For the sake of comparison, the computed K I presented in the graphs in this paper with the different simulation procedures are all corresponding to when U 0 = 1 .

2.4
Calculating the SIF using the J-integral method The J-integral is equivalent to the energy release rate in a body containing a crack.It relates the change in the potential energy to the crack growth, Equation (11). 2 where A is the crack area and  is the potential energy defined as: where U is the strain energy and P is the work done by external forces.The J-integral for an arbitrary counter-clockwise path  around the crack tip is written as: where W is the strain energy density given by: where  ij and  ij are the stress and strain tensors, respectively.T i (Equation 15) and u i are the components of the traction vector and the displacement vector, respectively, and ds is a length increment along the path  .
where n j are the components of the unit vector normal to the path  .The SIF's are extracted from the computed J-integral with Equation 16: where and B is the pre-logarithmic energy factor matrix, the diagonal for homogeneous and isotropic materials. 8

Calculating the SIF using the CTOD-method
The crack tip SIF is assumed proportional to the crack tip blunting and hence it can be estimated by the measurement of the displacement of the crack surfaces, that is, the CTOD-parameter, see Figure 8.The SIF's in mode I, mode II and mode III fracture at the crack lip at a distance r from the crack tip are computed according to Equations ( 18)-(20). 2 F I G U R E 8 Crack tip opening with (A) crack tip blunting, and (B) the effective crack length including the Irwin plastic zone correction. 2 F I G U R E 9 Linear extrapolation of computed ΔK derived from CTOD-method. 7,19ere G is the shear modulus and factor  = (3-4 ) for plane strain and 3D-models. 2Here the displacement u y is half the crack tip opening u opening .The CTOD's (u y , u x and u z ) are measured at the eight nodes near the crack tip (see Figure 6) and the SIF's at the crack tip are estimated by linear extrapolation, as illustrated in Figure 9.

RESULTS
Multiple series of simulations using the different simulation procedures, computation methods and frequency models were conducted on both specimen geometries.The results show a clear agreement between the two different computation models, the J-integral and the CTOD.Both the ΔK by the J-integral and ΔK by the CTOD method are acquired from the exact same simulations, transient dynamic simulations using experimentally measured frequencies.The SIF in the bar specimen increases significantly with crack length a/w, while in the sheet specimen the increase is not as large, Figure 10.All the same, the stress intensities computed by J-integral and CTOD methods agree perfectly.The static analysis results in significantly lower stress intensities than the dynamic analysis, Figure 11.However, the agreement between the J-integral and the CTOD methods was found even when computing by static simulations.
Regarding the difference between the static and transient dynamic analysis, it is in agreement with results achieved in previous work 10 using the J-integral method where the static simulations showed K I -values approximately 30% lower than dynamic transient simulations.In this work the same was computed with the CTOD method and the agreement between the two methods is again clear.
The results showed a clear difference between the SIF's computed with transient dynamic analysis with Rayleigh damping coefficients  R = 0.1 and  R = 10 −6 (according to earlier estimates 13,20 ) and the SIF's computed with undamped harmonic modal analysis, using the same 3D-model and boundary conditions.The two transient dynamic analysis using F I G U R E 10 K I -values normalized to U plotted against a/w (where a is the crack length and w is the specimen width) using J-integral and CTOD methods and bar and sheet specimen geometries.Transient dynamic analysis and damping values  R = 0.1 and  R = 10 −6 was used.two different frequency sets, experimentally measured and harmonic frequencies obtained by harmonic modal analysis, showed near perfect agreement.The difference between the two types of analysis, undamped harmonic modal and damped transient dynamic analysis, had opposite effect on the two specimen types.The SIF's from the harmonic modal analysis appears higher than from the transient analysis for the bar specimen but lower for the sheet specimen, Figure 12.

F I U E K
The SIF's were also computed using the direct steady-state simulation procedure.Here three different values of damping were used for both the steel and the titanium, undamped ( R = 0 and  R = 0), an experimentally measured values ( R = 0.1 and  R = 1.3 * 10 −9 ) measured with the impulse excitation technique (IET) 21 and the same values as used in the transient analysis ( R = 0.1 and  R = 10 −6 ).The effect of the small damping value ( R = 0.1 and  R = 1.3 * 10 −9 ) proofed to be negligible as the computed K I -values coincides with those computed with the undamped ( R = 0 and  R = 0) analysis.The higher damping value ( R = 0.1 and  R = 10 −6 ) however, showed a clear effect and the computed K I -values were in good agreement with those from the transient analysis (with the same damping values).The damping in the titanium horn were varied similarly to the damping in the steel specimen.The mixed mode SIF's, corresponding to U 0 = 1 μm specimen top displacement, are plotted in Figure 13.K I , K II and K III are all computed with the two different methods, J-integral and CTOD, for a series of crack lengths.The simulation procedure used was transient dynamic analysis and Rayleigh damping coefficients  R = 0.1 and  R = 10 −6 .The values computed using the J-integral method showed that for relatively short cracks (up to 20% of the specimen width) K II and K III were essentially non-existent.However, as the crack grows longer K II and K III starts to elevate to significantly high values (K III to higher values than K II ).
The CTOD method showed identical results regarding the K I and K II values.However, since there is no displacement in the mode III direction (i.e., u z = 0), the CTOD method (Equation 20) gives K III = 0.

DISCUSSION
The aim of this study is to acquire more knowledge on the different methods and simulation procedures used, as well as the effect of damping, on the computation of the SIF in cracked specimens loaded at 20 kHz frequency.The computation of SIF at 20 kHz is required in order to conduct crack growth rate testing at 20 kHz using the ultrasonic fatigue testing system.As described in Reference 10 the load train in such system consists of an oscillator, a magnifying horn and a specimen.The system vibrates at the resonance frequency of the whole load train (20 ± 0.5 kHz) which is decreasing with the growing crack due to loss in stiffness.It has been concluded that in order to correctly predict the drop in the resonance frequency by FEM, both the oscillator and the horn needs to be included.Also, the crack surface inter-penetration needed to be accounted for by calculation the effective natural frequency described in Reference 10 However, using different frequency inputs has very little effect on the computed SIF's, as demonstrated by the agreement between the (Transient -exp.freq.) and the (Transient -harm.freq.)curves in Figure 12.For the sake of simplicity, and since this is a comparative study where only the relationship between the SIF and the displacement of the specimen top surface (U 0 ) is of interest, only two different frequency inputs were used in this study and the FEM models only contains the specimen and the horn.The harmonic frequencies in Figure 7 are acquired from a FEM model including the specimen and the horn only, and with crack surface inter-penetration allowed, hence the deviation from the experimentally measured frequencies.
The obvious difference in the computed K I -values between the bar and sheet specimens (see Figure 12) is due to the higher dynamic stiffness in the thicker bar specimen.The lower dynamic stiffness of the sheet specimen results in lower stress intensities when loaded with the same displacement amplitude.
According to Anderson, 2 elastic-plastic fracture mechanics (EPFM) is required when large plastic zone are surrounding the crack tip.Two different methods, the J-integral and the CTOD, are recommended for the computation of the SIF in elastic-plastic materials.Both models were used in this study and the results show a clear agreement between the two.The J-integral is more practical and is here recommended when using transient dynamic analysis since the SIF values are computed directly by the FEM software and are easily extracted.The CTOD-method is a bit more tedious to process.The displacements of the 16 crack lip nodes (red nodes in Figure 6) are extracted and then extrapolated (according to Figure 9) to calculate the SIF's at the crack tip.However, it is not possible to acquire J-integral values from a harmonic modal or a direct steady-state analysis, hence the CTOD method is the recommended option when such simulation procedures are used.
The required fine mesh of the 3D-models used in this study yields a large number of degrees of freedom (DOF), 1 ± 0.2 * 10 6 DOF's, and increases the computation effort for the simulations.Consequently, running transient dynamic analysis for more than 20 cycles in many series was considered too computationally expensive.It makes it difficult to extract qualitative results when using very small damping values, that is, when the Rayleigh stiffness proportional damping coefficient  R < 10 −6 .As shown in Figure 1, the transient simulations reach a near steady-state oscillation during the first 20 cycles for  R ≥ 10 −6 while for the  R = 10 −8 significant distortion is present throughout.
The average of the last 8 cycles (where the oscillation starts to stabilize) from the  R = 10 −6 transient analysis is compared to the direct steady-state analysis with the same damping values, Figure 12.The comparison shows good agreement for both specimen shapes.The  R = 0 and the  R = 1.3 * 10 −9 analysis are, as expected, in good agreement to the undamped harmonic modal analysis.The result of the harmonic modal analysis in this study were similar to results where the same analysis type and sheet specimen shape have been used. 7n Figure 12, a noticeable difference between the K I -values using experimentally measured frequencies and harmonic modal frequencies at larger crack lengths is observed.This is in accordance with the increasing gap between the two frequency sets with increasing crack length, see Figure 7.However, the difference in the K I -values is small and hence indicates a low frequency effect.
Opposite to the bar specimen, the sheet specimen increased in K I when  R = 10 −6 was used compared to when the lower  R = 0 or  R = 1.3 * 10 −9 damping values were used.However, it is clear that the  R = 10 −6 damping value effects the results, and that the effect varies with the crack length (i.e., a/w).
The comparison of K I -values calculated with different damping values are based on finding the maximum K I -values during one load cycle.But as damping introduces a phase shift, the maximum K I -values will occur at different phase angles depending on the damping condition and increases with the distance from the applied displacement region.Therefore, the phase shift effect by damping was considered in both the transient and direct steady-state analyses and evaluation of the results.Examples and discussion on the matter follows below.
The applied harmonic displacement, the displacement at the specimen top and the K I -value at the specimen mid-section crack tip, are displayed during cycle 20 of the transient analysis, Figure 14.An evident phase shift in both specimen types can be seen, and is larger as the distance to the applied displacement region.It is also clear that the effect of the damping on the phase shift is more significant on the sheet specimen.
Thus, when extracting the relationship between the K I and the specimen top displacement from the transient analysis, the maximum peak value is used from both the K I and specimen top displacement curves in Figure 14.
In the direct steady-state analysis, however, it is a bit trickier to consider the phase shift.Here, the real and imaginary part of the displacements of the eight pair of nodes at the crack lips close to the tip, that is, red nodes in Figure 6, are extracted from the FEM software to calculate the maximum crack tip opening according to Equations ( 5)- (7).Half the distance between the node pair (u opening /2) is used to calculate K I with the CTOD-method at the eight contours.The calculated K I -values are then extrapolated to the crack tip according to Figure 9.

(A) (B)
F I U R E 14 U top displacement), displacement of the top surface of the horn (i.e., applied harmonic displacement) and K I normalized to U 0 plotted against time for (A) bar and (B) sheet specimen at (a/w = 0.33).Transient dynamic analysis, J-integral method and  R = 10 −6 damping were used.
TA B L E 1 Summary of K I -values normalized to U 0 (MPa √ m/μm) at a/w = 0.33 computed with the different simulation procedures, and using the CTOD-method.

Bar
Sheet The results imply that when using the same argument during the evaluation of the results, the direct steady-state analysis yields the same results as a converged transient analysis.In the present study, the maximum K I -values were extracted and normalized against the maximum displacement value at the specimen top under the same cycle.The same procedure was carried out when evaluating the results from both transient and direct steady-state simulation procedures, hence the perfect match displayed in Figure 12.A good agreement between the two different computation procedures have been proofed earlier for a very simple model of mechanical vibrations of the stator in an electrical machine. 14The size of the model in the present study was scaled up and included a vast amount of DOF's, still, two computation procedures are in good agreement.
The multi fracture mode analysis (combining K I , K II and K III ) revealed partly expected and partly unexplainable results.The analysis to investigate the K II and K computation by the two different methods, J-integral and CTOD.The computation of the mode I and mode II SIF's yielded the same results using the two different methods, as seen in Figure 13.In mode III SIF, however, the two methods ended up with mismatching results.Given the load mode and boundary conditions used, any crack tip opening displacements in the mode III direction are not expected.Consistent to this, it is confirmed by the CTOD method where the extracted u z magnitudes (Equation 20) are zero.Contradictory, the J-integral method resulted in non-zero values of the K III .In fact, they are significantly high, even higher than the K II -values.A possible explanation to this error is an anomaly in the computational procedure of the FEM software.Further investigation is required but remains outside the scope of this paper.
A summary of K I -values corresponding to U 0 = 1 μm at a/w = 0.33 computed with the different simulation procedures carried out in the present study, are presented in Table 1.Clearly, the direct steady-state analysis yields consistent results when comparing to the undamped harmonic modal and the damped transient analysis.
A method for computing the SIF at 20 kHz loading frequency was presented in the previous study, 10 and proposed a transient dynamic analysis of the whole ultrasonic system load train (specimen, horn and oscillator) at the effective natural frequency  0 , Equation (21), introduced by Chati et al. 22 taking into account the effect of a breathing crack on the eigenfrequency.
where  1 and  2 are the eigenfrequencies of an un-cracked and cracked specimen, respectively.Presently, the method is updated by recommending the direct steady-state dynamic analysis for the computation of the SIF, instead of the computationally expensive transient analysis.The remainder of the recommendations are the same.
As previously mentioned, the direct steady-state analysis does not include computation of the J-integral.Hence, the SIF computation is performed with the CTOD-method which proofed tedious and time consuming during the work of this study.However, there is room for development of the computer software to improve the data evaluation process required for the CTOD-method.

CONCLUSIONS
The SIF for fatigue crack growth testing at 20 kHz loading frequency was successfully computed with several different simulation procedures; harmonic modal, transient dynamic and direct steady-state dynamic analysis.Furthermore, two different methods for evaluation of the SIF, J-integral and CTOD method, were used.The existing recommendations for the SIF computation at 20 kHz was improved and the required computational cost was reduced significantly as of the results of this study.
1.The two stress intensity computational methods, J-integral and CTOD, showed perfectly good agreement in computing the K I in both specimen types and with different simulation procedures.2. Implementing damping into the dynamic simulations clearly affects the results.The effect is opposite in the bar and sheet specimen types.A decrease in the computed K I -value is noticed when using damping in the bar specimen and an increase in the sheet specimen.3. The transient dynamic, direct steady-state dynamic and harmonic modal analysis showed matching results with the different damping values used.4. The direct steady-state dynamic analysis procedure is hereinafter recommended as the most effective procedure computing the SIF at 20 kHz loading frequency with the damping taken into consideration.

NOMENCLATURE
μm) during the first 20 cycles in transient dynamic analysis using three different values of the Rayleigh damping stiffness proportional coefficient.Calculated in the present study with a/w = 0.29 and 20 kHz, for (A) bar specimen, and (B) sheet specimen.

F I G U R E 6
Focused hexagonal mesh around the crack tip containing eight contours, (A) closed, and (B) open crack in the sheet specimen.
I -values normalized to U plotted against a/w (where a is the crack length and w is the specimen width) using J-integral and CTOD methods in static and transient dynamic analysis of the bar specimen.Damping values  R = 0.1 and  R = 10 −6 was used for the transient dynamic analysis.FI G U R E 12 K I -values normalized to U 0 plotted against a/w for bar and sheet specimens computed by the CTOD method with transient, harmonic modal and direct steady-state (DSS) analysis with different damping values.The two curves yielded from the ( R = 0.1 and  R = 1.3*10 −9 ) and the ( R = 0 and  R = 0) DSS simulations are overlapping, for both the bar and the sheet specimens.
13 K I , K II , and K III normalized to U 0 computed with the two methods, J-integral and CTOD using transient dynamic analysis with  R = 0.1 and  R = 10 −6 .(A) Bar and (B) sheet specimen.