[BUG] CS fatigue model has a couple of bugs.

[CoronelBuendia]

Describe the bug

In making a similar model for another code, I came across a couple of errors in the CS fatigue model. As far as I can tell, only one of them will have been encountered in practice, and fortunately (for the random test cases I have checked) it appears to be on the conservative side.

Steps to reproduce

In the EFDA report and in the original reference (J. C. Newman, Jr. and I. S. Raju, https://core.ac.uk/download/pdf/42849129.pdf), equation 19 reads:

which corresponds to PROCESS L176 in cs_fatigue.py:

g = 1.0e0 + (0.1e0 + 0.35e0 * a_c**2.0e0) * (1.0e0 - sin_phi) ** 2.0e0

which should read:

g = 1.0e0 + (0.1e0 + 0.35e0 * a_t**2.0e0) * (1.0e0 - sin_phi) ** 2.0e0

equation 36 in Newman and Raju reads:

which corresponds to PROCESS L208:

+ (0.55e0 - 0.72e0 * c_a * 0.75e0 + 0.14e0 * c_a * 1.5e0) * a_t_2  # G22

which should read:

+ (0.55e0 - 0.72e0 * c_a ** 0.75e0 + 0.14e0 * c_a * 1.5e0) * a_t_2  # G22

Finally, AFAICT, although not explicitly mentioned in the original reference, values of stress are in MPa, not Pa. This is corrected for at the start of CsFatigue.ncycle for the hoop and residual stress, which are input in Pa and converted to MPa. What I am not sure about is how the fracture toughness is passed in. If in the IN.DAT (or whatever) the fracture toughness is specified in MPa.m^(1/2), then this is probably OK. But if it is specified in Pa.m^(1/2) then I reckon there should also be a conversion to MPa.m(1/2) in ncycle.

These last two points will not affect the result in a normal run, as the semi-elliptical crack is basically hard-coded down the a < c branch meaning the 2nd point will never be hit, and the last point appears not to be one of the break conditions of the while loop for inputs typical to the problem. (It also may not be an issue at all, depending on the input units of fracture toughness).

From my alternative implementation elsewhere, I reckon that the first error is actually more pessimistic in terms of number of cycles. For relatively random but realistic combinations of inputs, I get between 1.5 and 12.5 % more pulses out of ncycle with the corrected equation 19.

Expected behaviour

Evidence

Environment

• OS:
• Version:

Additional context

[ajpearcey]

The units of the fracture toughness through the code and the IN.DAT are MPa.m^1/2. So this should be okay, if a bit inconsistent with how this handled elsewhere.