tarletongroup / crystalplasticity Goto Github PK

Dear authors,

We try to a tensile test including only one element and one grain based on UMAT v1.33. We output STRESS and Sigma in the Umat.f. The displacement is correct, but all values are zero during tension. The material properties were used as the same as the testcase.inp.
Any comments are very appreciated.

Dear Eralp
Thank you very much for sharing this CrystalPlasticity code, it helps me a lot.
In creep.f, the second parameter you provided is expressed as bc = creepparam(5).
Have you ever used creep when using the code?
! Obtain creep parameters
! reference creep rate (1/s)
gammadotc = creepparam(1)
! creep stress multiplier (1/MPa)
bc = creepparam(5)
! activation energy for creep (J/mol)
Qc = creepparam(3)
! reference damage rate (1/s)
gammadotd = creepparam(4)
! damage stress multiplier (1/MPa)
bd = creepparam(5)
! activation energy for damage (J/mol)
Qd = creepparam(6)
!

Dear Erlap,

I am writing to bring to your attention a concern pertaining to the CPFEM Umat functionality. Specifically, I have observed successful operation when employing a single CPU, as evidenced by the command:

abaqus job=test1 user=UMAT.f double inter

However, my attempts to extend its execution to multiple CPUs have resulted in convergence issues (unable to solve even a single increment). The command employed for this purpose is as follows:

abaqus job=test1 user=UMAT.f double inter cpus=4

In light of these circumstances, I am inclined to inquire whether the current version of the subroutine necessitates the utilization of a Mutex mechanism, akin to the previous iteration, to ensure harmonious compatibility with parallel execution.

I kindly request your expert insights into this matter. Your guidance would be greatly appreciated.

Thank you for your attention and assistance.

Sincerely,
Vikram Roy

Dear Eralp,

I have noticed a correction that needs to be made in the documentation. I have attached a screenshot below, highlighting the issue:

Upon reviewing the equations marked in the screenshot, it appears that there might be an error. Specifically, it seems that the term "gamma_t" should be replaced with "gammadot _t".
I kindly request you to review this and if necessary correct it in future versions.

Best Regards,
Vikram Roy

Dear Eralp,.

We try to output the contour of SSD (SDV 27) during tension. But we found that all SSD are zero. In page 26 of the documentation, "In this case, the initial dislocation density (SSD density) shall have a non-zero value in order to end up with non-zero values for the activation volume and hence slip". This means that we should set a non-zero value for rho_0(1:nslip), for example "rho_0(1:nslip) = 1." in usermaterials.f Is it right.
thanks very much if any suggestions.
John WG.

Dear Eralp:

Thanks a lot for sharing this code, and it deepens my understanding to crystal plasticity.

I want to know, is it possible to add element C3D8R in this code?

When I run this code with element C3D8 (element number is 128X128X128=2097152), it requires a very large RAM (about 300GB) to operate, it is also a time-consuming work (46h with 48 cpus, Intel(R) Xeon(R) CPU E5-4640 v3 @ 1.90GHz) and the .odb file is quite large (about 33GB). Attached are my simulation files.

It would be great if the element C3D8R could be supported in this CrystalPlasticity code.

sim_file.zip

Kind regards,

Guofeng

Dear Erlap,

I would like to bring to your attention a correction that needs to be made in the CRSS file. In the file, the variable "sintmat2" is declared with a dimension of nslip x nslip. However, I believe the size of this variable should be nslip x maxnloop, based on its later usage in a nested loop. The outer loop iterates from 1 to nslip, while the inner loop iterates from 1 to maxnloop.

Or Dimension of 1, maxnloop will also be sufficient provided the interation coefficient is same for all slip systems

Here is the relevant code snippet:

! Number of defects
nloop = int(irradiationparam(1))
do is = 1, nslip
do il = 1, 3
Ploop(is) = Ploop(is) +
+ sintmat2(is,il)**2. * loop(il)
end do
end do

Could you please review this and provide a suitable correction or clarification?

Thank you.

Best regards,
Vikram

Dear Eralp

Thank you very much for sharing the CrystalPlasticity codes; it is beneficial.
When I followed your Video Tutorials, everything went well except the C3D4 model.
It went well when I simulated the C3D8 models from DREAM3D and NEPER.
When running the C3D4 model, I did change the element number and element type in the userinputs.f file.
The msg file shows the following error. Would you happen to know where the problem should be?
All the comments and suggestions are a tremendous help to me.

Sincerely,
Ta-Te

   Abaqus 2022                                  Date 10-8-2023   Time 10:31:09
   For use by TOKAI under license from Dassault Systemes or its subsidiary.

                                                                                 
                                                                                 
 STEP    1     INCREMENT     1     STEP TIME    0.00    


                        S T E P       1     S T A T I C   A N A L Y S I S


                                                                                          

     AUTOMATIC TIME CONTROL WITH -
          A SUGGESTED INITIAL TIME INCREMENT OF                1.000E-02
          AND A TOTAL TIME PERIOD OF                            100.    
          THE MINIMUM TIME INCREMENT ALLOWED IS                1.000E-05
          THE MAXIMUM TIME INCREMENT ALLOWED IS                 1.00    

     LINEAR EQUATION SOLVER TYPE         DIRECT SPARSE

 CONVERGENCE TOLERANCE PARAMETERS FOR FORCE    
     CRITERION FOR RESIDUAL FORCE     FOR A NONLINEAR PROBLEM          5.000E-03
     CRITERION FOR DISP.    CORRECTION IN A NONLINEAR PROBLEM          1.000E-02
     INITIAL VALUE OF TIME AVERAGE FORCE                               1.000E-02
     AVERAGE FORCE     IS TIME AVERAGE FORCE    
     ALTERNATE CRIT. FOR RESIDUAL FORCE     FOR A NONLINEAR PROBLEM    2.000E-02
     CRITERION FOR ZERO FORCE     RELATIVE TO TIME AVRG. FORCE         1.000E-05
     CRITERION FOR RESIDUAL FORCE     WHEN THERE IS ZERO FLUX          1.000E-05
     CRITERION FOR DISP.    CORRECTION WHEN THERE IS ZERO FLUX         1.000E-03
     CRITERION FOR RESIDUAL FORCE     FOR A LINEAR INCREMENT           1.000E-08
     FIELD CONVERSION RATIO                                             1.00    
     CRITERION FOR ZERO FORCE     REL. TO TIME AVRG. MAX. FORCE        1.000E-05
     CRITERION FOR ZERO DISP.    RELATIVE TO CHARACTERISTIC LENGTH     1.000E-08

     VOLUMETRIC STRAIN COMPATIBILITY TOLERANCE FOR HYBRID SOLIDS       1.000E-05
     AXIAL STRAIN COMPATIBILITY TOLERANCE FOR HYBRID BEAMS             1.000E-05
     TRANS. SHEAR STRAIN COMPATIBILITY TOLERANCE FOR HYBRID BEAMS      1.000E-05
     SOFT CONTACT CONSTRAINT COMPATIBILITY TOLERANCE FOR P>P0          5.000E-03
     SOFT CONTACT CONSTRAINT COMPATIBILITY TOLERANCE FOR P=0.0         0.100    
     CONTACT FORCE ERROR TOLERANCE FOR CONVERT SDI=YES                 1.00    
     DISPLACEMENT COMPATIBILITY TOLERANCE FOR DCOUP ELEMENTS           1.000E-05
     ROTATION COMPATIBILITY TOLERANCE FOR DCOUP ELEMENTS               1.000E-05

 EQUILIBRIUM WILL BE CHECKED FOR SEVERE DISCONTINUITY ITERATIONS

 TIME INCREMENTATION CONTROL PARAMETERS:
     FIRST EQUILIBRIUM ITERATION FOR CONSECUTIVE DIVERGENCE CHECK              4
     EQUILIBRIUM ITERATION AT WHICH LOG. CONVERGENCE RATE CHECK BEGINS         8
     EQUILIBRIUM ITERATION AFTER WHICH ALTERNATE RESIDUAL IS USED              9
     MAXIMUM EQUILIBRIUM ITERATIONS ALLOWED                                   16
     EQUILIBRIUM ITERATION COUNT FOR CUT-BACK IN NEXT INCREMENT               10
     MAXIMUM EQUILIB. ITERS IN TWO INCREMENTS FOR TIME INCREMENT INCREASE      4
     MAXIMUM ITERATIONS FOR SEVERE DISCONTINUITIES                            50
     MAXIMUM ATTEMPTS ALLOWED IN AN INCREMENT                                  5
     MAXIMUM DISCON. ITERS IN TWO INCREMENTS FOR TIME INCREMENT INCREASE      50
     MAXIMUM CONTACT AUGMENTATIONS FOR *SURFACE BEHAVIOR,AUGMENTED LAGRANGE   50
     CUT-BACK FACTOR AFTER DIVERGENCE                                 0.2500    
     CUT-BACK FACTOR FOR TOO SLOW CONVERGENCE                         0.5000    
     CUT-BACK FACTOR AFTER TOO MANY EQUILIBRIUM ITERATIONS            0.7500    
     CUT-BACK FACTOR AFTER TOO MANY SEVERE DISCONTINUITY ITERATIONS   0.2500    
     CUT-BACK FACTOR AFTER PROBLEMS IN ELEMENT ASSEMBLY               0.2500    
     INCREASE FACTOR AFTER TWO INCREMENTS THAT CONVERGE QUICKLY        1.500    
     MAX. TIME INCREMENT INCREASE FACTOR ALLOWED                       1.500    
     MAX. TIME INCREMENT INCREASE FACTOR ALLOWED (DYNAMICS)            1.250    
     MAX. TIME INCREMENT INCREASE FACTOR ALLOWED (DIFFUSION)           2.000    
     MINIMUM TIME INCREMENT RATIO FOR EXTRAPOLATION TO OCCUR          0.1000    
     MAX. RATIO OF TIME INCREMENT TO STABILITY LIMIT                   1.000    
     FRACTION OF STABILITY LIMIT FOR NEW TIME INCREMENT               0.9500    
     TIME INCREMENT INCREASE FACTOR BEFORE A TIME POINT                1.000    
     GLOBAL STABILIZATION CONTROL IS NOT USED

          PRINT OF INCREMENT NUMBER, TIME, ETC., EVERY    1  INCREMENTS

     THE MAXIMUM NUMBER OF INCREMENTS IN THIS STEP IS                   10000

          LARGE DISPLACEMENT THEORY WILL BE USED

     LINEAR EXTRAPOLATION WILL BE USED

     CHARACTERISTIC ELEMENT LENGTH     0.189    

     DETAILS REGARDING ACTUAL SOLUTION WAVEFRONT REQUESTED

     DETAILED OUTPUT OF DIAGNOSTICS TO DATABASE REQUESTED

     PRINT OF INCREMENT NUMBER, TIME, ETC., TO THE MESSAGE FILE EVERY     1  INCREMENTS

     COLLECTING MODEL CONSTRAINT INFORMATION FOR OVERCONSTRAINT CHECKS

     COLLECTING STEP CONSTRAINT INFORMATION FOR OVERCONSTRAINT CHECKS


  INCREMENT     1 STARTS. ATTEMPT NUMBER  1, TIME INCREMENT  1.000E-02
 
	SYMMETRIC DIRECT SPARSE SOLVER RUNNING ON
 	1 HOST:        1 MPI RANK  x 1 THREAD
 	NUMBER OF EQUATIONS:     2649
 	NUMBER OF RHS:           1
 	NUMBER OF FLOPS:         7.698e+07

     CHECK POINT   START OF SOLVER    

     CHECK POINT  END OF SOLVER       

       ELAPSED USER TIME (SEC)      =   0.0000    
       ELAPSED SYSTEM TIME (SEC)    =  0.10000    
       ELAPSED TOTAL CPU TIME (SEC) =  0.10000    
       ELAPSED WALLCLOCK TIME (SEC) =          0

               CONVERGENCE CHECKS FOR EQUILIBRIUM ITERATION     1


 AVERAGE FORCE                      7.999E-05   TIME AVG. FORCE       7.999E-05
 LARGEST RESIDUAL FORCE             6.741E-09   AT NODE        344   DOF  3
   INSTANCE: NEPER-1                                                                         
 LARGEST INCREMENT OF DISP.         1.000E-06   AT NODE         38   DOF  1
   INSTANCE: NEPER-1                                                                         
 LARGEST CORRECTION TO DISP.        1.000E-06   AT NODE         38   DOF  1
   INSTANCE: NEPER-1                                                                         
          DISP.    CORRECTION TOO LARGE COMPARED TO DISP.    INCREMENT
 
	SYMMETRIC DIRECT SPARSE SOLVER RUNNING ON
 	1 HOST:        1 MPI RANK  x 1 THREAD
 	NUMBER OF EQUATIONS:     2649
 	NUMBER OF RHS:           1
 	NUMBER OF FLOPS:         7.698e+07

     CHECK POINT   START OF SOLVER    

     CHECK POINT  END OF SOLVER       

       ELAPSED USER TIME (SEC)      =   0.0000    
       ELAPSED SYSTEM TIME (SEC)    =   0.0000    
       ELAPSED TOTAL CPU TIME (SEC) =   0.0000    
       ELAPSED WALLCLOCK TIME (SEC) =          0

               CONVERGENCE CHECKS FOR EQUILIBRIUM ITERATION     2


 AVERAGE FORCE                      7.999E-05   TIME AVG. FORCE       7.999E-05
 LARGEST RESIDUAL FORCE             8.755E-13   AT NODE        857   DOF  1
   INSTANCE: NEPER-1                                                                         
 LARGEST INCREMENT OF DISP.         1.000E-06   AT NODE         38   DOF  1
   INSTANCE: NEPER-1                                                                         
 LARGEST CORRECTION TO DISP.        8.147E-13   AT NODE          1   DOF  2
   INSTANCE: NEPER-1                                                                         
          THE FORCE     EQUILIBRIUM EQUATIONS HAVE CONVERGED

---------- RUNTIME EXCEPTION HAS OCCURRED ----------
*** ERROR: ABAQUS/standard rank 0 encountered a SEGMENTATION FAULT

Dear authors,

Thanks very much for your sharing. I remember the old code can consider the cohesive element and the effect of twin, but the new code deleted these functions. Is it ok for these functions in the old code?
Best regards,
Mike

Hi Erlap,

In the .log file, it seems there is a "STATEV_legend.txt" file for our reference. But I couldn't find the file in the Abaqus temporary folder. I ran the simulation via the Job model in Abaqus CAE, could you tell me where can I find the file?

Bests,
Feiyu

Subject: Discrepancies Found Between Documentation and CPSOLVER.for Code

Dear ED Tarleton Group,

I am writing to bring to your attention certain inconsistencies I have observed between the code and the formulae mentioned in the documentation file for CPSOLVER.for.

Specifically, these deviations pertain to the CP_forwardgradientPredictor subroutine. In the documentation, the formula is presented as C:P + WSigma - SigmaW, whereas in the code, it is programmed as C:SchmidVector + SigmaW - WSigma.

Additionally, I have noticed a variation in the calculation of P and SchmidVector. According to the documentation, P is defined as 0.5*(Schmid + Transpose(Schmid)). However, in the SchmidVector variable, the shear terms are doubled, indicated as SchmidVector(4, 5, 6) = (Schmid + Transpose(Schmid)).

I kindly request your assistance in clarifying which version is the accurate one. It would be greatly appreciated if you could spare some time to address this matter as soon as possible.

Thank you for your attention to this matter, and I look forward to your prompt reply.

Best regards,
Vikram Roy
India

Dear Eralp,

Recently，I generated a RVE with 3 grains and set different material property（hardeningparam1）to them，see following figure and “!!!Sum.xlsx” in attached link.

I found that the Mises stress extracted from Grain 3 (element 910) in Job8-1 and Job10-1 have the same value (see following figure and “!!!Sum.xlsx” in attached link), despite the difference of material property.

I don’t know whether it is the issue of the code. The sim files are attached in the following link.

https://doi.org/10.5281/zenodo.8385099

Bests,
Guofeng

Dear Eralp,

In the file globalvariables.f, lines 260-268 define the variables statev_ssd and statev_ssdtot as follows:

! SSD density per slip system at current time step
real(8), allocatable, public :: statev_ssd(:,:,:)

! SSD density per slip system at former time step
real(8), allocatable, public :: statev_ssd_t(:,:,:)

! Total SSD density at current time step
real(8), allocatable, public :: statev_ssdtot(:,:)

! Total SSD density at former time step
real(8), allocatable, public :: statev_ssdtot_t(:,:)

This leads me to believe that statev_ssdtot represents the total ssd density at the current time step. However, I found a discrepancy in the file Initializations.f, lines 573-579, where these variables are initialized:

! Initialize ssd density per slip system
statev_ssd(noel,npt,1:maxnslip) = rho_0
statev_ssd_t(noel,npt,1:maxnslip) = rho_0

! Initialize total ssd density
statev_ssdtot(noel,npt) = rho_0(1)
statev_ssdtot_t(noel,npt) = rho_0(1)

The variable rho_0 is read from props and represents the initial dislocation density for each slip system. It appears that the initialization for statev_ssdtot and statev_ssdtot_t is using only the dislocation density of the first slip system, which might not be correct.

I suggest a correction to the initialization of statev_ssdtot and statev_ssdtot_t to ensure it represents the total SSD density:

! Initialize ssd density per slip system
statev_ssd(noel, npt, 1:maxnslip) = rho_0
statev_ssd_t(noel, npt, 1:maxnslip) = rho_0

! Initialize total ssd density
statev_ssdtot(noel, npt) = sum(rho_0) ! Sum of dislocation density over all slip systems, representing total SSD density
statev_ssdtot_t(noel, npt) = sum(rho_0) ! Sum of dislocation density over all slip systems at former time step

By making this correction, statev_ssdtot will now be properly initialized as the sum of dislocation density over all slip systems, which represents the total SSD density.

Kindly review this suggestion and provide any necessary clarifications or corrections.

Thank you.

Sincerely,
Vikram Roy

Subject: Inquiry about GND Implementation in CPSolver.f

Dear Eralp,

I trust this message finds you well. It's been a while since we last corresponded, back in the months of May-August, discussing various aspects of CPFEM and its implementation. Unfortunately, due to other professional commitments, I had to take a hiatus from working on this topic. However, I have recently resumed my efforts in this area.

During my review of the implementation of GNDs and the subroutine totalandforest in the file cpsolver.f, I noticed a discrepancy in the handling of normal components of edge dislocations compared to the conventional approach outlined in many references. Specifically, your implementation appears to omit the definition of the normal components of edge dislocations.

In the typical approach, the forest component of the GND involves considering contributions from both edge dislocations (in transverse and normal directions) and screw dislocations.

However, in your code snippet, the normal components of edge dislocations are seemingly not defined.

I am reaching out to seek clarification on the rationale behind this implementation choice. Any insights you can provide will be immensely helpful in understanding the nuances of your approach.

Thank you for your time and assistance.

Best regards,
Vikram

Subject: Query Regarding Geometric Factor in CRSS.f Implementation

Dear Eralp,

I hope this message finds you well. I am currently reviewing the implementation of CRSS.f and have a specific question that I would appreciate your assistance in clarifying.

In the code snippet provided, there is a parameter referred to as the "geometric factor," and its value is typically set around 0.25, as indicated in the usermaterials.f file.

Traditionally, the calculation of CRSS from dislocations involves the use of an interaction matrix. I'm curious about why this geometric factor is employed in the current implementation and whether it is somehow related to the Taylor factor.

Additionally, I am interested in understanding how one determines the appropriate values for this geometric factor. Could you provide insights or guidance on this matter?

Your clarification on these points would be immensely helpful in deepening my understanding of the code.

Thank you for your time and consideration.

Best Regards,
Vikram Roy

Hello, excuse me to bother you:
You know that creep is not only related to temperature and stress, but also time is an indispensable factor. In the static analysis of ABAQUS, time is a variable with no practical significance, and no variables related to time are found in creep.f. So how do I use the UMAT subroutine to obtain a result that often appears in the literature!

Hi Eralp,
I hope this message finds you well. I'd like to express my gratitude for sharing the code; it has been incredibly helpful in my work. As I was going through the code and the equations related to hardeningmodel = 4, I came across a discrepancy between the equation presented in the code documentation and the one cited in the reference (https://doi.org/10.1016/j.actamat.2010.06.021). The equation in the documentation and code is:

The equation mentioned in the cited paper is:

Specifically, it appears that "b" is missing in the code documentation's equation when compared to the cited paper. I wanted to confirm whether I might be overlooking something or if there's an error or modification involved in the code.

Furthermore, I found another paper (https://doi.org/10.1016/j.ijplas.2019.09.002) where the authors have used the same hardening model, but in this case, "f" is missing for the same equation. The rest of the model appears to be quite similar to the one cited in the documentation.

My goal is to obtain a stress-strain curve for IN718. I have applied a strain of 1% in the x-direction, and I'm using slipmodel = 3 and hardeningmodel = 4. I have created a new custom material 11 by copying material 2 and commenting out Hintmat1 calculations. I've taken the material parameters from the paper (https://doi.org/10.1016/j.ijplas.2019.09.002). I got this as a response for S11 vs. LE11:

However, the response I'm getting is a straight line. I'm uncertain whether plotting S11 vs. LE11 is an accurate representation of the stress-strain curve. Is it happening due to parameter values or something else?

Thank you for all your help. I hope to hear back from you soon!

Best,
Ajay

I simulated uniaxial tensile according to the video introduction, but an error occurred.
The executable standard.exe aborted with system error code 1073741819. Please check the .dat, .msg, and .sta files for error messages if the files exist. If there are no error messages and you cannot resolve the problem, please run the command "abaqus job=support information=support" to report and save your system information. Use the same command to run Abaqus that you used when the problem occurred. Please contact your local Abaqus support office and send them the input file, the file support.log which you just created, the executable name, and the error code.
Could you help me? Thanks

Dear Eralp,

Thank you for providing clarification on the variable statev_ssdtot in the file Initializations.f. I have another question now, which relates to the sintmat1 and hintmat1 matrices used for dislocations.

I understand that the strength interaction matrix elements, represented by aij, are utilized to calculate the CRSS (Critical Resolved Shear Stress) due to forest hardening, as shown in the following formula:

tau_crss(i) = G12 * burgerv * sqrt(aij * rhoj)

where tau_crss(i) denotes the CRSS of slip system i, G12 is a material constant, burgerv is the magnitude of the Burgers vector, and rhoj corresponds to the dislocation density on the slip system.

While this formula and the values used by Monnet et. al. 2019 for austenitic stainless steel are clear, I am currently facing confusion regarding the hintmat1 matrix. Unfortunately, I couldn't find any details about its use, meaning, or application in the documentation or source code.

Could you please provide insight into the purpose of the hintmat1 matrix and any references where these hardening interaction matrices are employed?

Thank you for your assistance.

Best regards,
Vikram Roy

Hi Eralp,

I have a question regarding the coding of the "gammasum" variable in the latest version of the code.

In the previous versions of the code, it was implemented as follows:
! Total slip over time per slip system
gammasum = 0.
do is = 1, nslip
gammasum(is) = gammasum_t(is) + abs(gammadot(is)) * dt
enddo
However, in the current version of the code, it is written as:

! Total slip over time per slip system
gammasum = 0.
do is = 1, nslip
gammasum(is) = gammasum_t(is) + gammadot(is) * dt
enddo
I have concerns about the logic behind the current implementation. Since slip is an irreversible process, it seems unlikely that positive and negative slip would cancel each other out perfectly at the same location. Could you please clarify these queries and provide an explanation for the change in the code?

Thank you for your assistance.

Dear Eralp
Thank you very much for sharing this CrystalPlasticity code, it helps me a lot.
Q1：
When I use the HCP model you provided, that is, case(3) in usermaterials.f, the program will terminate and the following error message will appear.
INCREMENT 1 STARTS. ATTEMPT NUMBER 1, TIME INCREMENT 1.000E-02
***ERROR: Analysis is being terminated from a user subroutine

Q2：When I use case(2) like the tutorial you provided, everything works as usual,
but when I use the creep model provided in the program, the program will report an error as follows:
CONTACT FORCE ERROR TOLERANCE FOR CONVERT SDI=YES 1.00

Among them, the creep parameters I used are as follows:
! Obtain creep parameters
! reference creep rate (1/s)
creepparam(1) = 10.d-8
! creep stress multiplier (1/MPa)
creepparam(2) = 0.
! activation energy for creep (J/mol)
creepparam(3) = 650.d3
! reference damage rate (1/s)
creepparam(4) = 0.
! damage stress multiplier (1/MPa)
creepparam(5) = 0.
! activation energy for damage (J/mol)
creepparam(6) = 0.
Would you happen to know where the problem should be? All the comments and suggestions are a tremendous help to me.

Sincerely,
Xusheng

Hi Eralp,

I have identified some issues in the CP_Solver.for file that I would like to bring to your attention. Could you please respond to the following queries?

1. In the CP_Solver.for file, lines 801-803 update the elastic strains in the crystal reference system as follows:
   ! Elastic strains in the crystal reference
   dummy33_ = matmul(gmatinv,dummy33)
   dummy33 = matmul(dummy33_,transpose(gmatinv))

! Vectorize
call matvec6(dummy33,dummy6)

! Shear corrections
dummy6(4:6) = 2.0*dummy6(4:6)

Eec = Eec_t + dummy6

However, in the rest of the subroutines (CP_Explicit, CP_Implicit, and CP_ForwardGradientPlastic), the update is performed differently:
! Elastic strains in the crystal reference
dummy33_ = matmul(transpose(gmatinv),dstrane33)
dummy33 = matmul(dummy33_,gmatinv)

! Vectorize
call matvec6(dummy33,dummy6)

! Shear corrections
dummy6(4:6) = 2.0*dummy6(4:6)

Eec = Eec_t + dummy6

I believe the second version is correct. Could you please clarify this?

2. The crystal orientations are updated using two different approaches in different subroutines.
   Approach 1 (lines 813-823 and in the CP_ForwardGradientPlastic subroutine) updates the orientations as follows
   ! Update orientations
   ! All orientation changes are elastic rotations
   dW33=0.
   dW33(1,2) = domega(1)
   dW33(1,3) = -domega(2)
   dW33(2,1) = -domega(1)
   dW33(2,3) = domega(3)
   dW33(3,1) = domega(2)
   dW33(3,2) = -domega(3)

gmatinv = gmatinv_t + matmul(dW33,gmatinv_t)

Approach 2 (in the CP_Explicit and CP_Implicit subroutines) updates the orientations as follows:
! Orientation update
!
! Intermediate variable
dR = I3 - We*dt

! Invert or transpose since rotations are orthogonal
dR = transpose(dR)

! Update the crystal orientations
gmatinv = matmul(dR, gmatinv_t)

Can we follow a uniform approach for updating the crystal orientations in all parts of the code, as suggested below?
gmatinv = matmul(R_e, gmatinv_t)
Here, R_e is the elastic rotation tensor calculated as follows:
F_e = V_e * R_e
Or R_e = inv(V_e) * F_e

V_e is the elastic strench tensor calculated as V_e = sqrt( F_e*transpose(F_e))

And F_e is the elastic deformation gradient calculated as:
F = F_e * F_p
Or F_e = inv(F_p) * F

Since F_p and F are available in all subroutines, can we adopt this uniform approach for updating the crystal orientations throughout the CP_Solver.for file?

I hope to receive an early response.

Thank you.

Best regards,
Vikram Roy

Dear Eralp,

I wanted to bring to your attention a small correction that I believe needs to be made in the subroutine. The correction pertains to lines 1857 and 2915.
In the latest code line no are 1848 and 2906

In both of these lines, the current formula reads as follows:
dstranp33 = 0.5*(Lp - transpose(Lp))*dt

After carefully reviewing the code, I believe that the negative sign in these lines should be replaced with a positive sign, resulting in the corrected formula:
dstranp33 = 0.5*(Lp + transpose(Lp))*dt

I kindly request you to review this correction and provide your comments or suggestions. Your timely response would be greatly appreciated.

Thank you for your attention to this matter.

Best regards,
Vikram Roy

Dear Eralp
Thank you very much for sharing this CrystalPlasticity code, it helps me a lot.
The testcase.inp dosen't have boundary condition, so I add boundary condition according to your video tutorials, but I got an error:

Too many attempts made for this increment.
Abaqus/Standard 2021 DATE 14-8-2023 TIME 11:47:12
SUMMARY OF JOB INFORMATION:
STEP INC ATT SEVERE EQUIL TOTAL TOTAL STEP INC OF DOF IF
DISCON ITERS ITERS TIME/ TIME/LPF TIME/LPF MONITOR RIKS
ITERS FREQ
1 1 1U 0 1 1 0.00 0.00 0.01000
1 1 2U 0 1 1 0.00 0.00 0.005000
1 1 3U 0 1 1 0.00 0.00 0.002500
1 1 4U 0 1 1 0.00 0.00 0.001250
1 1 5U 0 1 1 0.00 0.00 0.0006250

THE ANALYSIS HAS NOT BEEN COMPLETED

Would you happen to know where the problem should be? All the comments and suggestions are a tremendous help to me.

Sincerely,
zz123

Dear Erlap,

Thank you for sharing the code.

In CPFEM, the elastic constants are necessary parameters for simulation, such as C11, C12, and C44 for cubic crystal, or C11, C12, C44, C13, and C33 for hcp. But in the usermaterials.f, you did a re-calculate elastic constants, for example: e1 = (C11**2 + C11*C12 - 2.*C12**2)/(C11 + C12), v12 = C12/(C11 + C12), g12 = C44. Moreover, in hcp material, no single C** is given; only e1, e3, and g12 are found.

What are those re-calculate elastic constants meaning in the code? Why the C** is not needed in the hcp material?

The next question is, if I want to define some fully elastic materials, i.e., α and β pure Titanium, should I only define the C11~C44 regardless of the so-called re-calculate elastic constants, and turn off all the additional models such as the slip and Harding model.

At last, how can I calculate the elastic strain energy through the UMAT?

Thanks for your kind help,
Best,
Feiyu

Hi Eralp,
When I'm printing "Crystal to Sample tranformation: statev_gmatinv" that is data statev_outputs(1),a problem was encountered . The 9 parameters of the rotation matrixwere indeed non-zero, corresponding to the last frame output after load. But the initial rotation matrix of the polycrystalline model with Euler Angle information had 9 zeros. I have the euler Angle and rotation matrix information of the initial polycrystalline model. How can I get the real final Euler Angle information through the SDV 1-9 output by abaqus? Is it through the rotation matrix output by abaqus plus the identity matrix, and then multiplied by the initial Euler Angle rotation matrix, or by other ways?
I would be extremely grateful for any insights, advice, or recommended resources you can offer.

Uploading initial transformation matrix output by abaqus.csv…

Best,
Mengchong Ren

Hi Erlap,

I wanted to share a suggestion regarding the documentation and source code. I believe it will significantly enhance the user-friendliness of the source code.

In the current implementation, certain variables like "KB" and "tolerance" depend on the unit system used for analysis. These variables are stored in different modules, mainly as global variables and user inputs. For a new user, it would be very helpful to include a note in both the documentation and the source code, explaining which variables need to be adjusted when changing units and what preferred units to use with this umat.

For instance, if a user is utilizing stress in MPa units, the present value of tolerance (1.d-8) implemented in user inputs is valid. However, if the user is employing stress in Pa units, adjusting the tolerance values to (1.d-2 or 1.d-1) may be more appropriate.

I became aware of these details while debugging the subroutine in VSCode. By incorporating these clarifications in the documentation and source code, we can minimize sources of error and significantly improve the code's user-friendliness.

Thank you for your attention to this matter.

Best regards,
Vikram Roy

Dear Eralp,

I would like to express my gratitude for your previous assistance regarding the sintmat2 interaction matrix.

I am currently working on the crss module and have some concerns regarding the code for irradiation model-2.

Specifically, I am uncertain about the interpretation of the values r1, r2, ..., rn, and 'a' mentioned in the documentation. In the file usermaterials.f, rn is defined as follows:
! defect vectors (crystallographic directions)
irradiationparam(8:10) = (/ 4.0, 6.0, 11.0 /)
Does this imply that r1 = 4.0, r2 = 6.0, and r3 = 11.0? If so, I am perplexed as to how we can define crystallographic directions using only a single real number.

Additionally, in the file, it is stated that 'a' represents the reaction segment and is defined as (a: reaction segment). Does a = 0 indicate that the screw dislocation entirely absorbs the reaction segment, resulting in helical jogs?
I would greatly appreciate your clarification on these matters.

Sincerely,
Vikram

Dear Erlap,

We intend to combine CP and XFEM to simulate crack propagation. For simplicity, the MAXPS criteria embedded in abaqus are used. When CPFE and XFEM were tested separately, the program ran normally. But when combining both, an error occurs.
abaqus/standard rank 0 encountered a segmentation violation referencing loacation.
---------- RUNTIME EXCEPTION HAS OCCURED ----------
ABAQUS/standard rank 0 terminated by signal 11
*** ERROR CATEGORY: INITIAL STRESS.
Any suggestions are very appreciated.

inp codes: (abaqus 2020)
...
*Elset, elset=_PICKEDSET24, internal, instance=TESS-1, generate
1, 50653, 1
** Constraint:
***Note: Eqn-1 for applying tension displacement
*Equation
2
SET-2, 1, 1.
SET-1, 1, -1.
*Enrichment, name=CRACK-1, type=PROPAGATION CRACK, elset=_PICKEDSET24, interaction=IntProp-1
*End Assembly
**
** MATERIALS
**
** MATERIALS
*Material, name=MATERIAL-GRAIN1
*Damage Initiation, criterion=MAXPS
1000.,
*Damage Evolution, type=DISPLACEMENT
1,
*Damage Stabilization
1e-06
*Depvar
125,
*User Material, constants=6
59.5614, 125.136, -21.5391, 1., 11., 0.
*Material, name=MATERIAL-GRAIN2
*Damage Initiation, criterion=MAXPS
1000.,
*Damage Evolution, type=DISPLACEMENT
1,
*Damage Stabilization
1e-06
*Depvar
125,
*User Material, constants=6
-55.9356, 57.2808, 199.753, 2., 11., 0.
...
** INTERACTION PROPERTIES
**
*Surface Interaction, name=IntProp-1
1.,
*Friction, slip tolerance=0.005
0.115,
*Surface Behavior, pressure-overclosure=LINEAR
100000.,
*Initial Conditions, type=ENRICHMENT
TESS-1.47933, 1, CRACK-1, -0.540541, 0.0777657
TESS-1.47933, 2, CRACK-1, 1e-06, 0.0777657
TESS-1.47933, 3, CRACK-1, 1e-06, 0.170132
TESS-1.47933, 4, CRACK-1, -0.540541, 0.170132
TESS-1.47933, 5, CRACK-1, -0.540541, -0.44067
TESS-1.47933, 6, CRACK-1, 1e-06, -0.44067
TESS-1.47933, 7, CRACK-1, 1e-06, -0.212088
TESS-1.47933, 8, CRACK-1, -0.540541, -0.212088
TESS-1.47970, 1, CRACK-1, -0.540541, 0.170132
TESS-1.47970, 2, CRACK-1, 1e-06, 0.170132
TESS-1.47970, 3, CRACK-1, 1e-06, 0.552352
TESS-1.47970, 4, CRACK-1, -0.540541, 0.552352
TESS-1.47970, 5, CRACK-1, -0.540541, -0.212088
TESS-1.47970, 6, CRACK-1, 1e-06, -0.212088
TESS-1.47970, 7, CRACK-1, 1e-06, 0.170132
TESS-1.47970, 8, CRACK-1, -0.540541, 0.170132
TESS-1.49339, 1, CRACK-1, -0.540541, -0.212088
TESS-1.49339, 2, CRACK-1, 1e-06, -0.212088
TESS-1.49339, 3, CRACK-1, 1e-06, 0.170132
TESS-1.49339, 4, CRACK-1, -0.540541, 0.170132
TESS-1.49339, 5, CRACK-1, -0.540541, -0.44067
TESS-1.49339, 6, CRACK-1, 1e-06, -0.44067
TESS-1.49339, 7, CRACK-1, 1e-06, 0.0777653
TESS-1.49339, 8, CRACK-1, -0.540541, 0.0777653
TESS-1.49302, 1, CRACK-1, -0.540541
TESS-1.49302, 2, CRACK-1, 1e-06
TESS-1.49302, 3, CRACK-1, 1e-06
TESS-1.49302, 4, CRACK-1, -0.540541
TESS-1.49302, 5, CRACK-1, -0.540541
TESS-1.49302, 6, CRACK-1, 1e-06
TESS-1.49302, 7, CRACK-1, 1e-06
TESS-1.49302, 8, CRACK-1, -0.540541
** ----------------------------------------------------------------
**
** STEP: Step-1
**
*Step, name=Step-1, nlgeom=YES, inc=10000
*Static
1., 10., 1e-06, 0.1
**
** BOUNDARY CONDITIONS
**
** Name: Disp-BC-1 Type: Displacement/Rotation
*Boundary
"X=0", 1, 1
** Name: Disp-BC-2 Type: Displacement/Rotation
*Boundary
"Y=0", 2, 2
** Name: Disp-BC-3 Type: Displacement/Rotation
*Boundary
"Z=0", 3, 3
** Name: Disp-BC-4 Type: Displacement/Rotation
*Boundary
SET-1, 1, 1, 2.
**
** INTERACTIONS
**
** Interaction: Int-1
*Enrichment Activation, name=CRACK-1, activate=ON
**
** OUTPUT REQUESTS

Hi Erlap,
Thanks for sharing the code, It really helps me a lot.
I have a problem when I test the creep.f in version 2.4. According to the defination of creep, we should test it under a constant load, therefore I apply a constant pressure on the model during Visco analysis step and run it with case(7) in version 2.4. From the result, it seems like there aren't any creep contribution to the slip deformation. I checked the code and the contribution of creep is indeed considered in the code.I'm confused as to why I don't see the corresponding contribution in the results.
Looking forward to your reply.
Best wishes
Yvfei Nie

Hi Eralp,

I appreciate your response to my previous inquiries. It's great to know that I can assist you in improving your code. I have some further questions regarding the Cp_solver.f file.

1. Please refer to Equation 2.112 in the documentation, which is displayed below:

In this equation, Hab is obtained by dividing delta_Tau_a by delta_gamma_b. However, in lines 3567-3577 of the code, only the diagonal terms of the Hab matrix are calculated as follows:
! Numerical calculation of derivative dtauc/dgamma
Hab = 0.
do is = 1, nslip
!
if (abs(gammadot(is))>sqrt(smallnum)) then
!
Hab(is,is) = (tauceff(is)-tauceff_t(is))/dt/
+ abs(gammadot(is))
!
!
end if
!
end do
I believe an additional loop should be included to calculate the full matrix, or could you please clarify why only the diagonal terms of this matrix are calculated while the non-diagonal terms are omitted?

2. In File UserMaterial.f Line49:
   ! number of slip systems should be corrected as !phase id

I look forward to your prompt response.

Thank you.

Best Regards,
Vikram Roy

Excuse me for my sudden mail.
I really appreciate your belief in sharing codes in crystal plasticity in detail as OSS. I would like to follow your research and make some papers by improving your code on the repository, and I'll share the improved code on it.
Especially because of the above context, I would like to specify this license.
If there are no licenses, I think it's good for you, the first author, to make some OSS licenses (ex., MIT) on the repository.

Best regards,

Hi eralp,
Thanks a lot for your helpful code. I am importing 2D Oxford EBSD Data (.ctf) into Dream3d , to realize the reconstruction and mesh of model, then generate Abaqus input files. Do you know how to set up dream3d to achieve this? I would be extremely grateful for any insights, advice, or recommended resources you can offer.

Best,
Mengchong Ren

Hello Eralp,

I hope you are doing well. I was trying the coupled temp-displacement step module using UMAT_v2.9 for one element case. However, this module uses the C3D8T element type, and the current UMAT version doesn't support this element type. I was interested in simulating behavior under variable temperature and variable strain loading conditions. Is there any way to perform this type of analysis using the current UMAT?
Also, I would appreciate it if you could provide references for CMSX-4 temperature-dependent properties used in the user material file.

Any leads would be very helpful. Thanks!

Dear Eralp,

I hope this message finds you well. I am reaching out to you regarding the use of your codes in my study.

Your program is immensely valuable and essential for my research. As a matter of academic integrity and to ensure proper acknowledgment of your work, I wanted to inquire if you have any specific citation requirements or if there are any particular articles of yours you would prefer me to reference in my paper regarding the use of your program. If there are any specific terms of use or licenses associated with the program, i would be glad to comply with them accordingly. Your guidance or any specific requirements you might have regarding citation in my paper would be immensely helpful for me.

Once again, thank you for your generosity in providing this resource to the academic community. Your assistance has been immensely appreciated and has significantly enriched our study.
I eagerly await your response. If you require any further information or have any queries, please feel free to reach out.

Best regards,
Ren Mengchong

Dear code developer,
We tired a serial of simulations to study the grain size effect (for exampel HP rule). Two cases were simulated firstly. The RVE includes 200 grains and the average grain size is about 30 micrometer or 300 micrometer (These two cases have exactly the same Euler Angle distributions). FCC copper is used (MatID=5 in the code) for verification. In my input and usermaterials.f, Units of stress is MPa. Dimension units are micrometers The RVE generators also generate models with micrometer dimensions.
According to HP rule, the strength of 30 micron grains should be much higher than that of 300 micron grains. But in our simulations, the calculated results do show that the macroscopic stress-strain curves of the two cases basically coincide.
Did we neglect to modify any parameters? Any suggestions are very appreciated.

Best regards,
Mike Jiang

Dear Eralp,

I want to express my gratitude for clarifying the doubts I had regarding the irradiation model - 2.

I am now facing difficultly in understanding the Dislocation Hardening formulation in the UMAT code.

In many research papers and book on crystal plasticity modeling, I have come across the formulation for Dislocation Hardening as follows:
τ = G12 * b * √(αij * ρj)

This formulation is depicted in the attached screenshot from the book "Dislocations, Mesoscale Simulations and Plastic Flow" by Kubin, published by Oxford.

However, I have noticed that the formulation you provided in the documentation differs from what I have previously encountered. The formulation you presented is:
τ = G12 * b * √(γ^2 * ρ_forest)

Could you please assist me in understanding the differences between these two models?

Here are the screenshots:

Thank you once again for your help.

Sincerely,
Vikram

Dear Erlap,I'd like to express my sincerest gratitude for your generosity in sharing this invaluable UMAT code. As someone new to the field of crystal plasticity theory, I've come across a few theoretical questions while delving into the documentation and code. I'm hopeful that you can provide some clarifications:

In "slipmodel2," a constant term referred to as the "reference slip rate" is utilized. However, as I've noticed in the relevant literature, this should be related to the evolution of dislocation density. So, I'm considering modifying it to "ραbv0." However, I'm uncertain about what adjustments need to be made in the slip increment calculation process. Additionally, there's some confusion in the documentation regarding dislocation density. For instance, "hardeningmodel3" and "hardeningmodel4" are both related to Kocks-Mecking. Despite reviewing the references you've provided, I've observed that "hardeningmodel4" calculates only the total substructure density as the total dislocation density, excluding forest dislocations. Whereas in "hardeningmodel3," the total dislocation density is equivalent to the total SSD density, and SSD calculations include forest dislocations (assuming GND OFF). I'm curious as to why forest density isn't included in the total dislocation density calculation for "hardeningmodel4".

Furthermore, I've encountered some challenges while using "usermaterials.f." Although the documentation mentions predefined parameters for "alpha-uranium," I've noticed that the "hardeningmodel" for alpha-uranium in "usermaterials.f" is consistently set to 0. Consequently, I've been manually inputting material properties into "PROPS" to use these models. Unfortunately, this approach has resulted in iterative failures in nearly every job run. I would greatly appreciate it if you could provide an example within "usermaterials.f" to demonstrate how to configure material parameters using "hardeningmodel3" and "hardeningmodel4."

Thank you for your assistance.

Warmest regards,
Soula

When I tried to use the UMAT subroutine(Abaqus2020 + Visual studio 2019 + Intel Fortran 2020 on the Windows), the function subroutine "lapverse" in the file "utils.f" called the external functions "DGETRF" and "DGETRI". And they also called more external functions, perplexingly.

I looked for these external functions and it seems that come from the function library "LAPACK", but I don't know how to call the functions in "LAPACK" when running the UMAT subroutine.

Looking forward to your reply, thanks.

I tried to run the single crystal case from the youtube video https://www.youtube.com/watch?v=T1bCw61qMLw
UMAT.f file compile successfully.
However, the job returns error as follows:

*** ABAQUS/standard rank 0 terminated by signal 6 (Aborted)

*** ERROR CATEGORY: ELEMENT LOOP

I attach here the files for troubleshooting.
Job-1.log

Following is the content of message file:

STEP 1 INCREMENT 1 STEP TIME 0.00

                    S T E P       1     S T A T I C   A N A L Y S I S


                                                                                      

 AUTOMATIC TIME CONTROL WITH -
      A SUGGESTED INITIAL TIME INCREMENT OF                1.000E-02
      AND A TOTAL TIME PERIOD OF                            1.00    
      THE MINIMUM TIME INCREMENT ALLOWED IS                1.000E-05
      THE MAXIMUM TIME INCREMENT ALLOWED IS                 1.00    

 LINEAR EQUATION SOLVER TYPE         DIRECT SPARSE

CONVERGENCE TOLERANCE PARAMETERS FOR FORCE
CRITERION FOR RESIDUAL FORCE FOR A NONLINEAR PROBLEM 5.000E-03
CRITERION FOR DISP. CORRECTION IN A NONLINEAR PROBLEM 1.000E-02
INITIAL VALUE OF TIME AVERAGE FORCE 1.000E-02
AVERAGE FORCE IS TIME AVERAGE FORCE
ALTERNATE CRIT. FOR RESIDUAL FORCE FOR A NONLINEAR PROBLEM 2.000E-02
CRITERION FOR ZERO FORCE RELATIVE TO TIME AVRG. FORCE 1.000E-05
CRITERION FOR RESIDUAL FORCE WHEN THERE IS ZERO FLUX 1.000E-05
CRITERION FOR DISP. CORRECTION WHEN THERE IS ZERO FLUX 1.000E-03
CRITERION FOR RESIDUAL FORCE FOR A LINEAR INCREMENT 1.000E-08
FIELD CONVERSION RATIO 1.00
CRITERION FOR ZERO FORCE REL. TO TIME AVRG. MAX. FORCE 1.000E-05
CRITERION FOR ZERO DISP. RELATIVE TO CHARACTERISTIC LENGTH 1.000E-08

 VOLUMETRIC STRAIN COMPATIBILITY TOLERANCE FOR HYBRID SOLIDS       1.000E-05
 AXIAL STRAIN COMPATIBILITY TOLERANCE FOR HYBRID BEAMS             1.000E-05
 TRANS. SHEAR STRAIN COMPATIBILITY TOLERANCE FOR HYBRID BEAMS      1.000E-05
 SOFT CONTACT CONSTRAINT COMPATIBILITY TOLERANCE FOR P>P0          5.000E-03
 SOFT CONTACT CONSTRAINT COMPATIBILITY TOLERANCE FOR P=0.0         0.100    
 CONTACT FORCE ERROR TOLERANCE FOR CONVERT SDI=YES                 1.00    
 DISPLACEMENT COMPATIBILITY TOLERANCE FOR DCOUP ELEMENTS           1.000E-05
 ROTATION COMPATIBILITY TOLERANCE FOR DCOUP ELEMENTS               1.000E-05

EQUILIBRIUM WILL BE CHECKED FOR SEVERE DISCONTINUITY ITERATIONS

TIME INCREMENTATION CONTROL PARAMETERS:
FIRST EQUILIBRIUM ITERATION FOR CONSECUTIVE DIVERGENCE CHECK 4
EQUILIBRIUM ITERATION AT WHICH LOG. CONVERGENCE RATE CHECK BEGINS 8
EQUILIBRIUM ITERATION AFTER WHICH ALTERNATE RESIDUAL IS USED 9
MAXIMUM EQUILIBRIUM ITERATIONS ALLOWED 16
EQUILIBRIUM ITERATION COUNT FOR CUT-BACK IN NEXT INCREMENT 10
MAXIMUM EQUILIB. ITERS IN TWO INCREMENTS FOR TIME INCREMENT INCREASE 4
MAXIMUM ITERATIONS FOR SEVERE DISCONTINUITIES 50
MAXIMUM ATTEMPTS ALLOWED IN AN INCREMENT 5
MAXIMUM DISCON. ITERS IN TWO INCREMENTS FOR TIME INCREMENT INCREASE 50
MAXIMUM CONTACT AUGMENTATIONS FOR *SURFACE BEHAVIOR,AUGMENTED LAGRANGE 50
CUT-BACK FACTOR AFTER DIVERGENCE 0.2500
CUT-BACK FACTOR FOR TOO SLOW CONVERGENCE 0.5000
CUT-BACK FACTOR AFTER TOO MANY EQUILIBRIUM ITERATIONS 0.7500
CUT-BACK FACTOR AFTER TOO MANY SEVERE DISCONTINUITY ITERATIONS 0.2500
CUT-BACK FACTOR AFTER PROBLEMS IN ELEMENT ASSEMBLY 0.2500
INCREASE FACTOR AFTER TWO INCREMENTS THAT CONVERGE QUICKLY 1.500
MAX. TIME INCREMENT INCREASE FACTOR ALLOWED 1.500
MAX. TIME INCREMENT INCREASE FACTOR ALLOWED (DYNAMICS) 1.250
MAX. TIME INCREMENT INCREASE FACTOR ALLOWED (DIFFUSION) 2.000
MINIMUM TIME INCREMENT RATIO FOR EXTRAPOLATION TO OCCUR 0.1000
MAX. RATIO OF TIME INCREMENT TO STABILITY LIMIT 1.000
FRACTION OF STABILITY LIMIT FOR NEW TIME INCREMENT 0.9500
TIME INCREMENT INCREASE FACTOR BEFORE A TIME POINT 1.000
GLOBAL STABILIZATION CONTROL IS NOT USED

      PRINT OF INCREMENT NUMBER, TIME, ETC., EVERY    1  INCREMENTS

 THE MAXIMUM NUMBER OF INCREMENTS IN THIS STEP IS                     100

      LARGE DISPLACEMENT THEORY WILL BE USED

 LINEAR EXTRAPOLATION WILL BE USED

 CHARACTERISTIC ELEMENT LENGTH     0.500    

 DETAILS REGARDING ACTUAL SOLUTION WAVEFRONT REQUESTED

 DETAILED OUTPUT OF DIAGNOSTICS TO DATABASE REQUESTED

 PRINT OF INCREMENT NUMBER, TIME, ETC., TO THE MESSAGE FILE EVERY     1  INCREMENTS

 COLLECTING MODEL CONSTRAINT INFORMATION FOR OVERCONSTRAINT CHECKS

 COLLECTING STEP CONSTRAINT INFORMATION FOR OVERCONSTRAINT CHECKS

INCREMENT 1 STARTS. ATTEMPT NUMBER 1, TIME INCREMENT 1.000E-02

NUMBER OF EQUATIONS = 81 NUMBER OF RHS = 1
NUMBER OF FLOATING PT. OPERATIONS = 4.83E+04

 CHECK POINT   START OF SOLVER    

 CHECK POINT  END OF SOLVER       

   ELAPSED USER TIME (SEC)      =   0.0000    
   ELAPSED SYSTEM TIME (SEC)    =   0.0000    
   ELAPSED TOTAL CPU TIME (SEC) =   0.0000    
   ELAPSED WALLCLOCK TIME (SEC) =          0

           CONVERGENCE CHECKS FOR EQUILIBRIUM ITERATION     1

AVERAGE FORCE 0.136 TIME AVG. FORCE 0.136
LARGEST RESIDUAL FORCE 4.407E-04 AT NODE 11 DOF 3
INSTANCE: PART-1-1
LARGEST INCREMENT OF DISP. 1.000E-04 AT NODE 1 DOF 1
INSTANCE: PART-1-1
LARGEST CORRECTION TO DISP. 1.000E-04 AT NODE 1 DOF 1
INSTANCE: PART-1-1
DISP. CORRECTION TOO LARGE COMPARED TO DISP. INCREMENT

NUMBER OF EQUATIONS = 81 NUMBER OF RHS = 1
NUMBER OF FLOATING PT. OPERATIONS = 4.83E+04

 CHECK POINT   START OF SOLVER    

 CHECK POINT  END OF SOLVER       

   ELAPSED USER TIME (SEC)      =   0.0000    
   ELAPSED SYSTEM TIME (SEC)    =   0.0000    
   ELAPSED TOTAL CPU TIME (SEC) =   0.0000    
   ELAPSED WALLCLOCK TIME (SEC) =          0

           CONVERGENCE CHECKS FOR EQUILIBRIUM ITERATION     2

AVERAGE FORCE 0.136 TIME AVG. FORCE 0.136
LARGEST RESIDUAL FORCE -1.858E-08 AT NODE 13 DOF 2
INSTANCE: PART-1-1
LARGEST INCREMENT OF DISP. 1.000E-04 AT NODE 1 DOF 1
INSTANCE: PART-1-1
LARGEST CORRECTION TO DISP. 5.996E-09 AT NODE 7 DOF 2
INSTANCE: PART-1-1
THE FORCE EQUILIBRIUM EQUATIONS HAVE CONVERGED

---------- RUNTIME EXCEPTION HAS OCCURED ----------

ABAQUS/standard rank 0 terminated by signal 6 (Aborted)

Hi Erlap,

Could you give me some reference to the backstress model in your code? including both the local and nonlocal backstress models.

Thank you for your help,

Best regards,
Feiyu

Dear Eralp,

I wanted to bring to your attention a pattern I've observed. It seems that when the ratio of rss/crss exceeds 1, the Newton Raphson loop of the Cp_explicit subroutine struggles to achieve convergence, particularly when employing the power law or double exponential slip model.
In light of this, I believe it would be beneficial to implement a condition that triggers a reduction in the time increment whenever this ratio surpasses 1.

I would greatly appreciate your thoughts and comments on this matter.

Thank you.

Best regards,
Vikram

Dear Erlap,
Thank you so much for sharing the code. Can you tell the references of the back stress model in the code. In addition, can the current version calculate the cyclic fatigue load response?

Thanks again.
mike John

Dear Eralp
When I use case(3) in usermaterials.f, the stress curve I extracted is as shown in the figure below. In my opinion, the normal stress curve should not have this period of decline.The boundary conditions are imposed consistent with the case you provided. I would like to know which parameter affects this stress result.

Hi Eralp,

I'm testing the periodic boundary conditions (PBC) by using a single crystal tensile case.

In the first case, I use the power law model, Voce hardening model, GND model (no backstress model). Everything goes well at the elastic and the small plastic period (Figure a); the stress field is uniform. However, the stress gradually becomes nonuniform with the increasing plastic strain. At first, the difference between the maximum and minimum stress is neglectable (<1e-3 MPa in Figure b), while the difference will increase to 1MPa with the increase of plastic strain. Obviously, this result is wrong because the stress field should be uniform under the PBC. I have also tried the Explicit/Implicit method, and the stress field are both nonuniform.

To understand what happened, I ran the second simulation with the same model setting, except the GND model was off. To my surprise, the result I obtained is a uniform stress field.

I'm curious about what happened to the GND model that led to a nonuniform stress field. Based on the model I selected, the GND is not involved in the slip and hardening model, so the calculation of the GND will not influence the slip and hardening process. Could you advise me on obtaining a uniform stress field under PBC with GND?

Bests,
Feiyu

Hi Erlap,

I'm trying to output the total GND with the power law slip model, Voce type hardening model, and the GND model. I set the 6th State-variable output (Total statistically-stored dislocation density) to 1, but it is strange that the SDV I got is all zero even when the GND is non-zero. I didn't use the SSD model because I only care about the GND, so the statev_ssdtot should be the total GND. Is that right? Is there something that I was missing? Or how can I export the total GND?

Bests,
Feiyu