about the codes for the mises plasticity model

[carlygao]

Hi, for the FEAP85, in the file mises.f, on the line 116, xin = sqrt(dot(xi,xi,ntm) + dot(xi(4),xi(4),ntm-3)) to compute the trial norm of stress - back stress. I don't understand the following two questions:

1. xin should be just sqrt(dot(xi,xi,ntm)). What is the second dot doing?
2. For the second dot, xi(4) is a scalar, so it should be dot(xi(4),xi(4),1). For a 3D problem, even though ntm-3 = 3 is used, the dot function only returns xi(4)*xi(4). Why is "ntm-3" used there?

Thanks.

[Prof. S. Govindjee]

 ntm takes values {6,4,1} depending on the space dimension of the problem {3,2,1} and the analysis type (defined by stype); see modlsd.f where ntm is set for small deformation computations.   Here is the stype key from inmate.f

Code: [Select]

!      Geometry types
!         stype     - Geometry identifier
!                     1: Plane stress
!                     2: Plane strain
!                     3: Axisymmetric without torsion
!                     4: Plate/shell
!                     5: Truss
!                     6: Thermal
!                     7: 3-d solid
!                     8: Axisymmetric with torsion
!                     9: Spherical symmetry


[carlygao]

Thanks Prof. Govindjee.

To computer the norm of xin in 3D, I thought it should be sqrt(dot(xi,xi,ntm)), which includes xi(4)*xi(4). Why is xi(4)*xi(4) is added again via the second dot? xi(4) should be the deviatoric stress s12. Is it a special one? Even it needs to be added, I think the second dot should be dot(xi(4),xi(4),1). Why is ntm-3=3 used?

Carly

[Prof. R.L. Taylor]

One is computing  s_ij * s_ij in Voigt notation -- you need to do all 9 term in 3-d less in other dimensions -- look at what Professor Govindjee said the size of ntm would be.

Do what the dot product says to see you the terms are done.

[carlygao]

Thanks Prof. Taylor. I tested the second dot(xi(4),xi(4),3). It returned xi(4)*xi(4)+xi(5)*xi(5)+xi(6)*xi(6), which is what I want.

By the way, I used the example case Iplast_2d in ex8 folder coming with FEAP85, but it gave me the error messages as below.

       T i e  - - -  N o d a l   C o o r d i n a t e s

          The following nodes have been deleted:

      No nodes merged by the tie command

  *ERROR* Residual norm is NaN or Inf
 --> ERRORS OCCURRED: For details see file: Oplast_2d

Lplast_2d:

  SOLUTION SUMMARY
  ----------------

  Load    Total     Solution      Time     Residual Norm        Energy Norm     CPU Time
  Step  Tang+Forms      Time      Incr.  Initial     Final   Initial     Final (Seconds)
  --------------------------------------------------------------------------------------
  *ERROR* Residual norm is NaN or Inf
       1   10   0  1.000E-01  1.00E-01  1.14E+03       NaN  6.21E+15  4.76E+11      0.14

 Total       10        0
  ---------------------------------- END OF FEAP LOG -----------------------------------

Oplast_2d:

     Solution Commands       Variable 1  Variable 2  Variable 3
       NOPRint               0.0000E+00  0.0000E+00  0.0000E+00
       PLOT      RANGe       0.0000E+00  4.4000E+02  0.0000E+00
       PLOT      SYMM        1.0000E+00  1.0000E+00  0.0000E+00
       PLOT      DEFOrm      0.0000E+00  0.0000E+00  0.0000E+00
       LOOP                  4.0000E+01  0.0000E+00  0.0000E+00
       TIME                  0.0000E+00  0.0000E+00  0.0000E+00
       LOOP                  3.0000E+01  0.0000E+00  0.0000E+00
       TANG                  1.0000E+00  0.0000E+00  0.0000E+00
       NEXT                  0.0000E+00  0.0000E+00  0.0000E+00
       PLOT      PSTRe       6.0000E+00  0.0000E+00  1.0000E+00
       NEXT                  0.0000E+00  0.0000E+00  0.0000E+00
 *Command   1 * nopr                v:   0.00       0.00       0.00   
                                                           t=     0.02     0.01
  -> Plots on deformed mesh with scale: 1.00000E+00
                     Eigenvector scale: 0.00000E+00
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of         2 equations.
                 Step:      1 Iteration:     0
 *D4TRI WARNING* Sign of diagonal changed when reducing         2 equations.
                 Step:      1 Iteration:     0
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of         2 equations.
                 Step:      1 Iteration:     1
 *D4TRI WARNING* Sign of diagonal changed when reducing         1 equations.
                 Step:      1 Iteration:     1
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of         2 equations.
                 Step:      1 Iteration:     2
 *D4TRI WARNING* Sign of diagonal changed when reducing        13 equations.
                 Step:      1 Iteration:     2
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of         2 equations.
                 Step:      1 Iteration:     3
 *D4TRI WARNING* Sign of diagonal changed when reducing       112 equations.
                 Step:      1 Iteration:     3
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of         2 equations.
                 Step:      1 Iteration:     4
 *D4TRI WARNING* Sign of diagonal changed when reducing       210 equations.
                 Step:      1 Iteration:     4
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of         1 equations.
                 Step:      1 Iteration:     5
 *D4TRI WARNING* Sign of diagonal changed when reducing       291 equations.
                 Step:      1 Iteration:     5
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of         8 equations.
                 Step:      1 Iteration:     6
 *D4TRI WARNING* Sign of diagonal changed when reducing       296 equations.
                 Step:      1 Iteration:     6
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of        77 equations.
                 Step:      1 Iteration:     7
 *D4TRI WARNING* Sign of diagonal changed when reducing       276 equations.
                 Step:      1 Iteration:     7
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of        81 equations.
                 Step:      1 Iteration:     8
 *D4TRI WARNING* Sign of diagonal changed when reducing       255 equations.
                 Step:      1 Iteration:     8
 *D4TRI WARNING* Lost at least 7 digits in reducing diagonals of        88 equations.
                 Step:      1 Iteration:     9
 *D4TRI WARNING* Sign of diagonal changed when reducing       269 equations.
                 Step:      1 Iteration:     9
  *ERROR* Residual norm is NaN or Inf

[Prof. S. Govindjee]

I just checked Iplast_2d with versions 8.5, 8.6, and the upcoming 8.7 and it runs properly to completion.
There must be a problem with your version of the code -- likely it has been modified or the input file has been changed.
Perhaps you can download a clean version of the code and rebuild.

One thing that I immediately see is wrong is the report that no nodes were merge by the tie.  31 nodes should get merged by the tie command.

[carlygao]

Hi,

I am trying to learn how the plasticity model based J2 works. I am trying to program it by myself in order to understand it deeply. I wrote down the code flow when the first plastic step (when the plasticity model in mises.f is firstly activated) starts. In the elastic stage, the results from my code can match FEAP well.

In the first iteration, it seems I have the right stress (sig), tangent matrix (dd) and lamba (lam) close to the intermediate results from FEAP. However, when the second iteration starts, the strains (eps) is 50 time bigger than the one FEAP provides at the same iteration step, which leads the simulation to diverge soon. I am trying to understand how it happens.

The strains (eps) will be updated for each iteration, and plastic strains will be reset the values at tn whenever the subroutine mises is called. After each iteration, will the plastic strain involve in anything else? Like formulating the right side of the equation? Or just is discard until converged? I don't understand I have the very close eps, very close tangent, stresses and lambda to those from FEAP at the first iteration after calling mises.f, but totally different strains for the second iteration. What did I do differently from FEAP after each iteration. Thanks.

At tn: the last load step for the elastic stage.

At tn+1:

       do niter=1,max_iteration
               a. strains (eps) will be updated after each iteration
               b. plastic strains (epsp, epp(1)) will be reset the values at tn whenever mises is called.
               call mises (3d)
                      1. compute stress (sig)
                      2. compute the tangent matrix
                      3. compute lambda
                      4. compute the plastic strains
       enddo

[Prof. R.L. Taylor]

1. What problem are you trying to solve?
2. Do you have hardening, or doing perfect plasticity (which is sensitive)
3. We always do an elastic solution for the first iteration, often line search is needed.  See how it looks at iteration number, etc.
4. Close is not good -- it should be exactly the same.

[carlygao]

Prof. Taylor,

1. What problem are you trying to solve?

    My input file is attached to this message. It is a steel rod in 3D. 9 with elements and 40 nodes. The force is applied on the 4 nodes on the right surface.

2. Do you have hardening, or doing perfect plasticity (which is sensitive)

    I didn't include any hardening. I thought this is the simplest case to test. If adding hardening will make the model stable, I will do it.

3. We always do an elastic solution for the first iteration, often line search is needed.  See how it looks at iteration number, etc.

    I also let the first iteration to be elastic always, and niter=0. The first iteration I mentioned in the previous message is the first one to activate plasticity, which should be the second iteration in the NR loop.

    I didn't use line search command in the batch. Will it be activated automatically? Since my code doesn't have line search, can I avoid line search in FEAP? I just want to have an apply-to-apply comparison.

4. Close is not good -- it should be exactly the same.

     From the 1st to the 5th step, the simulation is in the elastic stage. For FEAP, the displacement is 0.013321 at the 5th step, and my code gives 0.013359. If you think this difference has problems, I will check whether there is some differences when sig is computed in mises.f. The tangent matrix is exactly the same because it is computed by constants in the elastic stage.

     For FEAP, the index for 4~6 components in strain vector for FEAP is 12, 23, 31, but mine is 23, 13, 12 following the Viogt notation. In the first iteration when the plasticity is activated, the results for the first gaussian point of the first element are as follows:

Input: eps
FEAP: 
1.7999234077834408E-003  -5.3914677266568448E-004  -5.4082850624004493E-004  -2.5399900449410766E-006   9.6277926296760541E-008  -2.3496294179867525E-006

my code:

1.799946179032428E-003   -5.395370954740430E-004    -5.403998800673109E-004     5.420816834977165E-008     -2.563256898508662E-006   -2.655868762724535E-006

Output: before leaving the subroutine mises.f (plasticity is activated at the first time)

The tangent matrix: (dd)

FEAP:

   166697969233.50208        166597279608.37469        166704751158.12329        81159308.060911313       -3076330.9074173048        75076789.431979984     
   166597279608.37469        241517938176.75647        91884782214.868866       -40535906.106768847        1536507.1954880173       -37497925.498925500     
   166704751158.12329        91884782214.868866        241410466627.00790       -40623401.954142474        1539823.7119292878       -37578863.933054492     
   81159308.060911313       -40535906.106768847       -40623401.954142474        74789672984.859543        2504.5342714032417       -61122.290734704759     
  -3076330.9074173048        1536507.1954880173        1539823.7119292878        2504.5342714032417        74789738964.179153        2316.8308923752452     
   75076789.431979999       -37497925.498925507       -37578863.933054492       -61122.290734704766        2316.8308923752456        74789682517.658295

My code: (no damping)

  166697959616.526        166623445927.497        166678582199.493    -1732087.20802763        81902499.5670074        84861681.3693673   
  166623445927.497        241491673088.648        91884866670.8445     865564.618278256       -40928600.7339538       -42407373.3126696   
  166678582199.493        91884866670.8445        241436536815.893    866522.613575266       -40973899.9596712       -42454309.2240206   
  -1732087.20802763        865564.618278256        866522.613575266    74789626734.4098        1423.05974918902        1474.47567099524   
   81902499.5670074       -40928600.7339538       -40973899.9596712    1423.05974918902        74789559474.5035       -69721.2256089384   
   84861681.3693673       -42407373.3126696       -42454309.2240206    1474.47567099524       -69721.2256089384        74789554524.2135

Lamba:

FEAP:      5.2985514125994031E-005
My code: 5.298838355166770E-005

sig: (computed in mises.f)

FEAP:

   353326470.33873779        3449573.4652523547        3198020.6348656267       -189965.19267388814        7200.6009848872463       -175728.17105684505

My code:

   353336660.342127            3398498.49784575           3269443.82242841            4054.20867846687       -191705.026741005       -198631.433499675   

en: (normal direction vector)

FEAP:
  0.81649575272696640   -0.40780775440169792  -0.40868799832526853  -6.6473394985869043E-004   2.5196636639941707E-005  -6.1491518316511464E-004
My code:

  0.816495773510109      -0.408022095682386      -0.408473689059109       1.418665213892019E-005    -6.708220378743532E-004    -6.950592025224932E-004

Since my load is a force pull in x direction, I don't know how it is critical for the values in 4~6 components in strain/stress vector are off. I use the typical data for steel, so the numbers look big. For my code, in the second plasticity iteration, the eps11=1, while it is 0.002 for FEAP, and mine iteration diverge very quickly.

Anyway, if you could recommend some simpler case, with less elements and more stable behaviour, I would greatly appreciate. I have no experience to judge or handle the sensitivity of the plasticity model. Thank you very much!

Carly

[Prof. S. Govindjee]

Let me suggest starting with a 1-element problem with uniform extension (displacement control) in one direction, free to contract laterally, and supported to remove rigid body modes.  The big advantage will be that you know the exact answer by hand.  Thus you can see exactly what your element should be computing.  The numbers should be correct to 14-digits.

[carlygao]

Prof. Govindjee, thanks for your suggestions. I set up a one-element testing model, and used it to fix a few bugs in my code. Now I can get the results matching those from FEAP except for this case. The input file is attached to this message. The left surface can shrink, and all degrees of freedoms on the right surface are constrained. The displacement control is applied on the right surface nodes in x direction up to 0.002 via 10 loading steps. The first 8 steps are in elastic stage, and the plasticity is activated in the last two steps. I setup this case in Abaqus too, but got different results from FEAP.

The batch block is shown below. For the first gaussian point, after applying "tang,,1", the strain vector is (0,0,0,0,0,0); after another applying "tang,,1", the strain vector is (2.0000000000000004E-004  -1.5849364905389041E-005  -1.5849364905389034E-005   1.5849364905389041E-005   9.7408788934244539E-021   5.9150635094611017E-005).

I am confused about the strain vector after applying "tang,,1" in the first time. Should it be (0.0002,0,0,0,0,0) because of the displacement loading in x direction?

batch
  opti
  dt,,1
  LOOP,,10
    TIME
    LOOP,,30
      TANG,,1
      plot,mesh
    NEXT
    plot,pers
    plot,hide
    disp,node,1,1,1
    reac,coor,1,1.0
    PLOT,PSTRe,6,,1
  NEXT
end

[Prof. R.L. Taylor]

Did you look at the nodal displacements.  Why not output the results after every "tang,,1" to watch what is happening.

The program allows you to look at things by adding output commands.  If you try to learn, it is best to look at things iteration by iteration.  Learn to use interactive mode so you can do this more easily.

[carlygao]

Hi,

For the J2 plasticity model, "Plastic mises Y_0 Y_inf beta", is there any suggestion about setting up the value for "beta"? Is there any typical range? I used it as "100" sometimes, but not sure whether it is good or not. Thanks.

Carly

[Prof. S. Govindjee]

The value of beta is material dependent.  Look at your data and fit appropriate values to:

Code: [Select]

yinf + (yo - yinf) exp(-beta ep) + Hiso ep