-: 0:Colorization: line numbers: hotness: > 50% > 20% > 10%
-: 0:Source:fatigue2.f90
-: 0:Graph:a-fatigue2.gcno
-: 0:Data:a-fatigue2.gcda
-: 0:Runs:1
-: 1:module free_input
-: 2:!
-: 3:! Copyright (C) 1997 by Quetzal Computational Associates, Inc.
-: 4:! Address any questions about this copyright notice to:
-: 5:!
-: 6:! Dr. John K. Prentice
-: 7:! Quetzal Computational Associates, Incorporated
-: 8:! 3455 Main Avenue, Suite 4
-: 9:! Durango, CO 81301
-: 10:! USA
-: 11:!
-: 12:! Phone: 970-382-8979
-: 13:! e-mail: quetzal@quetzalcoatl.com
-: 14:!
-: 15:implicit none
-: 16:
-: 17:public :: next, value, convert_lower_case, check_eof, check_number
-: 18:
-: 19:integer, parameter, public :: nin = 10
-: 20:character (len=80), public :: card, field, lfield
-: 21:integer, public :: icpnt, free_format_error_flag
-: 22:logical, public :: eoff
-: 23:integer, parameter, private :: LONGreal = selected_real_kind(15,90)
-: 24:real (kind = LONGreal), public :: real_variable
-: 25:
-: 26:contains
-: 27:
38: 28: subroutine next ()
-: 29:!-----------------------------------------------
-: 30:! L o c a l P a r a m e t e r s
-: 31:!-----------------------------------------------
-: 32: integer, parameter :: icend = 80
-: 33: integer, parameter :: nodel = 3
-: 34: character (len=1), dimension(3), parameter :: delim = (/" ", ",", "="/)
-: 35:!-----------------------------------------------
-: 36:! L o c a l V a r i a b l e s
-: 37:!-----------------------------------------------
-: 38: integer :: i, j, k, iend, istart, kmin, err
-: 39: logical, save :: first = .true.
-: 40:!
-: 41:! on first pass, initialize icpnt pointer
-: 42:!
38: 43: if (first) then
1: 44: icpnt = icend + 1
1: 45: first = .false.
-: 46: end if
-: 47:!
-: 48:! if icpnt>icend, read the next record off unit 2 into
-: 49:! the string 'card'. next verify that this is a non-
-: 50:! blank card. if it is blank or an input comment card
-: 51:! (asterisk in column 1), skip it and get the next record
-: 52:!
-: 53:loop: do
61: 54: eoff = .false.
61: 55: if (icpnt > icend) then
24: 56: read (unit= nin, fmt="(a)", iostat = err) card
24: 57: if (err /= 0) then
1: 58: eoff = .true.
1: 59: icpnt = icend + 1
1: 60: return
-: 61:
-: 62: end if
23: 63: if (card(1:icend)==" " .or. card(1:1)=="*") then
-: 64: cycle loop
-: 65: end if
19: 66: icpnt = 1
-: 67: end if
-: 68:!
-: 69:! get the next sub-string. we do this as follows. we
-: 70:! look for the next delimeter. if it is as the same
-: 71:! position as the current pointer, we advance the pointer
-: 72:! and try again. if not, then the pointer is at the beginning
-: 73:! of a sub-string and the delimeter is trailing this sub-string.
-: 74:! note that we look for all the delimeters possible before
-: 75:! taking any action.
-: 76:!
931: 77: do i = icpnt, icend
-: 78: istart = i
-: 79: kmin = 0
3.6k: 80: do j = 1, nodel
2.7k: 81: k = index(card(i:icend),delim(j))
-: 82:!
-: 83:! index returns positions relative the beginning of the
-: 84:! sub-string. hence we add in the appropiate off-set to
-: 85:! give the index relative to the beginning of the string
-: 86:! card, not just the sub-string card(i:icend).
-: 87:!
3.6k: 88: if (k /= 0) then
948: 89: k = k + i - 1
948: 90: if (kmin == 0) then
-: 91: kmin = k
-: 92: else
36: 93: kmin = min(k,kmin)
-: 94: end if
-: 95: end if
-: 96: end do
-: 97:!
-: 98:! if kmin is not equal to the current pointer position, then
-: 99:! it must be pointing at the trailing delimeter of a valid
-: 100:! sub-string.
-: 101:!
931: 102: if (kmin /= i) then
37: 103: if (kmin > 0) then
37: 104: iend = kmin - 1
37: 105: exit loop
-: 106: end if
-: 107:!
-: 108:! if we fall through, there was no delimeter
-: 109:! found on the remainder of this record. this means
-: 110:! the entire remainder of this record is a valid sub-string
-: 111:!
-: 112: iend = icend
-: 113: exit loop
-: 114: end if
-: 115: end do
-: 116:!
-: 117:! if we fall through this loop, there were no more non-
-: 118:! delimeters on this record. go get next record
-: 119:!
23: 120: icpnt = icend + 1
-: 121: end do loop
-: 122:!
-: 123:! put the sub-string into the string 'field'. note
-: 124:! that fortran 77 pads the string with blanks
-: 125:!
37*: 126: field = card(istart:iend)
37: 127: icpnt = iend + 2
-: 128:!
-: 129: end subroutine next
-: 130:
18: 131: subroutine value (number, itype)
-: 132:!-----------------------------------------------
-: 133:! D u m m y A r g u m e n t s
-: 134:!-----------------------------------------------
-: 135: real (kind = LONGreal), intent(out) :: number
-: 136: integer, intent(out) :: itype
-: 137:!-----------------------------------------------
-: 138:! L o c a l V a r i a b l e s
-: 139:!-----------------------------------------------
-: 140: integer :: err
-: 141:!-----------------------------------------------
-: 142:!
-: 143:! get next field off unit 2
-: 144:!
18: 145: call next ()
18: 146: if (eoff) then
#####: 147: itype = -1
-: 148: else
-: 149:!
-: 150:! read field
-: 151:!
18: 152: read (unit = field, fmt="(bn,f20.0)", iostat = err) number
18: 153: if (err /= 0) then
#####: 154: itype = 0
-: 155: else
18: 156: itype = 1
-: 157: end if
-: 158: end if
-: 159:!
18: 160: end subroutine value
-: 161:
#####: 162: function convert_lower_case (input_string) result (output_string)
-: 163:!
-: 164: character (len=*), intent(in) :: input_string
-: 165: character (len=1) :: output_string
-: 166: integer :: collating_difference
-: 167:!
1.5k*: 168: if (ichar(input_string) >= ichar("A") .and. ichar(input_string) <= ichar("Z")) then
-: 169: collating_difference = ichar(input_string) - ichar("A")
1*: 170: output_string = char(ichar("a") + collating_difference)
-: 171: else
#####: 172: output_string = input_string
-: 173: end if
-: 174:!
#####: 175: end function convert_lower_case
-: 176:
18: 177: subroutine check_eof ()
-: 178:!
18: 179: if (free_format_error_flag == (-1)) then
#####: 180: print *," "
#####: 181: print *,"Abort. Unexpected end of file while reading input. "
#####: 182: print *,"Was reading the line:"
#####: 183: print *,card
#####: 184: print *," "
#####: 185: stop
-: 186: end if
-: 187:!
18: 188: end subroutine check_eof
-: 189:!
18: 190: subroutine check_number ()
-: 191:!
18: 192: if (free_format_error_flag == 0) then
#####: 193: print *," "
#####: 194: print *,"Abort. Expected a number on input and instead encountered the word: "
#####: 195: print *,field
#####: 196: print *," "
#####: 197: print *,"was reading the line:"
#####: 198: print *,card
#####: 199: print *," "
#####: 200: stop
-: 201: end if
-: 202:!
18: 203: end subroutine check_number
-: 204:
-: 205:
-: 206:
-: 207:end module free_input
-: 208:module read_input_m
-: 209:!
-: 210:! Copyright (C) 1997 by Quetzal Computational Associates, Inc.
-: 211:! Address any questions about this copyright notice to:
-: 212:!
-: 213:! Dr. John K. Prentice
-: 214:! Quetzal Computational Associates, Incorporated
-: 215:! 3455 Main Avenue, Suite 4
-: 216:! Durango, CO 81301
-: 217:! USA
-: 218:!
-: 219:! Phone: 970-382-8979
-: 220:! e-mail: quetzal@quetzalcoatl.com
-: 221:
-: 222:use free_input
-: 223:
-: 224:implicit none
-: 225:
-: 226:private
-: 227:public :: read_input
-: 228:
-: 229:integer, parameter, private :: LONGreal = selected_real_kind(15,90)
-: 230:
-: 231:
-: 232:contains
-: 233:
1: 234: subroutine read_input (input_file, dt, cmax, lambda, mu, yield_stress, R_infinity, b, &
-: 235: X_infinity, gamma, eta, plastic_strain_threshold, &
-: 236: failure_threshold, spin_frequency, wire_test, coil_test, &
-: 237: wire_radius, coil_radius, coil_pitch, radius_of_curvature, &
-: 238: time_steps_per_spin_cycle, number_of_sample_points, &
-: 239: dump_history, crack_closure_parameter, spin_cycle_dump_interval)
-: 240:!
-: 241:!=========== formal variables =============
-: 242:!
-: 243: character (len=*), intent(in) :: input_file
-: 244: real (kind = LONGreal), intent(out) :: dt, lambda, mu, yield_stress, R_infinity, b, &
-: 245: X_infinity, gamma, eta, &
-: 246: plastic_strain_threshold, wire_radius, &
-: 247: coil_radius, coil_pitch, radius_of_curvature, &
-: 248: failure_threshold, crack_closure_parameter
-: 249: integer, intent(out) :: spin_frequency, number_of_sample_points, cmax, &
-: 250: time_steps_per_spin_cycle, spin_cycle_dump_interval
-: 251: logical, intent(out) :: wire_test, coil_test, dump_history
-: 252:!
-: 253:!========== internal variables ============
-: 254:!
-: 255: real (kind = LONGreal) :: Poissons_ratio, Youngs_modulus, shear_modulus
-: 256: integer :: n
-: 257:!
1: 258: open (unit=nin,file=input_file,status="old",form="formatted",action="read")
1: 259: icpnt = 9999
-: 260:!
-: 261:! set defaults
-: 262:!
1: 263: cmax = huge(0)
1: 264: spin_cycle_dump_interval = huge(0)
1: 265: time_steps_per_spin_cycle = -999
1: 266: lambda = -1.0e20_LONGreal
1: 267: mu = -1.0e20_LONGreal
1: 268: yield_stress = -1.0e20_LONGreal
1: 269: R_infinity = -1.0e20_LONGreal
1: 270: b = -1.0e20_LONGreal
1: 271: X_infinity = -1.0e20_LONGreal
1: 272: gamma = -1.0e20_LONGreal
1: 273: eta = -1.0e20_LONGreal
1: 274: plastic_strain_threshold = -1.0e20_LONGreal
1: 275: spin_frequency = -999
-: 276: Poissons_ratio = -1.0e20_LONGreal
-: 277: Youngs_modulus = -1.0e20_LONGreal
-: 278: shear_modulus = -1.0e20_LONGreal
1: 279: wire_test = .false.
1: 280: coil_test = .false.
1: 281: wire_radius = -1.0e20_LONGreal
1: 282: coil_radius = -1.0e20_LONGreal
1: 283: coil_pitch = -1.0e20_LONGreal
1: 284: radius_of_curvature = -1.0e20_LONGreal
1: 285: number_of_sample_points = -999
1: 286: failure_threshold = -1.0e20_LONGreal
1: 287: dump_history = .false.
1: 288: crack_closure_parameter = -1.0e20_LONGreal
-: 289:!
-: 290:! parse input and read simulation parameters
-: 291:!
-: 292: do
20: 293: call next ()
20: 294: if (eoff) then
-: 295: exit
-: 296: else
-: 297: if (len(field) > 0) then
1.5k: 298: do n = 1, len(field)
1.5k: 299: field(n:n) = convert_lower_case(field(n:n))
-: 300: end do
-: 301: end if
-: 302:!
19: 303: if (field == "maximum_simulation_spin_cycles") then
1: 304: call value (real_variable, free_format_error_flag)
1: 305: call check_eof ()
1: 306: call check_number ()
1: 307: cmax = int(real_variable + 0.5_LONGreal)
18: 308: else if (field == "spin_cycle_dump_interval") then
1: 309: call value (real_variable, free_format_error_flag)
1: 310: call check_eof ()
1: 311: call check_number ()
1: 312: spin_cycle_dump_interval = int(real_variable + 0.5_LONGreal)
17: 313: else if (field == "time_steps_per_spin_cycle") then
1: 314: call value (real_variable, free_format_error_flag)
1: 315: call check_eof ()
1: 316: call check_number ()
1: 317: time_steps_per_spin_cycle = int(real_variable + 0.5_LONGreal)
16: 318: else if (field == "number_of_radial_sample_points") then
1: 319: call value (real_variable, free_format_error_flag)
1: 320: call check_eof ()
1: 321: call check_number ()
1: 322: number_of_sample_points = int(real_variable + 0.5_LONGreal)
15: 323: else if (field == "crack_closure_parameter") then
1: 324: call value (real_variable, free_format_error_flag)
1: 325: call check_eof ()
1: 326: call check_number ()
1: 327: crack_closure_parameter = real_variable
14: 328: else if (field == "poissons_ratio") then
#####: 329: call value (real_variable, free_format_error_flag)
#####: 330: call check_eof ()
#####: 331: call check_number ()
#####: 332: Poissons_ratio = real_variable
14: 333: else if (field == "youngs_modulus") then
1: 334: call value (real_variable, free_format_error_flag)
1: 335: call check_eof ()
1: 336: call check_number ()
1: 337: Youngs_modulus = real_variable
13: 338: else if (field == "shear_modulus") then
1: 339: call value (real_variable, free_format_error_flag)
1: 340: call check_eof ()
1: 341: call check_number ()
1: 342: shear_modulus = real_variable
12: 343: else if (field == "yield_stress") then
1: 344: call value (real_variable, free_format_error_flag)
1: 345: call check_eof ()
1: 346: call check_number ()
1: 347: yield_stress = real_variable
11: 348: else if (field == "asymptotic_stress_due_to_isotropic_strain_hardening") then
1: 349: call value (real_variable, free_format_error_flag)
1: 350: call check_eof ()
1: 351: call check_number ()
1: 352: R_infinity = real_variable
10: 353: else if (field == "rate_constant_for_isotropic_strain_hardening") then
1: 354: call value (real_variable, free_format_error_flag)
1: 355: call check_eof ()
1: 356: call check_number ()
1: 357: b = real_variable
9: 358: else if (field == "asymptotic_back_stress") then
1: 359: call value (real_variable, free_format_error_flag)
1: 360: call check_eof ()
1: 361: call check_number ()
1: 362: X_infinity = real_variable
8: 363: else if (field == "rate_constant_for_back_stress") then
1: 364: call value (real_variable, free_format_error_flag)
1: 365: call check_eof ()
1: 366: call check_number ()
1: 367: gamma = real_variable
7: 368: else if (field == "krajcinovic/lemaitre_damage_model_material_constant") then
1: 369: call value (real_variable, free_format_error_flag)
1: 370: call check_eof ()
1: 371: call check_number ()
1: 372: eta = real_variable
6: 373: else if (field == "accumulated_plastic_strain_threshold_for_damage") then
1: 374: call value (real_variable, free_format_error_flag)
1: 375: call check_eof ()
1: 376: call check_number ()
1: 377: plastic_strain_threshold = real_variable
5: 378: else if (field == "damage_threshold_for_material_failure") then
1: 379: call value (real_variable, free_format_error_flag)
1: 380: call check_eof ()
1: 381: call check_number ()
1: 382: failure_threshold = real_variable
4: 383: else if (field == "spin_test_frequency") then
1: 384: call value (real_variable, free_format_error_flag)
1: 385: call check_eof ()
1: 386: call check_number ()
1: 387: spin_frequency = int(real_variable + 0.5_LONGreal)
3: 388: else if (field == "simulate_wire_spin_test") then
1: 389: wire_test = .true.
2: 390: else if (field == "simulate_coil_spin_test") then
#####: 391: coil_test = .true.
2: 392: else if (field == "dump_constitutive_history") then
#####: 393: dump_history = .true.
2: 394: else if (field == "wire_radius" .or. field == "filar_radius") then
1: 395: call value (real_variable, free_format_error_flag)
1: 396: call check_eof ()
1: 397: call check_number ()
1: 398: wire_radius = real_variable
1: 399: else if (field == "coil_radius") then
#####: 400: call value (real_variable, free_format_error_flag)
#####: 401: call check_eof ()
#####: 402: call check_number ()
#####: 403: coil_radius = real_variable
1: 404: else if (field == "coil_pitch") then
#####: 405: call value (real_variable, free_format_error_flag)
#####: 406: call check_eof ()
#####: 407: call check_number ()
#####: 408: coil_pitch = real_variable
1: 409: else if (field == "radius_of_curvature") then
1: 410: call value (real_variable, free_format_error_flag)
1: 411: call check_eof ()
1: 412: call check_number ()
1: 413: radius_of_curvature = real_variable
-: 414: else
#####: 415: print *," "
#####: 416: print *,"unrecognized word in simulation input, abort."
#####: 417: print *," "
#####: 418: print *,"The unrecognized word was: ",field
#####: 419: print *," "
#####: 420: stop
-: 421: end if
-: 422: end if
-: 423: end do
-: 424:!
1: 425: close (unit = nin)
-: 426:!
-: 427:! compute time step
-: 428:!
1: 429: if (spin_frequency > 0 .and. time_steps_per_spin_cycle > 0) then
-: 430: dt = 1.0_LONGreal/(real(spin_frequency,LONGreal) * &
1: 431: real(time_steps_per_spin_cycle,LONGreal))
-: 432: end if
-: 433:!
-: 434:! check input to catch obvious errors
-: 435:!
1*: 436: if (.not. wire_test .and. .not. coil_test) then
#####: 437: print *," "
#####: 438: print *,"The type of test was not specified in the input."
#####: 439: print *,"You must either specify simulate_wire_spin_test or simulate_coil_spin_test."
#####: 440: print *,"Simulation aborted. Please correct the input and rerun."
#####: 441: print *," "
#####: 442: stop
1: 443: else if (wire_test .and. coil_test) then
#####: 444: print *," "
#####: 445: print *,"You have specified both simulate_wire_spin_test and simulate_coil_spin_test."
#####: 446: print *,"Only one of these can be specified for a given simulation."
#####: 447: print *,"Simulation aborted. Please correct the input and rerun."
#####: 448: print *," "
#####: 449: stop
-: 450: end if
-: 451:!
-: 452: n = 0
1: 453: if (Poissons_ratio > 0.0_LONGreal) then
-: 454: n = n + 1
-: 455: end if
1: 456: if (Youngs_modulus > 0.0_LONGreal) then
1: 457: n = n + 1
-: 458: end if
1: 459: if (shear_modulus > 0.0_LONGreal) then
1: 460: n = n + 1
-: 461: end if
-: 462:!
1: 463: if (n < 2) then
#####: 464: print *," "
#####: 465: print *,"You must specify two of the following items in the input."
#####: 466: print *," Poissons_ratio"
#####: 467: print *," Youngs_modulus"
#####: 468: print *," shear_modulus"
#####: 469: if (n == 0) then
#####: 470: print *,"None of these was specified."
-: 471: else if (n == 1) then
#####: 472: if (Poissons_ratio > 0.0_LONGreal) then
#####: 473: print *,"Only the Poissons_ratio was specified."
#####: 474: else if (Youngs_modulus > 0.0_LONGreal) then
#####: 475: print *,"Only the Youngs_modulus was specified."
#####: 476: else if (shear_modulus > 0.0_LONGreal) then
#####: 477: print *,"Only the shear_modulus was specified."
-: 478: end if
-: 479: end if
#####: 480: print *,"Simulation aborted. Please correct the input and rerun."
#####: 481: print *," "
#####: 482: stop
1: 483: else if (crack_closure_parameter < -1.0e19_LONGreal) then
#####: 484: print *," "
#####: 485: print *,"crack_closure_parameter was not set in the input."
#####: 486: print *,"Simulation aborted. Please correct the input and rerun."
#####: 487: print *," "
#####: 488: stop
1: 489: else if (crack_closure_parameter <= 0.0_LONGreal) then
#####: 490: print *," "
#####: 491: print *,"The crack_closure_parameter was less than or equal to zero."
#####: 492: print *,"This value must be greater than zero and less than or equal to one."
#####: 493: print *,"Simulation aborted. Please correct the input and rerun."
#####: 494: print *," "
#####: 495: stop
1: 496: else if (crack_closure_parameter > 1.0_LONGreal) then
#####: 497: print *," "
#####: 498: print *,"The crack_closure_parameter was greater than one."
#####: 499: print *,"This value must be greater than zero and less than or equal to one."
#####: 500: print *,"Simulation aborted. Please correct the input and rerun."
#####: 501: print *," "
#####: 502: stop
1: 503: else if (abs(Poissons_ratio - 0.5_LONGreal) < 1.0e-10_LONGreal) then
#####: 504: print *," "
#####: 505: print *,"The Poissons_ratio was set to one-half. This value is physically"
#####: 506: print *,"unrealistic (it is appropriate only to fluids)."
#####: 507: print *,"Simulation aborted. Please correct the input and rerun."
#####: 508: print *," "
#####: 509: stop
1: 510: else if (time_steps_per_spin_cycle == -999) then
#####: 511: print *," "
#####: 512: print *,"time_steps_per_spin_cycle was not set in the input."
#####: 513: print *,"Simulation aborted. Please correct the input and rerun."
#####: 514: print *," "
#####: 515: stop
1: 516: else if (time_steps_per_spin_cycle == 0) then
#####: 517: print *," "
#####: 518: print *,"The time_steps_per_spin_cycle in the input was set to zero."
#####: 519: print *,"This value must be an integer greater than zero."
#####: 520: print *,"Simulation aborted. Please correct the input and rerun."
#####: 521: print *," "
#####: 522: stop
1: 523: else if (number_of_sample_points == -999) then
#####: 524: print *," "
#####: 525: print *,"number_of_radial_sample_points was not set in the input."
#####: 526: print *,"Simulation aborted. Please correct the input and rerun."
#####: 527: print *," "
#####: 528: stop
1: 529: else if (number_of_sample_points == 0) then
#####: 530: print *," "
#####: 531: print *,"number_of_radial_sample_points in the input was set to zero."
#####: 532: print *,"This value must be an integer greater than zero."
#####: 533: print *,"Simulation aborted. Please correct the input and rerun."
#####: 534: print *," "
#####: 535: stop
1: 536: else if (spin_frequency == -999) then
#####: 537: print *," "
#####: 538: print *,"spin_test_frequency was not set in the input."
#####: 539: print *,"Simulation aborted. Please correct the input and rerun."
#####: 540: print *," "
#####: 541: stop
1: 542: else if (spin_frequency == 0) then
#####: 543: print *," "
#####: 544: print *,"The spin_test_frequency in the input was set to zero."
#####: 545: print *,"This value must be an integer greater than zero."
#####: 546: print *,"Simulation aborted. Please correct the input and rerun."
#####: 547: print *," "
#####: 548: stop
1: 549: else if (cmax <= 0) then
#####: 550: print *," "
#####: 551: print *,"The maximum_simulation_spin_cycles are negative or zero in the input."
#####: 552: print *,"Simulation aborted. Please correct the input and rerun."
#####: 553: print *," "
#####: 554: stop
1: 555: else if (spin_cycle_dump_interval <= 0) then
#####: 556: print *," "
#####: 557: print *,"The spin_cycle_dump_interval is negative or zero in the input."
#####: 558: print *,"Simulation aborted. Please correct the input and rerun."
#####: 559: print *," "
#####: 560: stop
1: 561: else if (yield_stress < -1.0e19_LONGreal) then
#####: 562: print *," "
#####: 563: print *,"yield_stress was not set in the input."
#####: 564: print *,"Simulation aborted. Please correct the input and rerun."
#####: 565: print *," "
#####: 566: stop
1: 567: else if (R_infinity < -1.0e19_LONGreal) then
#####: 568: print *," "
#####: 569: print *,"aymptotic_stress_due_to_isotropic_strain_hardening was not set in the input."
#####: 570: print *,"Simulation aborted. Please correct the input and rerun."
#####: 571: print *," "
#####: 572: stop
1: 573: else if (b < -1.0e19_LONGreal) then
#####: 574: print *," "
#####: 575: print *,"rate_constant_for_isotropic_strain_hardening was not set in the input."
#####: 576: print *,"Simulation aborted. Please correct the input and rerun."
#####: 577: print *," "
#####: 578: stop
1: 579: else if (X_infinity < -1.0e19_LONGreal) then
#####: 580: print *," "
#####: 581: print *,"asymptotic_back_stress was not set in the input."
#####: 582: print *,"Simulation aborted. Please correct the input and rerun."
#####: 583: print *," "
#####: 584: stop
1: 585: else if (gamma < -1.0e19_LONGreal) then
#####: 586: print *," "
#####: 587: print *,"rate_constant_for_back_stress was not set in the input."
#####: 588: print *,"Simulation aborted. Please correct the input and rerun."
#####: 589: print *," "
#####: 590: stop
1: 591: else if (eta < -1.0e19_LONGreal) then
#####: 592: print *," "
#####: 593: print *,"krajcinovic/lemaitre_damage_model_material_constant was not set in the input."
#####: 594: print *,"Simulation aborted. Please correct the input and rerun."
#####: 595: print *," "
#####: 596: stop
1: 597: else if (plastic_strain_threshold < -1.0e19_LONGreal) then
#####: 598: print *," "
#####: 599: print *,"accumulated_plastic_strain_threshold_for_damage was not set in the input."
#####: 600: print *,"Simulation aborted. Please correct the input and rerun."
#####: 601: print *," "
#####: 602: stop
1: 603: else if (wire_radius < -1.0e19_LONGreal) then
#####: 604: print *," "
#####: 605: if (wire_test) then
#####: 606: print *,"wire_radius (or equivalently, filar_radius) was not set in the input."
#####: 607: else if (coil_test) then
#####: 608: print *,"filar_radius (or equivalently, wire_radius) was not set in the input."
-: 609: end if
#####: 610: print *,"Simulation aborted. Please correct the input and rerun."
#####: 611: print *," "
#####: 612: stop
1*: 613: else if (coil_test .and. coil_radius < -1.0e19_LONGreal) then
#####: 614: print *," "
#####: 615: print *,"coil_radius was not set in the input."
#####: 616: print *,"Simulation aborted. Please correct the input and rerun."
#####: 617: print *," "
#####: 618: stop
1*: 619: else if (coil_test .and. coil_pitch < -1.0e19_LONGreal) then
#####: 620: print *," "
#####: 621: print *,"coil_pitch was not set in the input."
#####: 622: print *,"Simulation aborted. Please correct the input and rerun."
#####: 623: print *," "
#####: 624: stop
1: 625: else if (radius_of_curvature < -1.0e19_LONGreal) then
#####: 626: print *," "
#####: 627: print *,"radius_of_curvature was not set in the input."
#####: 628: print *,"Simulation aborted. Please correct the input and rerun."
#####: 629: print *," "
#####: 630: stop
1: 631: else if (failure_threshold < -1.0e19_LONGreal) then
#####: 632: print *," "
#####: 633: print *,"damage_threshold_for_material_failure was not set in the input."
#####: 634: print *,"Simulation aborted. Please correct the input and rerun."
#####: 635: print *," "
#####: 636: stop
1: 637: else if (failure_threshold <= 0.0_LONGreal .or. failure_threshold > 1.0_LONGreal) then
#####: 638: print *," "
#####: 639: print *,"damage_threshold_for_material_failure was invalid in the input."
#####: 640: print *,"The value must be greater than zero and less than or equal to one."
#####: 641: print *,"In the input, the value was ", failure_threshold
#####: 642: print *,"Simulation aborted. Please correct the input and rerun."
#####: 643: print *," "
#####: 644: stop
-: 645: end if
-: 646:!
-: 647:! compute Lame constants
-: 648:!
1*: 649: if (Poissons_ratio > 0.0_LONGreal .and. Youngs_modulus > 0.0_LONGreal) then
-: 650: lambda = -(Youngs_modulus * Poissons_ratio)/(2.0_LONGreal * Poissons_ratio**2 + &
#####: 651: Poissons_ratio - 1.0_LONGreal)
#####: 652: mu = 0.5_LONGreal * Youngs_modulus / (1.0_LONGreal + Poissons_ratio)
1*: 653: else if (Poissons_ratio > 0.0_LONGreal .and. shear_modulus > 0.0_LONGreal) then
#####: 654: mu = shear_modulus
-: 655: lambda = (2.0_LONGreal * mu * Poissons_ratio)/(1.0_LONGreal - 2.0_LONGreal * &
#####: 656: Poissons_ratio)
1: 657: else if (Youngs_modulus > 0.0_LONGreal .and. shear_modulus > 0.0_LONGreal) then
1: 658: mu = shear_modulus
-: 659: lambda = mu * (2.0_LONGreal * mu - Youngs_modulus)/(3.0_LONGreal * mu - &
1: 660: Youngs_modulus)
-: 661: end if
-: 662:!
1: 663: close (unit=nin)
-: 664:!
1: 665: end subroutine read_input
-: 666:
-: 667:end module read_input_m
-: 668:module computer_time_m
-: 669:
-: 670:!
-: 671:! copyright (c) 1997 by quetzal computational associates, inc.
-: 672:! address any questions about this copyright notice to:
-: 673:!
-: 674:! john k. prentice
-: 675:! quetzal computational associates, incorporated
-: 676:! 3455 main avenue, suite 4
-: 677:! durango, co 81301
-: 678:! usa
-: 679:!
-: 680:! phone: 970-382-8979
-: 681:! e-mail: quetzal@quetzalcoatl.com
-: 682:!
-: 683:
-: 684:implicit none
-: 685:
-: 686:public :: computer_time
-: 687:
-: 688:!
-: 689:integer, parameter, private :: LONGreal = selected_real_kind(15,90)
-: 690:
-: 691:contains
-: 692:
-: 693:
2: 694: subroutine computer_time (tnow)
-: 695:!
-: 696:! purpose: return the elapsed system time for this
-: 697:! process on a sun microsystems computer.
-: 698:! coded: f90 version coded 1 may 1993
-: 699:! author: john k. prentice
-: 700:!
-: 701:! input: none
-: 702:!
-: 703:! output:
-: 704:!
-: 705:! tnow real elapsed system time (seconds)
-: 706:!
-: 707:!
-: 708: real (kind = LONGreal), intent(out) :: tnow
-: 709:!
-: 710: logical, save :: first = .true.
-: 711: logical, save :: first_flip
-: 712: integer :: counted, count_rate, count_max
-: 713: real (kind = LONGreal) :: trate, tmax
-: 714: real (kind = LONGreal), save :: tfirst
-: 715:!
2: 716: call system_clock (counted, count_rate, count_max)
2: 717: if (counted < 0 .or. count_rate == 0) then
#####: 718: tnow = 0.0_LONGreal
-: 719: else
2: 720: tnow = real(counted,LONGreal)
2: 721: trate = real(count_rate,LONGreal)
2: 722: tnow = tnow/trate
-: 723:!
2: 724: if (first) then
1: 725: first = .false.
1: 726: tfirst = tnow
1: 727: first_flip = .true.
1: 728: else if (tnow < tfirst) then
#####: 729: if (.not. first_flip) then
#####: 730: tmax = real(count_max,LONGreal)/trate
#####: 731: tfirst = tfirst - tmax
-: 732: else
#####: 733: tmax = real(count_max,LONGreal)/trate
#####: 734: tfirst = -(tmax - tfirst)
#####: 735: first_flip = .false.
-: 736: end if
-: 737: end if
-: 738:!
2: 739: tnow = tnow - tfirst
-: 740: end if
-: 741:!
2: 742: end subroutine computer_time
-: 743:
-: 744:end module computer_time_m
-: 745:
-: 746:module perdida_m
-: 747:!
-: 748:! Copyright (C) 1997 by Quetzal Computational Associates, Inc.
-: 749:! Address any questions about this copyright notice to:
-: 750:!
-: 751:! Dr. John K. Prentice
-: 752:! Quetzal Computational Associates, Incorporated
-: 753:! 3455 Main Avenue, Suite 4
-: 754:! Durango, CO 81301
-: 755:! USA
-: 756:!
-: 757:! Phone: 970-382-8979
-: 758:! e-mail: quetzal@quetzalcoatl.com
-: 759:!
-: 760:
-: 761:implicit none
-: 762:
-: 763:public :: perdida
-: 764:private :: generalized_hookes_law, damage_rate
-: 765:
-: 766:integer, parameter, private :: LONGreal = selected_real_kind(15,90)
-: 767:
-: 768:contains
-: 769:
585.9M: 770: subroutine perdida (dt, lambda, mu, yield_stress, R_infinity, b, X_infinity, gamma, &
585.9M: 771: eta, plastic_strain_threshold, stress_tensor, strain_tensor, &
585.9M: 772: plastic_strain_tensor, strain_rate_tensor, &
585.9M: 773: accumulated_plastic_strain, back_stress_tensor, &
-: 774: isotropic_hardening_stress, damage, failure_threshold, &
-: 775: crack_closure_parameter)
-: 776:!
-: 777:! Author: Dr. John K. Prentice
-: 778:! Affiliation: Quetzal Computational Associates, Inc.
-: 779:! Dates: 28 November 1997
-: 780:!
-: 781:! Purpose: Cumulative damage constitutive model for ductile materials.
-: 782:!
-: 783:!############################################################################################
-: 784:!
-: 785:! Input:
-: 786:!
-: 787:! dt [selected_real_kind(15,90)]
-: 788:! time step (sec)
-: 789:!
-: 790:! lambda [selected_real_kind(15,90)]
-: 791:! Lame constant Lambda (dynes/cm**2)
-: 792:!
-: 793:! mu [selected_real_kind(15,90)]
-: 794:! Lame constant mu (dynes/cm**2)
-: 795:!
-: 796:! yield_stress [selected_real_kind(15,90)]
-: 797:! yield stress of material (dynes/cm**2)
-: 798:!
-: 799:! R_infinity [selected_real_kind(15,90)]
-: 800:! asymptotic stress due to isotropic strain hardening
-: 801:! (dynes/cm**2)
-: 802:!
-: 803:! b [selected_real_kind(15,90)]
-: 804:! rate constant for stress due to isotropic strain
-: 805:! hardening (dimensionless)
-: 806:!
-: 807:! X_infinity [selected_real_kind(15,90)]
-: 808:! asymptotic back stress (dynes/cm**2)
-: 809:!
-: 810:! gamma [selected_real_kind(15,90)]
-: 811:! rate constant for back stress
-: 812:! (dimensionless)
-: 813:!
-: 814:! eta [selected_real_kind(15,90)]
-: 815:! material constant for the Krajcinovic/Lemaitre ductile
-: 816:! damage model (ergs/gm)
-: 817:!
-: 818:! plastic_strain_threshold [selected_real_kind(15,90)]
-: 819:! accumulated plastic strain threshold for damage to
-: 820:! occur (dimensionless)
-: 821:!
-: 822:!
-: 823:! strain_tensor [selected_real_kind(15,90), dimension(3,3)]
-: 824:! total strain tensor at end of time step
-: 825:! (dimensionless)
-: 826:!
-: 827:! plastic_strain_tensor [selected_real_kind(15,90), dimension(3,3)]
-: 828:! total plastic strain tensor at beginning of time step
-: 829:! (dimensionless)
-: 830:!
-: 831:! strain_rate_tensor [selected_real_kind(15,90), dimension(3,3)]
-: 832:! total strain rate tensor at beginning of time step
-: 833:! (1/sec)
-: 834:!
-: 835:! accumulated_plastic_strain [selected_real_kind(15,90)]
-: 836:! accumulated plastic strain at beginning of time step
-: 837:! (dimensionless)
-: 838:!
-: 839:! back_stress_tensor [selected_real_kind(15,90)]
-: 840:! back stress tensor at the beginning of time step
-: 841:! (dynes/cm**2)
-: 842:!
-: 843:! isotropic_hardening_stress [selected_real_kind(15,90)]
-: 844:! stress due to isotropic strain hardening at beginning
-: 845:! of time step (dynes/cm**2)
-: 846:!
-: 847:! damage [selected_real_kind(15,90)]
-: 848:! cumulative damage at beginning of time step
-: 849:! (dimensionless)
-: 850:!
-: 851:! failure_threshold [selected_real_kind(15,90)]
-: 852:! damage at which material failure occurs.
-: 853:! (dimensionless)
-: 854:!
-: 855:! crack_closure_parameter [selected_real_kind(15,90)]
-: 856:! if the material is in compression, this parameter
-: 857:! accounts for partial closure of cracks when
-: 858:! calculating the stress
-: 859:!
-: 860:!
-: 861:! Output:
-: 862:!
-: 863:! stress_tensor [selected_real_kind(15,90), dimension(3,3)]
-: 864:! total stress tensor at end of time step (dynes/cm**2)
-: 865:!
-: 866:! plastic_strain_tensor [selected_real_kind(15,90), dimension(3,3)]
-: 867:! total plastic strain tensor at end of time step
-: 868:! (dimensionless)
-: 869:!
-: 870:! accumulated_plastic_strain [selected_real_kind(15,90)]
-: 871:! accumulated plastic strain at end of time step
-: 872:! (dimensionless)
-: 873:!
-: 874:! back_stress_tensor [selected_real_kind(15,90)]
-: 875:! back stress tensor at the end of time step
-: 876:! (dynes/cm**2)
-: 877:!
-: 878:! isotropic_hardening_stress [selected_real_kind(15,90)]
-: 879:! stress due to isotropic strain hardening at end of
-: 880:! time step (dynes/cm**2)
-: 881:!
-: 882:! damage [selected_real_kind(15,90)]
-: 883:! cumulative damage at end of time step (dimensionless)
-: 884:!
-: 885:!############################################################################################
-: 886:!
-: 887:!
-: 888:!=========== formal variables =============
-: 889:!
-: 890: real (kind = LONGreal), intent(in) :: dt, yield_stress, lambda, mu, R_infinity, b, &
-: 891: X_infinity, gamma, eta, failure_threshold, &
-: 892: plastic_strain_threshold, &
-: 893: crack_closure_parameter
-: 894: real (kind = LONGreal), dimension(:,:), intent(in) :: strain_rate_tensor, &
-: 895: strain_tensor
-: 896: real (kind = LONGreal), dimension(:,:), intent(inout) :: plastic_strain_tensor, &
-: 897: back_stress_tensor
-: 898: real (kind = LONGreal), dimension(:,:), intent(out) :: stress_tensor
-: 899: real (kind = LONGreal), intent(inout) :: damage, accumulated_plastic_strain, &
-: 900: isotropic_hardening_stress
-: 901:!
-: 902:!========== internal variables ============
-: 903:!
-: 904: integer :: i, j
-: 905: real (kind = LONGreal), dimension(3,3) :: deviatoric_stress_tensor, &
-: 906: deviatoric_strain_rate_tensor, &
-: 907: damaged_dev_stress_tensor, &
-: 908: plastic_strain_rate_tensor, &
-: 909: back_stress_rate_tensor, &
-: 910: deviatoric_strain_tensor, &
-: 911: old_plastic_strain_tensor
-: 912: real (kind = LONGreal) :: equivalent_stress, isotropic_hardening_stress_rate, &
-: 913: triaxiality_ratio, Poissons_ratio, yield_function, &
-: 914: Youngs_modulus, term1, term2, term3, xi_dot, &
-: 915: hydrostatic_stress, hydrostatic_strain_rate, &
-: 916: yield_function_dot, plastic_time_step, damage_dot, &
-: 917: discriminant, strain_energy_density, frac, &
-: 918: accumulated_plastic_strain_rate, hydrostatic_strain, &
-: 919: crack_parameter
-: 920:!
585.9M: 921: Poissons_ratio = 0.5_LONGreal * lambda / (lambda + mu)
-: 922: Youngs_modulus = (1.0_LONGreal - 2.0_LONGreal * Poissons_ratio) * (3.0_LONGreal * &
585.9M: 923: lambda + 2.0_LONGreal * mu)
-: 924:!
-: 925:! make an initial calculation of the stress
-: 926:!
-: 927: stress_tensor(:,:) = (1.0_LONGreal - damage) * &
-: 928: generalized_hookes_law (strain_tensor(:,:) - &
14.6G: 929: plastic_strain_tensor(:,:), lambda, mu)
-: 930:!
-: 931:! compute hydrostatic stress at beginning of time step
-: 932:!
-: 933: hydrostatic_stress = (stress_tensor(1,1) + stress_tensor(2,2) + &
585.9M: 934: stress_tensor(3,3))/3.0_LONGreal
-: 935:!
-: 936:! if in compression, modify the stress to account for partial closure of cracks
-: 937:!
585.9M: 938: if (hydrostatic_stress < -1.0e-10_LONGreal .and. damage > 0.0_LONGreal) then
99.7M: 939: crack_parameter = crack_closure_parameter
-: 940: stress_tensor(:,:) = (1.0_LONGreal - crack_parameter * damage) / &
1.3G: 941: (1.0_LONGreal - damage) * stress_tensor(:,:)
-: 942: hydrostatic_stress = (stress_tensor(1,1) + stress_tensor(2,2) + &
99.7M: 943: stress_tensor(3,3))/3.0_LONGreal
-: 944: else
-: 945: crack_parameter = 1.0_LONGreal
-: 946: end if
-: 947:!
-: 948:! compute deviatoric stress tensor at beginning of time step
-: 949:!
7.6G: 950: deviatoric_stress_tensor(:,:) = stress_tensor(:,:)
585.9M: 951: deviatoric_stress_tensor(1,1) = stress_tensor(1,1) - hydrostatic_stress
585.9M: 952: deviatoric_stress_tensor(2,2) = stress_tensor(2,2) - hydrostatic_stress
585.9M: 953: deviatoric_stress_tensor(3,3) = stress_tensor(3,3) - hydrostatic_stress
-: 954:!
-: 955: damaged_dev_stress_tensor(:,:) = deviatoric_stress_tensor(:,:)/(1.0_LONGreal - &
7.6G: 956: crack_parameter * damage)
-: 957:!
-: 958:! compute the yield function
-: 959:!
-: 960: equivalent_stress = sqrt(1.5_LONGreal* sum((damaged_dev_stress_tensor(:,:) &
7.6G: 961: - back_stress_tensor(:,:))**2))
585.9M: 962: yield_function = equivalent_stress - isotropic_hardening_stress - yield_stress
-: 963:!
-: 964:! if the yeild function is non-negative, compute the time derivative also. If
-: 965:! the yield function is negative, arbitrarily set the time derivative to a negative
-: 966:! number. It doesn't matter what the value is, we just want to set a flag to
-: 967:! identify that no plastic strain is being accumulated in this time step.
-: 968:!
-: 969:! For the purposes of computing the time derivative of the yield function, for
-: 970:! computational convenience, we ignore the time derivative of damage and also
-: 971:! of the back stress.
-: 972:!
585.9M: 973: if (yield_function < 0.0_LONGreal .or. equivalent_stress < 1.0e-10_LONGreal) then
-: 974: yield_function_dot = -1.0_LONGreal
-: 975: else
-: 976: hydrostatic_strain_rate = (strain_rate_tensor(1,1) + strain_rate_tensor(2,2) + &
-: 977: strain_rate_tensor(3,3))/3.0_LONGreal
-: 978: deviatoric_strain_rate_tensor(:,:) = strain_rate_tensor(:,:)
-: 979: deviatoric_strain_rate_tensor(1,1) = strain_rate_tensor(1,1) - &
-: 980: hydrostatic_strain_rate
-: 981: deviatoric_strain_rate_tensor(2,2) = strain_rate_tensor(2,2) - &
-: 982: hydrostatic_strain_rate
-: 983: deviatoric_strain_rate_tensor(3,3) = strain_rate_tensor(3,3) - &
-: 984: hydrostatic_strain_rate
-: 985:!
-: 986: yield_function_dot = 1.5_LONGreal / equivalent_stress * &
-: 987: sum((damaged_dev_stress_tensor(:,:) - back_stress_tensor(:,:)) * &
872.4M: 988: generalized_hookes_law (strain_rate_tensor(:,:), lambda, mu))
-: 989: end if
-: 990:!
-: 991:! if the yield function and/or its time derivative are negative, then we are done.
-: 992:! In this case, the strain is all elastic and the stress we have calculated is
-: 993:! correct. However, if the yield function and its time derivative are non-negative,
-: 994:! then we have exceeded the yield condition. In this case, plastic work is being
-: 995:! done and we need to yield the material.
-: 996:!
585.9M: 997: if (yield_function >= 0.0_LONGreal .and. yield_function_dot >= 0.0_LONGreal) then
766.6M: 998: old_plastic_strain_tensor(:,:) = plastic_strain_tensor(:,:)
-: 999:!
-: 1000:! modify the stress tensor to account for plasticity
-: 1001:!
766.6M: 1002: term1 = sum(damaged_dev_stress_tensor(:,:)**2)
-: 1003: term2 = -2.0_LONGreal * sum(damaged_dev_stress_tensor(:,:) * &
766.6M: 1004: back_stress_tensor(:,:))
-: 1005: term3 = sum(back_stress_tensor(:,:)**2) - (isotropic_hardening_stress + &
766.6M: 1006: yield_stress)**2 / 1.5_LONGreal
59.0M: 1007: discriminant = term2**2 - 4.0_LONGreal * term1 * term3
-: 1008:!
59.0M: 1009: if (discriminant < 0.0_LONGreal) then
#####: 1010: print *,"discriminant is negative in perdida, abort."
#####: 1011: stop
-: 1012: end if
-: 1013:!
59.0M: 1014: frac = (-term2 + sqrt(discriminant))/(2.0_LONGreal * term1)
59.0M: 1015: if (frac < 0.0_LONGreal .or. frac > 1.0_LONGreal) then
#####: 1016: print *,"bad plastic strain fraction in perdida, abort.",frac
#####: 1017: stop
-: 1018: end if
-: 1019:!
766.6M: 1020: deviatoric_stress_tensor(:,:) = frac * deviatoric_stress_tensor(:,:)
766.6M: 1021: stress_tensor(:,:) = deviatoric_stress_tensor(:,:)
59.0M: 1022: stress_tensor(1,1) = stress_tensor(1,1) + hydrostatic_stress
59.0M: 1023: stress_tensor(2,2) = stress_tensor(2,2) + hydrostatic_stress
59.0M: 1024: stress_tensor(3,3) = stress_tensor(3,3) + hydrostatic_stress
-: 1025:!
-: 1026:! somewhat inconsistently, compute the real accumulated plastic strain from the
-: 1027:! total strain tensor using the yield function to pick out the elastic piece of
-: 1028:! the strain
-: 1029:!
-: 1030: hydrostatic_strain = (strain_tensor(1,1) + strain_tensor(2,2) + &
59.0M: 1031: strain_tensor(3,3)) / 3.0_LONGreal
766.6M: 1032: deviatoric_strain_tensor(:,:) = strain_tensor(:,:)
59.0M: 1033: deviatoric_strain_tensor(1,1) = deviatoric_strain_tensor(1,1) - hydrostatic_strain
59.0M: 1034: deviatoric_strain_tensor(2,2) = deviatoric_strain_tensor(2,2) - hydrostatic_strain
59.0M: 1035: deviatoric_strain_tensor(3,3) = deviatoric_strain_tensor(3,3) - hydrostatic_strain
-: 1036: plastic_strain_tensor(:,:) = frac * plastic_strain_tensor(:,:) + (1.0_LONGreal - &
766.6M: 1037: frac) * deviatoric_strain_tensor(:,:)
-: 1038:!
-: 1039:! The time over which the plasticity and damage parameters are to be integrated
-: 1040:! may be shorter than the current time step since we may have been elastic during
-: 1041:! part of this time step. So compute the plastic time step.
-: 1042:!
59.0M: 1043: plastic_time_step = (1.0_LONGreal - frac) * dt
-: 1044:!
-: 1045:! compute the plastic strain rate tensor
-: 1046:!
59.0M: 1047: if (plastic_time_step <= 0.0_LONGreal) then
#####: 1048: plastic_strain_rate_tensor(:,:) = 0.0_LONGreal
-: 1049: else
-: 1050: plastic_strain_rate_tensor(:,:) = (plastic_strain_tensor(:,:) - &
766.6M: 1051: old_plastic_strain_tensor(:,:)) / plastic_time_step
-: 1052: end if
-: 1053:!
-: 1054:! Compute the accumulated plastic strain rate and the time derivative of
-: 1055:! the Lagrange multiplier.
-: 1056:!
-: 1057: accumulated_plastic_strain_rate = sqrt(sum(plastic_strain_rate_tensor(:,:)**2) / &
766.6M: 1058: 1.5_LONGreal)
-: 1059: xi_dot = (1.0_LONGreal - crack_parameter * damage) * &
59.0M: 1060: accumulated_plastic_strain_rate
-: 1061:!
-: 1062:! Next integrate the accumulated plastic strain rate to get the
-: 1063:! new accumulated plastic strain
-: 1064:!
-: 1065: accumulated_plastic_strain = accumulated_plastic_strain + plastic_time_step * &
59.0M: 1066: accumulated_plastic_strain_rate
-: 1067:!
-: 1068:! Calculate the new back stress tensor and the isotropic stress
-: 1069:!
-: 1070: back_stress_rate_tensor(:,:) = gamma * (X_infinity * (1.0_LONGreal - &
-: 1071: crack_parameter * damage) / 1.5_LONGreal * &
-: 1072: plastic_strain_rate_tensor(:,:) - &
766.6M: 1073: back_stress_tensor(:,:) * xi_dot)
-: 1074: isotropic_hardening_stress_rate = b * (R_infinity - isotropic_hardening_stress) * &
59.0M: 1075: xi_dot
-: 1076: back_stress_tensor(:,:) = back_stress_tensor(:,:) + &
766.6M: 1077: plastic_time_step * back_stress_rate_tensor(:,:)
-: 1078: isotropic_hardening_stress = isotropic_hardening_stress + &
59.0M: 1079: plastic_time_step * isotropic_hardening_stress_rate
-: 1080:!
-: 1081:! Finally, calculate the damage if the hydrostatic stress is positive (we assume
-: 1082:! negative is compression and positive is tension)
-: 1083:!
59.0M: 1084: if (hydrostatic_stress > -1.0e-10_LONGreal) then
383.3M: 1085: equivalent_stress = sqrt(1.5_LONGreal* sum(damaged_dev_stress_tensor(:,:)**2))
-: 1086: triaxiality_ratio = 2.0_LONGreal / 3.0_LONGreal * (1.0_LONGreal + &
-: 1087: Poissons_ratio) + 3.0_LONGreal * (1.0_LONGreal - &
-: 1088: 2.0_LONGreal * Poissons_ratio) * (hydrostatic_stress / &
29.5M: 1089: equivalent_stress)**2
-: 1090: strain_energy_density = 0.5_LONGreal * equivalent_stress**2 * &
29.5M: 1091: triaxiality_ratio / Youngs_modulus
-: 1092: damage_dot = damage_rate (eta, accumulated_plastic_strain_rate, &
-: 1093: strain_energy_density, accumulated_plastic_strain, &
-: 1094: plastic_strain_threshold)
29.5M: 1095: damage = min(1.0_LONGreal, damage + damage_dot * plastic_time_step)
-: 1096:!
-: 1097:! if the damage exceeds the failure threshold, fail the material completely
-: 1098:!
29.5M: 1099: if (damage >= failure_threshold) then
51: 1100: damage = 1.0_LONGreal
-: 1101: end if
-: 1102:!
-: 1103: end if
-: 1104:!
-: 1105: end if
-: 1106:!
585.9M: 1107: end subroutine perdida
-: 1108:
-: 1109:
-: 1110:!---------------------------- FUNCTION GENERALIZED_HOOKES_LAW -------------------------------
-: 1111:
653.0M: 1112: function generalized_hookes_law (strain_tensor, lambda, mu) result (stress_tensor)
-: 1113:!
-: 1114:! Author: Dr. John K. Prentice
-: 1115:! Affiliation: Quetzal Computational Associates, Inc.
-: 1116:! Dates: 28 November 1997
-: 1117:!
-: 1118:! Purpose: Apply the generalized Hooke's law for elasticity to the strain tensor
-: 1119:! (or strain rate tensor) to compute the stress tensor (or stress rate
-: 1120:! tensor)
-: 1121:!
-: 1122:!############################################################################################
-: 1123:!
-: 1124:! Input:
-: 1125:!
-: 1126:! strain_tensor [selected_real_kind(15,90), dimension(3,3)]
-: 1127:! stress tensor
-: 1128:!
-: 1129:! lambda [selected_real_kind(15,90)]
-: 1130:! Lame constant Lambda
-: 1131:!
-: 1132:! mu [selected_real_kind(15,90)]
-: 1133:! Lame constant mu
-: 1134:!
-: 1135:! Output:
-: 1136:!
-: 1137:! stress_tensor [selected_real_kind(15,90), dimension(3,3)]
-: 1138:! stress tensor
-: 1139:!
-: 1140:!############################################################################################
-: 1141:!
-: 1142:!
-: 1143:!=========== formal variables =============
-: 1144:!
-: 1145: real (kind = LONGreal), dimension(:,:), intent(in) :: strain_tensor
-: 1146: real (kind = LONGreal), intent(in) :: lambda, mu
-: 1147: real (kind = LONGreal), dimension(3,3) :: stress_tensor
-: 1148:!
-: 1149:!========== internal variables ============
-: 1150:!
-: 1151: real (kind = LONGreal), dimension(6) ::generalized_strain_vector, &
-: 1152: generalized_stress_vector
-: 1153: real (kind = LONGreal), dimension(6,6) :: generalized_constitutive_tensor
-: 1154: integer :: i
-: 1155:!
-: 1156:! construct the generalized constitutive tensor for elasticity
-: 1157:!
653.0M: 1158: generalized_constitutive_tensor(:,:) = 0.0_LONGreal
653.0M: 1159: generalized_constitutive_tensor(1,1) = lambda + 2.0_LONGreal * mu
653.0M: 1160: generalized_constitutive_tensor(1,2) = lambda
653.0M: 1161: generalized_constitutive_tensor(1,3) = lambda
653.0M: 1162: generalized_constitutive_tensor(2,1) = lambda
653.0M: 1163: generalized_constitutive_tensor(2,2) = lambda + 2.0_LONGreal * mu
653.0M: 1164: generalized_constitutive_tensor(2,3) = lambda
653.0M: 1165: generalized_constitutive_tensor(3,1) = lambda
653.0M: 1166: generalized_constitutive_tensor(3,2) = lambda
653.0M: 1167: generalized_constitutive_tensor(3,3) = lambda + 2.0_LONGreal * mu
653.0M: 1168: generalized_constitutive_tensor(4,4) = mu
653.0M: 1169: generalized_constitutive_tensor(5,5) = mu
653.0M: 1170: generalized_constitutive_tensor(6,6) = mu
-: 1171:!
-: 1172:! construct the generalized strain vector (using double index notation)
-: 1173:!
653.0M: 1174: generalized_strain_vector(1) = strain_tensor(1,1)
653.0M: 1175: generalized_strain_vector(2) = strain_tensor(2,2)
653.0M: 1176: generalized_strain_vector(3) = strain_tensor(3,3)
653.0M: 1177: generalized_strain_vector(4) = strain_tensor(2,3)
653.0M: 1178: generalized_strain_vector(5) = strain_tensor(1,3)
653.0M: 1179: generalized_strain_vector(6) = strain_tensor(1,2)
-: 1180:!
-: 1181:! compute the generalized stress vector
-: 1182:!
4.6G: 1183: do i = 1, 6
-: 1184: generalized_stress_vector(i) = dot_product(generalized_constitutive_tensor(i,:), &
28.1G: 1185: generalized_strain_vector(:))
-: 1186: end do
-: 1187:!
-: 1188:! update the stress tensor
-: 1189:!
653.0M: 1190: stress_tensor(1,1) = generalized_stress_vector(1)
653.0M: 1191: stress_tensor(2,2) = generalized_stress_vector(2)
653.0M: 1192: stress_tensor(3,3) = generalized_stress_vector(3)
653.0M: 1193: stress_tensor(2,3) = generalized_stress_vector(4)
653.0M: 1194: stress_tensor(1,3) = generalized_stress_vector(5)
653.0M: 1195: stress_tensor(1,2) = generalized_stress_vector(6)
653.0M: 1196: stress_tensor(3,2) = stress_tensor(2,3)
653.0M: 1197: stress_tensor(3,1) = stress_tensor(1,3)
653.0M: 1198: stress_tensor(2,1) = stress_tensor(1,2)
-: 1199:!
653.0M: 1200: end function generalized_hookes_law
-: 1201:
-: 1202:!--------------------------- FUNCTION DAMAGE_RATE --------------------------------------------
-: 1203:
-: 1204: function damage_rate (eta, accumulated_plastic_strain_rate, strain_energy_density, &
-: 1205: accumulated_plastic_strain, plastic_strain_threshold) &
-: 1206: result (damage_dot)
-: 1207:!
-: 1208:! Author: Dr. John K. Prentice
-: 1209:! Affiliation: Quetzal Computational Associates, Inc.
-: 1210:! Dates: 28 November 1997
-: 1211:!
-: 1212:! Purpose: Compute the time derivative of the damage. This routine uses the damage
-: 1213:! potential function of Krajcinovic and Lemaitre (1987).
-: 1214:!
-: 1215:!############################################################################################
-: 1216:!
-: 1217:! Input:
-: 1218:!
-: 1219:! eta [real, selected_real_kind(15,90)]
-: 1220:! material specific constant (kg-sec/J)
-: 1221:!
-: 1222:! accumulated_plastic_strain_rate [real, selected_real_kind(15,90)]
-: 1223:! accumulated plastic strain rate (1/sec)
-: 1224:!
-: 1225:! strain_energy_density [real, selected_real_kind(15,90)]
-: 1226:! strain energy density
-: 1227:!
-: 1228:! accumulated_plastic_strain [real, selected_real_kind(15,90)]
-: 1229:! accumulated plastic strain (dimensionless)
-: 1230:!
-: 1231:! plastic_strain_threshold [real, selected_real_kind(15,90)]
-: 1232:! plastic strain threshold. Damage will nucleate only
-: 1233:! if the accumulated plastic strain exceeds this
-: 1234:! threshold (dimensionless)
-: 1235:!
-: 1236:!
-: 1237:! Output:
-: 1238:!
-: 1239:! damage_dot [real, selected_real_kind(15,90)]
-: 1240:! damage rate (1/sec)
-: 1241:!
-: 1242:!############################################################################################
-: 1243:!
-: 1244:!
-: 1245:!=========== formal variables =============
-: 1246:!
-: 1247: real (kind = LONGreal), intent(in) :: accumulated_plastic_strain, &
-: 1248: accumulated_plastic_strain_rate, &
-: 1249: plastic_strain_threshold, eta, &
-: 1250: strain_energy_density
-: 1251: real (kind = LONGreal) :: damage_dot
-: 1252:!
-: 1253:!========== internal variables ============
-: 1254:!
29.5M: 1255: if (accumulated_plastic_strain < plastic_strain_threshold) then
-: 1256: damage_dot = 0.0_LONGreal
-: 1257: else
29.3M: 1258: damage_dot = strain_energy_density / eta * accumulated_plastic_strain_rate
-: 1259: end if
-: 1260:!
-: 1261: end function damage_rate
-: 1262:
-: 1263:end module perdida_m
-: 1264:
1: 1265: program iztaccihuatl
-: 1266:!
-: 1267:!############################################################################################
-: 1268:!
-: 1269:! Author: Dr. John K. Prentice
-: 1270:! Affiliation: Quetzal Computational Associates, Inc.
-: 1271:! Dates: 28 November 1997
-: 1272:!
-: 1273:! Purpose: Driver for modeling ductile metal fatigue.
-: 1274:!
-: 1275:! Copyright (C) 1997 by Quetzal Computational Associates, Inc.
-: 1276:! Address any questions about this copyright notice to:
-: 1277:!
-: 1278:! Dr. John K. Prentice
-: 1279:! Quetzal Computational Associates, Incorporated
-: 1280:! 3455 Main Avenue, Suite 4
-: 1281:! Durango, CO 81301
-: 1282:! USA
-: 1283:!
-: 1284:! Phone: 970-382-8979
-: 1285:! e-mail: quetzal@quetzalcoatl.com
-: 1286:!
-: 1287:!############################################################################################
-: 1288:!
1: 1289: use read_input_m
-: 1290: use perdida_m
-: 1291: use computer_time_m
-: 1292:!
-: 1293: implicit none
-: 1294:!
-: 1295: integer, parameter :: LONGreal = selected_real_kind(15,90)
-: 1296: real (kind = LONGreal), parameter :: pi = 3.1415926535_LONGreal
-: 1297: real (kind = LONGreal), parameter :: twopi = 2.0_LONGreal * pi
-: 1298:!
-: 1299: real (kind = LONGreal) :: dt, yield_stress, lambda, mu, R_infinity, b, X_infinity, &
-: 1300: gamma, eta, plastic_strain_threshold, wire_radius, &
-: 1301: coil_radius, coil_pitch, radius_of_curvature, time, &
-: 1302: failure_threshold, dr, Poissons_ratio, bulk_modulus, &
-: 1303: coefficient, total_damage, total_surface_area, tnow, &
-: 1304: crack_closure_parameter
-: 1305: real (kind = LONGreal), dimension(:,:,:), allocatable :: strain_rate_tensor, &
-: 1306: strain_tensor, &
-: 1307: stress_tensor, &
-: 1308: plastic_strain_tensor, &
-: 1309: back_stress_tensor
-: 1310: real (kind = LONGreal), dimension(:), allocatable :: damage, radius, &
-: 1311: accumulated_plastic_strain, &
-: 1312: isotropic_hardening_stress, &
-: 1313: max_strain, max_stress, &
-: 1314: surface_area
-: 1315: integer :: number_of_sample_points, spin_frequency, spin_cycles, cmax, n, &
-: 1316: time_steps_per_spin_cycle, next_dump, iteration_count, &
-: 1317: spin_cycle_dump_interval, nn
-: 1318: character (len=80) :: input_file, output_file, file_name
-: 1319: character (len=1) :: answer
-: 1320: logical :: wire_test, coil_test, exists, failure, print_max, dump_history
-: 1321: logical, dimension(:), allocatable :: failed
-: 1322:!
1: 1323: input_file = "fatigue.in"
-: 1324:!
1: 1325: inquire (file=input_file, exist=exists)
1: 1326: if (.not. exists) then
#####: 1327: print *," "
#####: 1328: print *,"Input file does not exist, abort."
#####: 1329: print *," "
#####: 1330: stop
-: 1331: end if
-: 1332:!
1: 1333: output_file = "FATIGUE.OUT"
-: 1334:!
1: 1335: open (unit=11,file=output_file,status="unknown",form="formatted",action="write")
-: 1336:!
-: 1337:! read input file
-: 1338:!
-: 1339: call read_input (input_file, dt, cmax, lambda, mu, yield_stress, R_infinity, b, &
-: 1340: X_infinity, gamma, eta, plastic_strain_threshold, failure_threshold, &
-: 1341: spin_frequency, wire_test, coil_test, wire_radius, coil_radius, &
-: 1342: coil_pitch, radius_of_curvature, time_steps_per_spin_cycle, &
-: 1343: number_of_sample_points, dump_history, crack_closure_parameter, &
1: 1344: spin_cycle_dump_interval)
-: 1345:!
-: 1346:! allocate necessary memory
-: 1347:!
1*: 1348: allocate (strain_rate_tensor(3,3,number_of_sample_points))
1*: 1349: allocate (strain_tensor(3,3,number_of_sample_points))
1*: 1350: allocate (stress_tensor(3,3,number_of_sample_points))
1*: 1351: allocate (plastic_strain_tensor(3,3,number_of_sample_points))
1*: 1352: allocate (back_stress_tensor(3,3,number_of_sample_points))
1*: 1353: allocate (damage(number_of_sample_points))
1*: 1354: allocate (accumulated_plastic_strain(number_of_sample_points))
1*: 1355: allocate (isotropic_hardening_stress(number_of_sample_points))
1*: 1356: allocate (failed(number_of_sample_points))
1*: 1357: allocate (radius(number_of_sample_points))
1*: 1358: allocate (max_strain(number_of_sample_points))
1*: 1359: allocate (max_stress(number_of_sample_points))
1*: 1360: allocate (surface_area(number_of_sample_points))
-: 1361:!
-: 1362:! compute the radius from the center of the wire for each of the sample points
-: 1363:!
1: 1364: dr = wire_radius / real(number_of_sample_points, LONGreal)
201: 1365: do n = 1, number_of_sample_points
200: 1366: radius(n) = 0.5_LONGreal * dr + real(n-1,LONGreal) * dr
-: 1367: surface_area(n) = pi * ((radius(n)+0.5_LONGreal*dr)**2 - &
201: 1368: (radius(n)-0.5_LONGreal*dr)**2)
-: 1369: end do
1: 1370: total_surface_area = pi * wire_radius**2
-: 1371:!
-: 1372:! initialize constitutive values
-: 1373:!
2.6k: 1374: strain_tensor = 0.0_LONGreal
2.6k: 1375: strain_rate_tensor = 0.0_LONGreal
2.6k: 1376: stress_tensor = 0.0_LONGreal
2.6k: 1377: plastic_strain_tensor = 0.0_LONGreal
2.6k: 1378: back_stress_tensor = 0.0_LONGreal
201: 1379: damage = 0.0_LONGreal
201: 1380: accumulated_plastic_strain = 0.0_LONGreal
201: 1381: isotropic_hardening_stress = 0.0_LONGreal
201: 1382: failed = .false.
1: 1383: bulk_modulus = (3.0_LONGreal * lambda + 2.0_LONGreal * mu) / 3.0_LONGreal
-: 1384: Poissons_ratio = 0.5_LONGreal * (3.0_LONGreal * bulk_modulus - 2.0_LONGreal * mu) / &
1: 1385: (3.0_LONGreal * bulk_modulus + mu)
-: 1386: next_dump = 1
-: 1387:!
-: 1388:! loop over time steps and perform calculation
-: 1389:!
1: 1390: time = 0.0_LONGreal
1: 1391: spin_cycles = 0
201: 1392: max_strain = -huge(1.0_LONGreal)
201: 1393: max_stress = -huge(1.0_LONGreal)
-: 1394: print_max = .true.
-: 1395: iteration_count = 9999
-: 1396: do
3.3M: 1397: if (iteration_count < time_steps_per_spin_cycle) then
3.2M: 1398: iteration_count = iteration_count + 1
-: 1399: else
-: 1400: iteration_count = 1
65.9k: 1401: spin_cycles = spin_cycles + 1
-: 1402: end if
3.3M: 1403: time = time + dt
-: 1404:!
-: 1405:! define the total strain and strain rate tensors for this time step
-: 1406:!
-: 1407:!########
-: 1408:! In order to remove proprietary information, the calculations done for the coil
-: 1409:! test have been removed. Consequently, this code only works for a wire
-: 1410:! test.
-: 1411:!########
-: 1412:!
3.3M: 1413: if (wire_test) then
662.5M: 1414: do n = 1, number_of_sample_points
-: 1415: coefficient = (radius(n) / radius_of_curvature) * sin(twopi * &
659.2M: 1416: real(spin_frequency,LONGreal) * time)
659.2M: 1417: strain_tensor(1,1,n) = - coefficient * Poissons_ratio
659.2M: 1418: strain_tensor(2,2,n) = - coefficient * Poissons_ratio
659.2M: 1419: strain_tensor(3,3,n) = coefficient
-: 1420:!
-: 1421: coefficient = twopi * real(spin_frequency,LONGreal) * (radius(n) / &
-: 1422: radius_of_curvature) * cos(twopi * &
659.2M: 1423: real(spin_frequency,LONGreal) * time)
-: 1424:
659.2M: 1425: strain_rate_tensor(1,1,n) = -coefficient * Poissons_ratio
659.2M: 1426: strain_rate_tensor(2,2,n) = - coefficient * Poissons_ratio
662.5M: 1427: strain_rate_tensor(3,3,n) = coefficient
-: 1428: end do
#####: 1429: else if (coil_test) then
#####: 1430: print *,"Abort, coil test not supported"
-: 1431: end if
-: 1432:!
-: 1433:! call the damage model to update the constitutive values for this material
-: 1434:!
662.5M: 1435: do n = 1, number_of_sample_points
662.5M: 1436: if (.not. failed(n)) then
-: 1437: call perdida (dt, lambda, mu, yield_stress, R_infinity, b, X_infinity, &
-: 1438: gamma, eta, plastic_strain_threshold, stress_tensor(:,:,n), &
-: 1439: strain_tensor(:,:,n), plastic_strain_tensor(:,:,n), &
-: 1440: strain_rate_tensor(:,:,n), accumulated_plastic_strain(n), &
-: 1441: back_stress_tensor(:,:,n), isotropic_hardening_stress(n), &
585.9M: 1442: damage(n), failure_threshold, crack_closure_parameter)
-: 1443: end if
-: 1444: end do
-: 1445:!
3.3M: 1446: if (spin_cycles == 1 .and. print_max) then
10.1k: 1447: do n = 1, number_of_sample_points
130.0k*: 1448: max_strain(n) = max(max_strain(n), maxval(strain_tensor(:,:,n)))
130.1k*: 1449: max_stress(n) = max(max_stress(n), maxval(stress_tensor(:,:,n)))
-: 1450: end do
3.3M: 1451: else if (print_max) then
-: 1452: print_max = .false.
1: 1453: write (unit=11,fmt="(a)") " "
-: 1454: write (unit=11,fmt="(a)") &
1: 1455: "Maximum strains and stresses imposed on the undamaged sample: "
1: 1456: write (unit=11,fmt="(a)") " "
201: 1457: do n = 1, number_of_sample_points
200: 1458: write (unit=11,fmt="(a,i4,a,1pe14.5,a,e14.5,a,e14.5)") "region ", n, &
401: 1459: " r=",radius(n), " emax=",max_strain(n), " smax=",max_stress(n)
-: 1460: end do
1: 1461: write (unit=11,fmt="(a)") " "
202*: 1462: if (minval(accumulated_plastic_strain(:)) <= 0.0_LONGreal) then
1: 1463: total_damage = 0.0_LONGreal
201: 1464: do n = 1, number_of_sample_points
201: 1465: if (accumulated_plastic_strain(n) > 0.0_LONGreal) then
84: 1466: total_damage = total_damage + surface_area(n)
-: 1467: end if
-: 1468: end do
1: 1469: total_damage = total_damage / total_surface_area
1: 1470: if (total_damage < failure_threshold) then
#####: 1471: write (unit=*,fmt="(a)") "Simulation terminated."
#####: 1472: write (unit=*,fmt="(a)") "This sample will never fracture."
-: 1473: write (unit=*,fmt="(a)") &
#####: 1474: "The imposed stress is below the fatigue limit."
#####: 1475: write (unit=*,fmt="(a)") "See the output file for more information."
#####: 1476: print *," "
#####: 1477: write (unit=11,fmt="(a)") "Simulation terminated."
#####: 1478: write (unit=11,fmt="(a)") "This sample will never fracture."
-: 1479: write (unit=11,fmt="(a)") &
#####: 1480: "The imposed stress is below the fatigue limit."
#####: 1481: do n = 1, number_of_sample_points
#####: 1482: if (accumulated_plastic_strain(n) > 0.0_LONGreal) then
-: 1483: write (unit=11,fmt="(a,1pe15.5)") &
#####: 1484: "Damage can nucleate at radius ",radius(n)
-: 1485: end if
-: 1486: end do
#####: 1487: stop
-: 1488: end if
-: 1489: end if
-: 1490: end if
-: 1491:!
-: 1492:! update the failed array
-: 1493:!
662.5M: 1494: do n = 1, number_of_sample_points
662.5M: 1495: if (.not. failed(n) .and. damage(n) >= 1.0_LONGreal) then
51: 1496: failed(n) = .true.
51: 1497: print *," "
51: 1498: write (unit=*,fmt="(a)") " "
-: 1499: write (unit=*,fmt="(a,i4,a,1pe15.5,a,i9)") &
51: 1500: "-------> Material failure in region ", n," at r= ", &
102: 1501: radius(n)," at cycle ",spin_cycles
51: 1502: write (unit=11,fmt="(a)") " "
-: 1503: write (unit=11,fmt="(a,i4,a,1pe15.5,a,i9)") &
51: 1504: "-------> Material failure in region ", n," at r= ", &
102: 1505: radius(n)," at cycle ",spin_cycles
-: 1506: end if
-: 1507: end do
-: 1508:!
-: 1509:! check the damage to see if all fatigue criteria have been satisfied. If so,
-: 1510:! exit
-: 1511:!
662.5M: 1512: total_damage = sum(surface_area(:) * damage(:)) / total_surface_area
3.3M: 1513: if (total_damage >= failure_threshold) then
-: 1514: failure = .true.
-: 1515: else
-: 1516: failure = .false.
-: 1517: end if
-: 1518:!
-: 1519:! decide whether to take another time step or to stop simulation now
-: 1520:!
1: 1521: if (failure) then
1: 1522: write (unit=*,fmt="(a)") " "
1: 1523: write (unit=*,fmt="(a)") "Simulation terminated."
1: 1524: write (unit=*,fmt="(a,i10)") " Spin cycle = ",spin_cycles
1: 1525: write (unit=*,fmt="(a,1pe15.5)") " Time (sec) = ", time
1: 1526: write (unit=11,fmt="(a)") " "
1: 1527: write (unit=11,fmt="(a)") "Simulation terminated."
1: 1528: write (unit=11,fmt="(a,i10)") " Spin cycle = ",spin_cycles
1: 1529: write (unit=11,fmt="(a,1pe15.5)") " Time (sec) = ", time
1: 1530: if (wire_test) then
-: 1531: write (unit=*,fmt="(a,1pe15.5)") &
1: 1532: "The wire has fractured, total damage = ", total_damage
-: 1533: write (unit=11,fmt="(a,1pe15.5)") &
1: 1534: "The wire has fractured, total damage = ", total_damage
#####: 1535: else if (coil_test) then
-: 1536: write (unit=11,fmt="(a,1pe15.5)") &
#####: 1537: "The coil has fractured, total damage = ", total_damage
#####: 1538: print *,"The coil has fractured, total damage = ", total_damage
-: 1539: end if
1: 1540: write (unit=11,fmt="(a)") " "
1: 1541: write (unit=11,fmt="(a)") "Final damage results:"
201: 1542: do n = 1, number_of_sample_points
200: 1543: write (unit=11,fmt="(a,1pe14.5,a,e14.5,a,e14.5)") " r = ",radius(n), &
401: 1544: " p = ", accumulated_plastic_strain(n)," D = ", damage(n)
-: 1545: end do
1: 1546: write (unit=11,fmt="(a)") " "
1: 1547: print *," "
1: 1548: call computer_time (tnow)
1: 1549: write (unit=*,fmt="(a,1pe15.5)") "================= Elasped CPU time = ", tnow
1: 1550: write (unit=*,fmt="(a)") " "
1: 1551: write (unit=11,fmt="(a,1pe15.5)") "================= Elasped CPU time = ", tnow
1: 1552: write (unit=11,fmt="(a)") " "
1: 1553: exit
3.3M*: 1554: else if (spin_cycles >= cmax .and. &
-: 1555: iteration_count == time_steps_per_spin_cycle) then
#####: 1556: write (unit=11,fmt="(a)") " "
#####: 1557: write (unit=11,fmt="(a)") "Simulation terminated."
#####: 1558: write (unit=11,fmt="(a,i10)") " Spin cycle = ",spin_cycles
#####: 1559: write (unit=11,fmt="(a,1pe15.5)") " Time (sec) = ", time
-: 1560: write (unit=11,fmt="(a)") &
#####: 1561: "The maximum number of simulation spin cycles was exceeded."
#####: 1562: write (unit=*,fmt="(a)") " "
#####: 1563: write (unit=*,fmt="(a)") "Simulation terminated."
#####: 1564: write (unit=*,fmt="(a,i10)") " Spin cycle = ",spin_cycles
#####: 1565: write (unit=*,fmt="(a,1pe15.5)") " Time (sec) = ", time
-: 1566: write (unit=*,fmt="(a)") &
#####: 1567: "The maximum number of simulation spin cycles was exceeded."
#####: 1568: write (unit=11,fmt="(a)") "Final damage results:"
#####: 1569: do n = 1, number_of_sample_points
#####: 1570: write (unit=11,fmt="(a,1pe14.5,a,e14.5,a,e14.5)") " r = ",radius(n), &
#####: 1571: " p = ", accumulated_plastic_strain(n)," D = ", damage(n)
-: 1572: end do
#####: 1573: write (unit=11,fmt="(a,f6.2)") "Total damage = ",total_damage
#####: 1574: write (unit=11,fmt="(a)") " "
#####: 1575: print *," "
#####: 1576: call computer_time (tnow)
#####: 1577: write (unit=*,fmt="(a,1pe15.5)") "================= Elasped CPU time = ", tnow
#####: 1578: write (unit=*,fmt="(a)") " "
#####: 1579: write (unit=11,fmt="(a,1pe15.5)") "================= Elasped CPU time = ", tnow
#####: 1580: write (unit=11,fmt="(a)") " "
#####: 1581: exit
3.3M: 1582: else if (spin_cycles >= next_dump .and. &
-: 1583: iteration_count == time_steps_per_spin_cycle) then
-: 1584:!
1: 1585: if (next_dump == 1 .and. spin_cycle_dump_interval > 1) then
-: 1586: next_dump = spin_cycle_dump_interval
-: 1587: else
#####: 1588: next_dump = spin_cycles + spin_cycle_dump_interval
-: 1589: end if
-: 1590:!
1: 1591: write (unit=11,fmt="(a)") " "
1: 1592: write (unit=11,fmt="(a,i10)") "================= Spin cycle number = ", &
2: 1593: spin_cycles
1: 1594: write (unit=11,fmt="(a,1pe15.5)") "================= Simulation time = ", &
2: 1595: time
1: 1596: write (unit=*,fmt="(a)") " "
1: 1597: write (unit=*,fmt="(a,i10)") "================= Spin cycle number = ", &
2: 1598: spin_cycles
1: 1599: write (unit=*,fmt="(a,1pe15.5)") "================= Simulation time = ", &
2: 1600: time
1: 1601: call computer_time (tnow)
1: 1602: write (unit=*,fmt="(a,1pe15.5)") "================= Elasped CPU time = ", tnow
1: 1603: write (unit=11,fmt="(a,1pe15.5)") "================= Elasped CPU time = ", tnow
1: 1604: write (unit=11,fmt="(a)") "================= Current results:"
1: 1605: print *," "
1: 1606: write (unit=11,fmt="(a)") " "
202*: 1607: if (maxval(damage(:)) > 0.0_LONGreal) then
#####: 1608: do n = 1, number_of_sample_points
#####: 1609: if (accumulated_plastic_strain(n) > 0.0_LONGreal) then
#####: 1610: write (unit=11,fmt="(a,1pe14.5,a,e14.5,a,e14.5)") " r = ", &
#####: 1611: radius(n)," p = ", accumulated_plastic_strain(n), &
#####: 1612: " D = ",damage(n)
-: 1613: end if
-: 1614: end do
1: 1615: else if (spin_cycles /= 1) then
#####: 1616: write (unit=11,fmt="(a)") "No damage has nucleated at this time."
-: 1617: end if
1: 1618: write (unit=11,fmt="(a)") " "
-: 1619:!
-: 1620:! dump history points, if necessary
-: 1621:!
1: 1622: if (dump_history) then
#####: 1623: if (spin_cycles < 10) then
#####: 1624: write(unit=file_name,fmt="(a,i1,a)") "cycle_",spin_cycles,".dat"
#####: 1625: else if (spin_cycles < 100) then
#####: 1626: write(unit=file_name,fmt="(a,i2,a)") "cycle_",spin_cycles,".dat"
#####: 1627: else if (spin_cycles < 1000) then
#####: 1628: write(unit=file_name,fmt="(a,i3,a)") "cycle_",spin_cycles,".dat"
#####: 1629: else if (spin_cycles < 10000) then
#####: 1630: write(unit=file_name,fmt="(a,i4,a)") "cycle_",spin_cycles,".dat"
#####: 1631: else if (spin_cycles < 100000) then
#####: 1632: write(unit=file_name,fmt="(a,i5,a)") "cycle_",spin_cycles,".dat"
#####: 1633: else if (spin_cycles < 1000000) then
#####: 1634: write(unit=file_name,fmt="(a,i6,a)") "cycle_",spin_cycles,".dat"
#####: 1635: else if (spin_cycles < 10000000) then
#####: 1636: write(unit=file_name,fmt="(a,i7,a)") "cycle_",spin_cycles,".dat"
#####: 1637: else if (spin_cycles < 100000000) then
#####: 1638: write(unit=file_name,fmt="(a,i8,a)") "cycle_",spin_cycles,".dat"
#####: 1639: else if (spin_cycles < 1000000000) then
#####: 1640: write(unit=file_name,fmt="(a,i9,a)") "cycle_",spin_cycles,".dat"
-: 1641: else
#####: 1642: print *," "
#####: 1643: print *,"Warning: Number of spin cycles is too large to continue "
#####: 1644: print *," to dump history files."
#####: 1645: print *," "
#####: 1646: dump_history = .false.
-: 1647: end if
-: 1648: end if
-: 1649:!
1: 1650: if (dump_history) then
-: 1651: open (unit=12,file=file_name,status="unknown",form="formatted", &
#####: 1652: action="write")
#####: 1653: do nn = 1, number_of_sample_points
#####: 1654: write (unit=12,fmt="(i4,2x,1p4e15.5)") nn, radius(nn), &
#####: 1655: accumulated_plastic_strain(nn), damage(nn)
-: 1656: end do
#####: 1657: close (unit = 12)
-: 1658: end if
-: 1659: end if
-: 1660: end do
-: 1661:!
-: 1662:! deallocate memory
-: 1663:!
1*: 1664: deallocate (strain_rate_tensor)
1*: 1665: deallocate (strain_tensor)
1*: 1666: deallocate (stress_tensor)
1*: 1667: deallocate (plastic_strain_tensor)
1*: 1668: deallocate (back_stress_tensor)
1*: 1669: deallocate (damage)
1*: 1670: deallocate (accumulated_plastic_strain)
1*: 1671: deallocate (isotropic_hardening_stress)
1*: 1672: deallocate (failed)
1*: 1673: deallocate (max_strain)
1*: 1674: deallocate (max_stress)
1*: 1675: deallocate (surface_area)
-: 1676:!
1: 1677: close (unit = 11)
-: 1678:!
1: 1679: end program iztaccihuatl
-: 1680: