FEAP User Forum
FEAP => Programming => Topic started by: carlygao on June 30, 2022, 08:54:23 PM
-
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.
-
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
! 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
-
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
-
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.
-
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
-
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.
-
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
-
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.
-
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
-
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.
-
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
-
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.
-
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
-
The value of beta is material dependent. Look at your data and fit appropriate values to:
yinf + (yo - yinf) exp(-beta ep) + Hiso ep