MultiPhysicsVault/.raw/FiniteElementsinPlasticityTheoryandPractice/FiniteElementsinPlasticityTheoryandPractice_030.md

\boldsymbol {M} _ {1} = \left[ \begin{array}{c c c c} \frac {2}{3} & - \frac {1}{3} & 0 & - \frac {1}{3} \\ & \frac {2}{3} & 0 & - \frac {1}{3} \\ \text { Symmetric } & 2 & 0 \\ & & & \frac {2}{3} \end{array} \right], \tag {8.44}

and

\boldsymbol {M} _ {2} = \left[ \begin{array}{c c c c} \left(\sigma_ {x} ^ {\prime}\right) ^ {2} & \sigma_ {x} ^ {\prime} \sigma_ {y} ^ {\prime} & 2 \sigma_ {x} ^ {\prime} \tau_ {x y} & \sigma_ {x} ^ {\prime} \sigma_ {z} ^ {\prime} \\ & \left(\sigma_ {y} ^ {\prime}\right) ^ {2} & 2 \sigma_ {y} ^ {\prime} \tau_ {x y} & \sigma_ {y} ^ {\prime} \sigma_ {z} ^ {\prime} \\ \text { Symmetric } & & 4 \left(\tau_ {x y}\right) ^ {2} & 2 \tau_ {x y} \sigma_ {z} ^ {\prime} \\ - & - & - & - & - \\ & & & & \left(\sigma_ {z} ^ {\prime}\right) ^ {2} \end{array} \right], \tag {8.45}

and J_{2}' is given by (7.76). For plane stress and plane strain problems only the upper 3 \times 3 partition is employed while for axisymmetric situations the complete matrices are utilised with x and y being replaced by r and z respectively.

Similar expressions can be derived for the Tresca, Mohr–Coulomb and Drucker–Prager yield criteria by employing the appropriate expression for F in (8.36) and repeating the above calculations. The form of F is given in (7.63), (7.65) and (7.66) for the Tresca, Mohr–Coulomb and Drucker–Prager laws respectively.

8.6 Program structure

The computation sequence for the program is shown in Fig. 8.1. The program structure follows closely that for static elasto-plastic analysis described in Chapter 7. In fact, the majority of the subroutines utilised are common to both applications and it is only the additional subroutines required that are described in this chapter. For the viscoplastic program the time stepping loop replaces the nonlinear solution iteration loop for conventional plasticity and subroutine STEPVP, whose main role is to evaluate quantities at the end of a timestep, replaces the plasticity subroutine RESIDU. In this chapter we need to describe in detail subroutines STIFVP, TANGVP, STEPVP, FLOWVP and STEADY. The descriptions of all other subroutines required for assembly of a working viscoplastic program have been given in Chapters 6 and 7. The version described is restricted to the case of infinitesimal strains. The modifications required to include large deformation effects are straightforward and are left as an exercise to the reader. Furthermore, for implicit schemes, only the Von Mises yield criterion is considered.

The list of material properties accepted in subroutine INPUT described in Section 6.5.1 must be extended beyond those required for elasto-plastic analysis, since additional material parameters are required to define the

flowchart

graph TD
    A["START"] --> B["DIMEN<br>Presents the variables associated with the dynamic dimensioning process"]
    B --> C["INPUT<br>Inputs data defining geometry, boundary conditions and material properties"]
    C --> D["LOADPS<br>Evaluates the equivalent nodal forces for pressure loading, gravity loading, etc."]
    D --> E["ZERO<br>Sets to zero arrays required for accumulation of data"]
    E --> F["INCREM<br>Increments the applied loads according to specified load factors"]
    F --> G["STIFVP<br>Calculates the element stiffnesses as K_T^n(σ^n) (Eq. (8.24))"]
    F --> H["TANGVP<br>Evaluates D̂^n according to (8.18)"]
    G --> I["FRONT<br>Solves the simultaneous equation system by the frontal method, i.e. Δd^n = [K_T^n"]^-1ΔV^n d^{n+1} = d^n + Δd^n]
    H --> I
    I --> J["STEPVP<br>Evaluates quantities at the end of the timestep<br>a) Δσ^n = D̂^n(B^nΔd - ε̇_vp^nΔt_n) b) σ^{n+1} = σ^n + Δσ^n<br>c) ε̇_vp^{n+1} = ε̇_vp^n + ε̇_vp^nΔt_n d) Δt_{n+1}"]
    I --> K["INVAR<br>Evaluates the effective stress level"]
    I --> L["YIELDF & FLOWVP<br>Determines:—<br>a) The flow vector, a<br>b) ε̇_vp^{n+1} = γ(Φ)a^{n+1}"]
    J --> M["OUTPUT<br>Prints the results for the current timestep"]
    L --> M
    M --> N["END"]
    O["LOAD INCREMENT LOOP"] --> F
    P["TIME STEPPING LOOP"] --> F
    Q["END"] --> M

Fig. 8.1 Flow sequence for the two-dimensional elasto-viscoplastic stress analysis program.

viscoplastic flow. This is accomplished by specifying the value of NPROP as 10 in subroutine DIMEN, described in Section 7.8.1, and inputting the following properties for each different material.

PROPS(NUMAT, 1)	Elastic modulus, $E$ .
PROPS(NUMAT, 2)	Poissons ratio, $\nu$ .
PROPS(NUMAT, 3)	Material thickness, $t$ .
PROPS(NUMAT, 4)	Material mass density, $\rho$ .
PROPS(NUMAT, 5)	Uniaxial yield stress $\sigma_{Y}$ (Tresca and Von Mises solids); Cohesion $c$ (Mohr–Coulomb and Drucker–Prager materials).
PROPS(NUMAT, 6)	Hardening parameter $H'$ for linear strain hardening.
PROPS(NUMAT, 7)	Angle of internal friction for Mohr–Coulomb and Drucker–Prager materials only.
PROPS(NUMAT, 8)	The fluidity parameter, $\gamma$ .
PROPS(NUMAT, 9)	The coefficient $M$ in (8.8) or coefficient $N$ in (8.9).
PROPS(NUMAT, 10)	Indicator specifying type of flow function to be employed:0 – Flow function (8.8)1 – Flow function (8.9)

8.7 Formulation of the tangential stiffness matrix

The role of the subroutines described in this section is to calculate the tangential stiffness matrix for each element according to (8.24) . The complete operation is shared between three subroutines which will now be described.

8.7.1 Subroutine STIFVP

This subroutine controls the overall formulation of the tangential stiffness matrix for each element and is very similar to subroutine STIFFP, described in Section 7.8.5, which performs the same task for conventional plasticity. For the case of small deformations, matrix B^{n} is constant and equal to B_{0} the usual infinitesimal elastic value. Matrix B_{0} is given by subroutine BMATPS described in Section 6.4.7. To evaluate K_{T^{n}} it is necessary to find \hat{D}^{n} whose precise form is given by (8.18). With the normal elastic material matrix D replaced by \hat{D}^{n} , the stiffness evaluation follows the standard procedure described in Section 7.8.5. Subroutine STIFVP can now be presented and described.

SUBROUTINE STIFVP(COORD,IINCS,LNODS,MATNO,MEVAB,MMATS,	STVP	1
MPOIN,MTOTV,NELEM,NEVAB,NGAUS,NNODE,NSTRE,	STVP	2
NSTR1,POSGP,PROPS,WEIGP,MELEM,MTOTG,	STVP	3
STRSG,NTYPE,NCRIT,TIMEX,DTIME)	STVP	4
C***************	STVP	5
C	STVP	6
C**** THIS SUBROUTINE EVALUATES THE STIFFNESS MATRIX FOR EACH ELEMENT	STVP	7
C IN TURN	STVP	8
C	STVP	9
C***************	STVP	10
DIMENSION BMATX(4,18),CARTD(2,9),COORD(MPOIN,2),DBMAT(4,18),	STVP	11
DERIV(2,9),DEVIA(4),DMATX(4,4),	STVP	12
ELCOD(2,9),EPSTN(MTOTG),ESTIF(18,18),LNODS(MELEM,9),	STVP	13
MATNO(MELEM),POSGP(4),PROPS(MMATS,10),SHAPE(9),	STVP	14
WEIGP(4),STRES(4),STRSG(4,MTOTG),	STVP	15

DVECT(4), AVECT(4), GPCOD(2,9)
TWOPI=6.283185308
REWIND 1
KGAUS=0
STVP 16
STVP 17
STVP 18
STVP 19
STVP 20
STVP 21
STVP 22
STVP 23
STVP 24
STVP 25
STVP 26
STVP 27
STVP 28
STVP 29
STVP 30
STVP 31
STVP 32
STVP 33
STVP 34
STVP 35
STVP 36
STVP 37
STVP 38
STVP 39
STVP 40
STVP 41
STVP 42
STVP 43
STVP 44
STVP 45
STVP 46
STVP 47
STVP 48
STVP 49
STVP 50
STVP 51
STVP 52
STVP 53
STVP 54
STVP 55
STVP 56
STVP 57
STVP 58
STVP 59
STVP 60
STVP 61
STVP 62
STVP 63
STVP 64
STVP 65
STVP 66
STVP 67
STVP 68
STVP 69
STVP 70
STVP 71
STVP 72
STVP 73
STVP 74
STVP 75
STVP 76
STVP 77
STVP 78
STVP 79
STVP 80
STVP 81
STVP 82
STVP 83
STVP 84
STVP 85
STVP 86
STVP 87
STVP 88
STVP 89
STVP 90
STVP 91
STVP 92
STVP 93
STVP 94
STVP 95
STVP 96
STVP 97
STVP 98
STVP 99
STVP 100
STVP 101
STVP 102
STVP 103
STVP 104
STVP 105
STVP 106
STVP 107
STVP 108
STVP 109
STVP 110
STVP 111
STVP 112
STVP 113
STVP 114
STVP 115
STVP 116
STVP 117
STVP 118
STVP 119
STVP 120
STVP 121
STVP 122
STVP 123
STVP 124
STVP 125
STVP 126
STVP 127
STVP 128
STVP 129
STVP 130
STVP 131
STVP 132
STVP 133
STVP 134
STVP 135
STVP 136
STVP 137
STVP 138
STVP 139
STVP 140
STVP 141
STVP 142
STVP 143
STVP 144
STVP 145
STVP 146
STVP 147
STVP 148
STVP 149
STVP 150
STVP 151
STVP 152
STVP 153
STVP 154
STVP 155
STVP 156
STVP 157
STVP 158
STVP 159
STVP 160
STVP 161
STVP 162
STVP 163
STVP 164
STVP 165
STVP 166
STVP 167
STVP 168
STVP 169
STVP 170
STVP 171
STVP 172
STVP 173
STVP 174
STVP 175
STVP 176
STVP 177
STVP 178
STVP 179
STVP 180
STVP 181
STVP 182
STVP 183
STVP 184
STVP 185
STVP 186
STVP 187
STVP 188
STVP 189
STVP 190
STVP 191
STVP 192
STVP 193
STVP 194
STVP 195
STVP 196
STVP 197
STVP 198
STVP 199
STVP 200
STVP 201
STVP 202
STVP 203
STVP 204
STVP 205
STVP 206
STVP 207
STVP 208
STVP 209
STVP 210
STVP 211
STVP 212
STVP 213
STVP 214
STVP 215
STVP 216
STVP 217
STVP 218
STVP 219
STVP 220
STVP 221
STVP 222
STVP 223
STVP 224
STVP 225
STVP 226
STVP 227
STVP 228
STVP 229
STVP 230
STVP 231
STVP 232
STVP 233
STVP 234
STVP 235
STVP 236
STVP 237
STVP 238
STVP 239
STVP 240
STVP 241
STVP 242
STVP 243
STVP 244
STVP 245
STVP 246
STVP 247
STVP 248
STVP 249
STVP 250
STVP 251
STVP 252
STVP 253
STVP 254
STVP 255
STVP 256
STVP 257
STVP 258
STVP 259
STVP 260
STVP 261
STVP 262
STVP 263
STVP 264
STVP 265
STVP 266
STVP 267
STVP 268
STVP 269
STVP 270
STVP 271
STVP 272
STVP 273
STVP 274
STVP 275
STVP 276
STVP 277
STVP 278
STVP 279
STVP 280
STVP 281
STVP 282
STVP 283
STVP 284
STVP 285
STVP 286
STVP 287
STVP 288
STVP 289
STVP 290
STVP 291
STVP 292
STVP 293
STVP 294
STVP 295
STVP 296
STVP 297
STVP 298
STVP 299
STVP 300
STVP 301
STVP 302
STVP 303
STVP 304
STVP 305
STVP 306
STVP 307
STVP 308
STVP 309
STVP 310
STVP 311
STVP 312
STVP 313
STVP 314
STVP 315
STVP 316
STVP 317
STVP 318
STVP 319
STVP 320
STVP 321
STVP 322
STVP 323
STVP 324
STVP 325
STVP 326
STVP 327
STVP 328
STVP 329
STVP 330
STVP 331
STVP 332
STVP 333
STVP 334
STVP 335
STVP 336
STVP 337
STVP 338
STVP 339
STVP 340
STVP 341
STVP 342
STVP 343
STVP 344
STVP 345
STVP 346
STVP 347
STVP 348
STVP 349
STVP 350
STVP 351
STVP 352
STVP 353
STVP 354
STVP 355
STVP 356
STVP 357
STVP 358
STVP 359
STVP 360
STVP 361
STVP 362
STVP 363
STVP 364
STVP 365
STVP 366
STVP 367
STVP 368
STVP 369
STVP 370
STVP 371
STVP 372
STVP 373
STVP 374
STVP 375
STVP 376
STVP 377
STVP 378
STVP 379
STVP 380
STVP 381
STVP 382
STVP 383
STVP 384
STVP 385
STVP 386
STVP 387
STVP 388
STVP 389
STVP 390
STVP 391
STVP 392
STVP 393
STVP 394
STVP 395
STVP 396
STVP 397
STVP 398
STVP 399
STVP 400
STVP 401
STVP 402
STVP 403
STVP 404
STVP 405
STVP 406
STVP 407
STVP 408
STVP 409
STVP 410
STVP 411
STVP 412
STVP 413
STVP 414
STVP 415
STVP 416
STVP 417
STVP 418
STVP 419
STVP 420
STVP 421
STVP 422
STVP 423
STVP 424
STVP 425
STVP 426
STVP 427
STVP 428
STVP 429
STVP 430
STVP 431
STVP 432
STVP 433
STVP 434
STVP 435
STVP 436
STVP 437
STVP 438
STVP 439
STVP 440
STVP 441
STVP 442
STVP 443
STVP 444
STVP 445
STVP 446
STVP 447
STVP 448
STVP 449
STVP 450
STVP 451
STVP 452
STVP 453
STVP 454
STVP 455
STVP 456
STVP 457
STVP 458
STVP 459
STVP 460
STVP 461
STVP 462
STVP 463
STVP 464
STVP 465
STVP 466
STVP 467
STVP 468
STVP 469
STVP 470
STVP 471
STVP 472
STVP 473
STVP 474
STVP 475
STVP 476
STVP 477
STVP 478
STVP 479
STVP 480
STVP 481
STVP 482
STVP 483
STVP 484
STVP 485
STVP 486
STVP 487
STVP 488
STVP 489
STVP 490
STVP 491
STVP 492
STVP 493
STVP 494
STVP 495
STVP 496
STVP 497
STVP 498
STVP 499
STVP 500
STVP 501
STVP 502
STVP 503
STVP 504
STVP 505
STVP 506
STVP 507
STVP 508
STVP 509
STVP 510
STVP 511
STVP 512
STVP 513
STVP 514
STVP 515
STVP 516
STVP 517
STVP 518
STVP 519
STVP 520
STVP 521
STVP 522
STVP 523
STVP 524
STVP 525
STVP 526
STVP 527
STVP 528
STVP 529
STVP 530
STVP 531
STVP 532
STVP 533
STVP 534
STVP 535
STVP 536
STVP 537
STVP 538
STVP 539
STVP 540
STVP 541
STVP 542
STVP 543
STVP 544
STVP 545
STVP 546
STVP 547
STVP 548
STVP 549
STVP 550
STVP 551
STVP 552
STVP 553
STVP 554
STVP 555
STVP 556
STVP 557
STVP 558
STVP 559
STVP 560
STVP 561
STVP 562
STVP 563
STVP 564
STVP 565
STVP 566
STVP 567
STVP 568
STVP 569
STVP 570
STVP 571
STVP 572
STVP 573
STVP 574
STVP 575
STVP 576
STVP 577
STVP 578
STVP 579
STVP 580
STVP 581
STVP 582
STVP 583
STVP 584
STVP 585
STVP 586
STVP 587
STVP 588
STVP 589
STVP 590
STVP 591
STVP 592
STVP 593
STVP 594
STVP 595
STVP 596
STVP 597
STVP 598
STVP 599
STVP 600
STVP 601
STVP 602
STVP 603
STVP 604
STVP 605
STVP 606
STVP 607
STVP 608
STVP 609
STVP 610
STVP 611
STVP 612
STVP 613
STVP 614
STVP 615
STVP 616
STVP 617
STVP 618
STVP 619
STVP 620
STVP 621
STVP 622
STVP 623
STVP 624
STVP 625
STVP 626
STVP 627
STVP 628
STVP 629
STVP 630
STVP 631
STVP 632
STVP 633
STVP 634
STVP 635
STVP 636
STVP 637
STVP 638
STVP 639
STVP 640
STVP 641
STVP 642
STVP 643
STVP 644
STVP 645
STVP 646
STVP 647
STVP 648
STVP 649
STVP 650
STVP 651
STVP 652
STVP 653
STVP 654
STVP 655
STVP 656
STVP 657
STVP 658
STVP 659
STVP 660
STVP 661
STVP 662
STVP 663
STVP 664
STVP 665
STVP 666
STVP 667
STVP 668
STVP 669
STVP 670
STVP 671
STVP 672
STVP 673
STVP 674
STVP 675
STVP 676
STVP 677
STVP 678
STVP 679
STVP 680
STVP 681
STVP 682
STVP 683
STVP 684
STVP 685
STVP 686
STVP 687
STVP 688
STVP 689
STVP 690
STVP 691
STVP 692
STVP 693
STVP 694
STVP 695
STVP 696
STVP 697
STVP 698
STVP 699
STVP 700
STVP 701
STVP 702
STVP 703
STVP 704
STVP 705
STVP 706
STVP 707
STVP 708
STVP 709
STVP 710
STVP 711
STVP 712
STVP 713
STVP 714
STVP 715
STVP 716
STVP 717
STVP 718
STVP 719
STVP 720
STVP 721
STVP 722
STVP 723
STVP 724
STVP 725
STVP 726
STVP 727
STVP 728
STVP 729
STVP 730
STVP 731
STVP 732
STVP 733
STVP 734
STVP 735
STVP 736
STVP 737
STVP 738
STVP 739
STVP 740
STVP 741
STVP 742
STVP 743
STVP 744
STVP 745
STVP 746
STVP 747
STVP 748
STVP 749
STVP 750
STVP 751
STVP 752
STVP 753
STVP 754
STVP 755
STVP 756
STVP 757
STVP 758
STVP 759
STVP 760
STVP 761
STVP 762
STVP 763
STVP 764
STVP 765
STVP 766
STVP 767
STVP 768
STVP 769
STVP 770
STVP 771
STVP 772
STVP 773
STVP 774
STVP 775
STVP 776
STVP 777
STVP 778
STVP 779
STVP 780
STVP 781
STVP 782
STVP 783
STVP 784
STVP 785
STVP 786
STVP 787
STVP 788
STVP 789
STVP 790
STVP 791
STVP 792
STVP 793
STVP 794
STVP 795
STVP 796
STVP 797
STVP 798
STVP 799
STVP 800
STVP 801
STVP 802
STVP 803
STVP 804
STVP 805
STVP 806
STVP 807
STVP 808
STVP 809
STVP 810
STVP 811
STVP 812
STVP 813
STVP 814
STVP 815
STVP 816
STVP 817
STVP 818
STVP 819
STVP 820
STVP 821
STVP 822
STVP 823
STVP 824
STVP 825
STVP 826
STVP 827
STVP 828
STVP 829
STVP 830
STVP 831
STVP 832
STVP 833
STVP 834
STVP 835
STVP 836
STVP 837
STVP 838
STVP 839
STVP 840
STVP 841
STVP 842
STVP 843
STVP 844
STVP 845
STVP 846
STVP 847
STVP 848
STVP 849
STVP 850
STVP 851
STVP 852
STVP 853
STVP 854
STVP 855
STVP 856
STVP 857
STVP 858
STVP 859
STVP 860
STVP 861
STVP 862
STVP 863
STVP 864
STVP 865
STVP 866
STVP 867
STVP 868
STVP 869
STVP 870
STVP 871
STVP 872
STVP 873
STVP 874
STVP 875
STVP 876
STVP 877
STVP 878
STVP 879
STVP 880
STVP 881
STVP 882
STVP 883
STVP 884
STVP 885
STVP 886
STVP 887
STVP 888
STVP 889
STVP 890
STVP 891
STVP 892
STVP 893
STVP 894
STVP 895
STVP 896
STVP 897
STVP 898
STVP 899
STVP 900
STVP 901
STVP 902 

50 CONTINUE	STVP	81
C	STVP	82
C*** CONSTRUCT THE LOWER TRIANGLE OF THE STIFFNESS MATRIX	STVP	83
C	STVP	84
DO 60 IEVAB=1,NEVAB	STVP	85
DO 60 JEVAB=1,NEVAB	STVP	86
60 ESTIF(JEVAB,IEVAB)=ESTIF(IEVAB,JEVAB)	STVP	87
C	STVP	88
C*** STORE THE STIFFNESS MATRIX,STRESS MATRIX AND SAMPLING POINT	STVP	89
C COORDINATES FOR EACH ELEMENT ON DISC FILE	STVP	90
C	STVP	91
WRITE(1) ESTIF	STVP	92
70 CONTINUE	STVP	93
RETURN	STVP	94
END	STVP	95

STVP 17	Compute the value of $2\pi$ .
STVP 18	Rewind the disc file on which the element stiffness matrices will be stored in turn.
STVP 19	Set to zero the counter which indicates the overall Gauss point location.
STVP 23	Enter the loop over each element in the structure.
STVP 24	Identify the material property type of the current element.
STVP 28–33	Store the element nodal coordinates in the local array ELCOD for convenient use later.
STVP 34	Identify the element thickness.
STVP 38–40	Zero the element stiffness array.
STVP 41	Set to zero the element Gauss point counter.
STVP 45–48	Enter the numerical integration loops and locate the position $(\xi, \eta)$ of the current point.
STVP 49–50	Increment the local and global Gauss point counters.
STVP 54	Call subroutine MODPS to evaluate the elasticity matrix, $D$ .
STVP 58	Evaluate the shape functions $N_i$ and $\partial N_i/\partial \xi$ , $\partial N_i/\partial \eta$ for the current Gauss point.
STVP 59–60	Evaluate the Gauss point coordinates, GPCOD(IDIME, KGASP), the determinant of the Jacobian matrix $|J|$ and the Cartesian derivatives of the shape functions $\partial N_i/\partial x$ , $\partial N_i/\partial y$ (or $\partial N_i/\partial r$ , $\partial N_i/\partial z$ for axisymmetric problems).
STVP 61–63	Calculate the elemental volume for numerical integration as $|J|W_\xi W_\eta$ taking care to multiply by the appropriate element thickness or by $2\pi r$ for axisymmetric problems.
STVP 67	Evaluate the $B$ matrix.
STVP 68–69	Store the current stresses in a local array.
STVP 70–71	For an implicit or semi-implicit timestepping scheme $(\Theta \neq 0)$ , call subroutine TANGVP to evaluate $\hat{D}^n$ which is stored as DMATX.
STVP 72	Evaluate $DB$ (or $\hat{D}^n B$ for implicit schemes).
STVP 76–80	Compute the upper triangle of the element stiffness matrix as

\int_ {\Omega} \boldsymbol {B} ^ {T} \hat {\boldsymbol {D}} ^ {n} \boldsymbol {B} d \Omega .

STVP 81 End of loop for numerical integration.

STVP 85–87 Complete the lower triangle of the element stiffness matrix by symmetry.

STVP 92 Store the element stiffness matrix on disc file 1.

STVP 93 Return to process the next element.

8.7.2 Subroutine TANGVP

The function of this subroutine is to evaluate \hat{D}^{n} for use in (8.24). Matrix \hat{D}^{n} , which is defined in (8.18), is stress dependent and therefore must be calculated for each Gaussian integrating point in turn. The computational sequence followed is:

a) Evaluate H^n according to (8.42)
b) Calculate C^n according to (8.14)
c) Evaluate \tilde{D}^n according to (8.18)

Two forms of the flow function \Phi are considered as defined in (8.8) and (8.9). Thus, for use in (8.43), we have

\frac {d \Phi}{d F} = \frac {M}{F _ {0}} e ^ {M \left(\frac {F - F _ {0}}{F _ {0}}\right)}

or

\frac {d \Phi}{d F} = \frac {N}{F _ {0}} \left(\frac {F - F _ {0}}{F _ {0}}\right) ^ {N - 1}. \tag {8.46}

Array DMATX which originally contains the elastic matrix D is used to finally store \hat{D}^{n} . The matrix inversions required in (8.18) are performed by a separate subroutine, INVERT.

Subroutine TANGVP is now presented and described.

SUBROUTINE TANGVP(LPROP,STRES,PROPS,TIMEX,DTIME,NTSTRE,NTYPE,MMATS,NCRIT,DMATX) TGVP 1
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
D
D
E
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
C
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
II
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
I
No
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
P
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
V
P
P
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
1
T
G
W
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
No
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
C
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
T
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
F
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
P
N
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
M
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
N
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
S
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
M
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
M
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
S
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
N
s
I
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
N
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
№
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
0
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
Yes
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
NO
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
 No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
S
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
No
# 

10 CMULT=FNORM**DELTA TGVP 24
GRADP=DELTA*(FNORM**(DELTA-1.0))/FDATM TGVP 25
20 FACT1=GAMMA*ROOT3*CMULT/(2.0*STEFF) TGVP 26
FACT2=GAMMA*(0.75*GRADP/VARJ2-3.0*CMULT/(4.0*ROOT3*STEFF*VARJ2)) TGVP 27
C TGVP 28
C*** MATRICES M1 AND M2 FOR A VON MISES MATERIAL TGVP 29
C TGVP 30
TRIX1(1,1)=0.666666667 TGVP 31
TRIX1(1,2)=-0.333333333 TGVP 32
TRIX1(1,3)=0.0 TGVP 33
TRIX1(2,2)=0.666666667 TGVP 34
TRIX1(2,3)=0.0 TGVP 35
TRIX1(3,3)=2.0 TGVP 36
IF(NTYPE.NE.3) GO TO 30 TGVP 37
TRIX1(1,4)=-0.333333333 TGVP 38
TRIX1(2,4)=-0.333333333 TGVP 39
TRIX1(3,4)=0.0 TGVP 40
TRIX1(4,4)=0.666666667 TGVP 41
30 TRIX2(1,1)=DEVIA(1)*DEVIA(1) TGVP 42
TRIX2(1,2)=DEVIA(1)*DEVIA(2) TGVP 43
TRIX2(1,3)=2.0*DEVIA(1)*DEVIA(3) TGVP 44
TRIX2(2,2)=DEVIA(2)*DEVIA(2) TGVP 45
TRIX2(2,3)=2.0*DEVIA(2)*DEVIA(3) TGVP 46
TRIX2(3,3)=4.0*DEVIA(3)*DEVIA(3) TGVP 47
IF(NTYPE.NE.3) GO TO 40 TGVP 48
TRIX2(1,4)=DEVIA(1)*DEVIA(4) TGVP 49
TRIX2(2,4)=DEVIA(2)*DEVIA(4) TGVP 50
TRIX2(3,4)=2.0*DEVIA(3)*DEVIA(4) TGVP 51
TRIX2(4,4)=DEVIA(4)*DEVIA(4) TGVP 52
40 DO 50 ISTRE=1,NSTRE TGVP 53
DO 50 JSTRE=1,NSTRE TGVP 54
TRIX1(JSTRE,ISTRE)=TRIX1(ISTRE,JSTRE) TGVP 55
50 TRIX2(JSTRE,ISTRE)=TRIX2(ISTRE,JSTRE) TGVP 56
DO 60 ISTRE=1,NSTRE TGVP 57
DO 60 JSTRE=1,NSTRE TGVP 58
60 CMATX(ISTRE,JSTRE)=TIMEX*DTIME*(FACT1*TRIX1(ISTRE,JSTRE) TGVP 59
.+FACT2*TRIX2(ISTRE,JSTRE)) TGVP 60
CALL INVERT(DMATX,TMATX,NSTRE) TGVP 61
DO 70 ISTRE=1,NSTRE TGVP 62
DO 70 JSTRE=1,NSTRE TGVP 63
70 TMATX(ISTRE,JSTRE)=TMATX(ISTRE,JSTRE)+CMATX(ISTRE,JSTRE) TGVP 64
CALL INVERT(TMATX,DMATX,NSTRE) TGVP 65
RETURN TGVP 66
END TGVP 67 

TGVP 10 Evaluate \sqrt{(3)} .

TGVP 11 Identify the yield stress F as FDATM.

TGVP 12 Identify the fluidity parameter \gamma as GAMMA.

TGVP 13 For flow law (8.8) store the index M as DELTA, or for flow law (8.9) store the index N as DELTA.

TGVP 14 Identify the type of flow function to be used as governed by material property PROPS(LPROP,10) supplied as input:

NFLOW = 0 - Flow function (8.8) to be used,

NFLOW = 1 - Flow function (8.9) to be used.

TGVP 15–16 Call subroutine INVAR to evaluate the effective stress components, the effective stress level and J_{2}' .

TGVP 17-18 Evaluate F-F_{0}/F_{0} as FNORM.

TGVP 21-22 Evaluate \Phi and d\Phi / dF for flow function (8.8).
TGVP 24–25 Evaluate \Phi and d\Phi/dF for flow function (8.9).
TGVP 26–27 Compute p_{1} and p_{2} according to (8.43).
TGVP 31–41 Evaluate M_{1} according to (8.44) taking the full 4 \times 4 matrix for axisymmetric situations.
TGVP 42–52 Evaluate M_{2} according to (8.45) taking the full 4 \times 4 matrix for axisymmetric situations.
TGVP 53–56 Complete the lower triangle of M_{1} and M_{2} by symmetry.
TGVP 57–60 Compute matrix C^{n} according to (8.14) and (8.42).
TGVP 61 Call subroutine INVERT to evaluate D^{-1} and store as TMATX.
TGVP 62-64 Compute D^{-1} + C^n
TGVP 65 Call subroutine INVERT to evaluate (D^{-1}+C^{n})^{-1} and store as DMATX.

8.7.3 Subroutine INVERT

The function of this subroutine is to determine the inverse of any arbitrary square matrix. In particular, the subroutine accepts a matrix AMATX with dimensions NARAY×NARAY and evaluates the inverse as BMATX. The procedure employed is the standard method of reduction in which starting from the original matrix AMATX and assuming an identity matrix for BMATX, an elimination process is followed until AMATX is reduced to an identity form. Then at this stage BMATX is the inverse of AMATX.

The subroutine is presented below without further comment.

SUBROUTINE INVERT(AMATX,BMATX,NARAY) INVT 1
C*******************************
C INVT 2
C INVT 3
C*** TO PROVIDE THE INVERSE OF AMATX AS BMATX INVT 4
C INVT 5
C*******************************
INVT 6
DIMENSION AMATX(4,4),BMATX(4,4) INVT 7
DO 10 IARAY=1,NARAY INVT 8
DO 10 JARAY=1,NARAY INVT 9
BMATX(IARAY,JARAY)=0.0 INVT 10
10 IF(IARAY.EQ.JARAY) BMATX(IARAY,JARAY)=1.0 INVT 11
DO 20 IARAY=1,NARAY INVT 12
DENOM=AMATX(IARAY,IARAY) INVT 13
DO 30 JARAY=1,NARAY INVT 14
AMATX(IARAY,JARAY)=AMATX(IARAY,JARAY)/DENOM INVT 15
30 BMATX(IARAY,JARAY)=BMATX(IARAY,JARAY)/DENOM INVT 16
KARAY=IARAY+1 INVT 17
IF(KARAY.GT.NARAY) GO TO 40 INVT 18
DO 20 JARAY=KARAY,NARAY INVT 19
CONST=AMATX(JARAY,IARAY) INVT 20
DO 20 LARAY=IARAY,NARAY INVT 21
AMATX(JARAY,LARAY)=AMATX(JARAY,LARAY)-AMATX(IARAY,LARAY) INVT 22
.*CONST INVT 23
20 BMATX(JARAY,LARAY)=BMATX(JARAY,LARAY)-BMATX(IARAY,LARAY) INVT 24
.*CONST INVT 25
40 CONTINUE INVT 26
DO 50 IARAY=2,NARAY INVT 27
KARAY=NARAY-IARAY+2 INVT 28 

LIMIT=KARAY-1 INVT 29 DO 50 LARAY=1,LIMIT INVT 30 CONST=AMATX(LARAY,KARAY) INVT 31 DO 50 JARAY=1,KARAY INVT 32 AMATX(LARAY,JARAY)=AMATX(LARAY,JARAY)-AMATX(KARAY,JARAY) INVT 33 .*CONST INVT 34 50 BMATX(LARAY,JARAY)=BMATX(LARAY,JARAY)-BMATX(KARAY,JARAY) INVT 35 .*CONST INVT 36 RETURN INVT 37 END INVT 38

8.8 Subroutine STEPVP for the evaluation of end of time step quantities and equilibrium correction terms

With reference to Fig. 8.1, this subroutine evaluates quantities, such as stresses and viscoplastic strains, at the end of the current timestep and also calculates the loading to be applied during the next timestep. The subroutine is structured to perform the following operations sequentially:

(a) All quantities at the end of timestep n are calculated as ( )^{n+1} .
(b) Subroutine INVAR, YIELDF and FLOWVP are called to evaluate the current viscoplastic flow rate, \dot{\epsilon}_{vp}^{n+1} .
(c) The maximum permissible interval length, \Delta t_{n + 1} , for the next timestep as governed by (8.29) and (8.32) is calculated.
(d) The residual forces, \psi^{n+1} , are evaluated and the loads, \Delta V^{n+1} , for the next timestep then calculated.

In the program presented we restrict ourselves to loads applied in discrete increments. An increment of load is applied and the time stepping process is followed until either steady state conditions are achieved, or a specified number of timesteps is reached. Then a further increment of load is applied and the process repeated. Thus in (8.23), \Delta f^{n} = 0 for all stages other than the first timestep of a particular load increment.

The attainment of steady state conditions can be monitored by accumulating some measure of the viscoplastic strain rate for all Gauss points in the structure. At steady state this quantity will become zero. The degree of total viscoplastic flow at any point is best monitored by evaluating the total effective viscoplastic strain rate at all Gauss points according to

\bar {\dot {\epsilon}} _ {v p} = (\sqrt {\frac {2}{3}}) \{(\dot {\epsilon} _ {i j}) _ {v p} (\dot {\epsilon} _ {i j}) _ {v p} \} ^ {1 / 2}. \tag {8.47}

Subroutine STEPVP is now presented and described.

SUBROUTINE STEPVP(ASDIS,COORD,ELOAD,ISTEP,LNODS,LPROP,TIMEX, SPVP 1 . MATNO,MELEM,MMATS,MPOIN,MTOTG,TAUFT,DTIME, SPVP 2 . MTOTV,NDOFN,NELEM,NEVAB,NGAUS,NNODE,NSTR1, SPVP 3 . NTYPE,POSGP,PROPS,NSTRE,NCRIT,STRSG,WEIGP, SPVP 4 . TDISP,VISTN,VIVEL,TLOAD,FTIME,DTINT,IINCS) SPVP 5 C*************** SPVP 6 C C**** EVALUATES QUANTITIES AT END OF TIME STEP AND CALCULATES THE SPVP 7 C RESIDUAL FORCES AND PSEUDO FORCES FOR THE NEXT STEP SPVP 8 C C*************** SPVP 9 C C*************** SPVP 10 C

DIMENSION ASDIS(MTOTV), AVECT(4), CARTD(2,9), COORD(MPOIN,2), SPVP 12
. DEVIA(4), ELCOD(2,9), ELDIS(2,9), ELOAD(MELEM,18), SPVP 13
. LNODS(MELEM,9), POSGP(4), PROPS(MMATS,10), STRAN(4), SPVP 14
. STRES(4), STRSG(4, MTOTG), VIVEL(5, MTOTG), SPVP 15
. VISTN(4, MTOTG), WEIGP(4), DMATX(4,4), TLDIS(2,9), SPVP 16
. DERIV(2,9), SHAPE(9), GPCOD(2,9), TDISP(MTOTV), SPVP 17
. MATNO(MELEM), DJACM(2,2), BMATX(4,18), DESTN(4), SPVP 18
. TLOAD(MELEM,18), SVECT(4) SPVP 19
TWOPI=6.283185308 SPVP 20
DO 10 IELEM=1, NELEM SPVP 21
DO 10 IEVAB=1, NEVAB SPVP 22
10 ELOAD(IELEM,IEVAB)=0.0 SPVP 23
KGAUS=0 SPVP 24
DNEXT=FTIME*DTIME SPVP 25
DO 80 IELEM=1, NELEM SPVP 26
LPROP=MATNO(IELEM) SPVP 27
C
C*** STORE COORDINATES AND INCREMENTAL DISPLACEMENTS OF THE SPVP 29
C ELEMENT NODAL POINTS SPVP 30
C
DO 20 INODE=1, NNODE SPVP 31
LNODE=IABS(LNODS(IELEM,INODE)) SPVP 32
NPOSN=(LNODE-1)*NDOFN SPVP 33
DO 20 IDOFN=1, NDOFN SPVP 34
NPOSN=NPOSN+1 SPVP 35
ELCOD(IDOFN,INODE)=COORD(LNODE,IDOFN) SPVP 36
TLDIS(IDOFN,INODE)=TDISP(NPOSN) SPVP 37
20 ELDIS(IDOFN,INODE)=ASDIS(NPOSN) SPVP 38
THICK=PROPS(LPROP,3) SPVP 39
KGASP=0 SPVP 40
DO 70 IGAUS=1, NGAUS SPVP 41
DO 70 JGAUS=1, NGAUS SPVP 42
EXISP=POSGP(IGAUS) SPVP 43
ETASP=POSGP(JGAUS) SPVP 44
KGAUS=KGAUS+1 SPVP 45
KGASP=KGASP+1 SPVP 46
CALL MODPS(DMATX, LPROP, MMATS, NTYPE, PROPS) SPVP 47
DO 30 ISTR1=1, NSTR1 SPVP 48
30 STRES(ISTR1)=STRSG(ISTR1, KGAUS) SPVP 49
CALL INVAR(DEVIA, LPROP, MMATS, NCRIT, PROPS, SINT3, STEFF, STRES, THETA, SPVP 50
. VARJ2, YIELD) SPVP 51
IF(TIMEX.GT.0.0) CALL TANGVP(LPROP, STRES, PROPS, TIMEX, DTIME, SPVP 52
. NSTRE, NTYPE, MMATS, NCRIT, DMATX) SPVP 53
CALL SFR2(DERIV, ETASP, EXISP, NNODE, SHAPE) SPVP 54
CALL JACOB2(CARTD, DERIV, DJACB, ELCOD, GPCOD, IELEM, KGASP, NNODE, SHAPE) SPVP 55
DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)- SPVP 56
IF(NTYPE.EQ.3) DVOLU=DVOLU*TWOPI*GPCOD(1, KGASP) SPVP 57
IF(THICK.NE.0.0) DVOLU=DVOLU*THICK SPVP 58
CALL STRESS(DMATX, LPROP, NTYPE, PROPS, NDOFN, CARTD, ELDIS, SHAPE, SPVP 59
. GPCOD, NSTRE, VIVEL, DTIME, STRSG, KGASP, MTOTG, MMATS, SPVP 60
. SVECT, NNODE, NSTR1, KGAUS, TLDIS) SPVP 61
DO 60 ISTR1=1, NSTR1 SPVP 62
DESTN(ISTR1)=VIVEL(ISTR1, KGAUS)*DTIME SPVP 63
60 VISTN(ISTR1, KGAUS)=VISTN(ISTR1, KGAUS)+DESTN(ISTR1) SPVP 64
DEBAR=SORT((2.0*(DESTN(1)*DESTN(1)+DESTN(2)*DESTN(2)+DESTN(4)* SPVP 65
. DESTN(4))+DESTN(3)*DESTN(3))/3.0) SPVP 66
DO 65 ISTR1=1, NSTR1 SPVP 67
65 STRES(ISTR1)=STRSG(ISTR1, KGAUS) SPVP 68
VIVEL(5, KGAUS)=VIVEL(5, KGAUS)+DEBAR SPVP 69
CALL INVAR(DEVIA, LPROP, MMATS, NCRIT, PROPS, SINT3, STEFF, STRES, THETA, SPVP 70
. VARJ2, YIELD) SPVP 71
CALL YIELDF(AVECT, DEVIA, LPROP, MMATS, NCRIT, NSTR1, SPVP 72
. PROPS, SINT3, STEFF, THETA, VARJ2) SPVP 73
CALL FLOWVP(AVECT, PROPS, LPROP, STEFF, NSTR1, MTOTG, VIVEL, SPVP 74
. YIELD, KGAUS, MMATS, NCRIT, FNORM, ALLOW) SPVP 75
SPVP 76