5. B modeling#
5.1. Characteristics of modeling#
We do three calculations on a model using plane constraints (C_PLAN):
Elastic calculation: we charge up to \(p=10\mathrm{MPa}\);
Elasto-plastic calculation: we load up to \(p=230\mathrm{MPa}\);
Elasto-plastic calculation then discharge: we charge to \(p=230\mathrm{MPa}\) then we discharge to \(p=0\).
5.2. Characteristics of the mesh#
We use the same mesh as modeling A (which includes 315 TRIA6 and 686 nodes).
5.3. Elasto-plastic calculation with STAT_NON_LINE#
The aim is to carry out the elasto-plastic calculation with isotropic work hardening given by a tensile curve such that the uniaxial stress tends to a constant value (\(270\mathrm{MPa}\)).
There is therefore a limit load for this structure for which an upper bound is known \({p}_{\text{lim}}<243\mathit{MPa}\). In this modeling, we only charge up to \(\mathrm{230MPa}\) and we proceed with elastic feedback. The limit load case will be dealt with in the next paragraph.
5.3.1. Preparing the order file#
Start the AsterStudy module.
Then in the left column, click on the Case View tab.
We define the command file for the calculation case (right click on CurrentCase and choose Add Stage).
Note: Add commands via Menu Commands → Show All.
The command file for the calculation case is defined. The command file is very similar to the previous model, below, in bold, we indicate the differences:
Read the mesh in MED format: Command LIRE_MAILLAGE.
Orient the normal of the edge on which the traction load will be applied: Category Mesh/Command MODI_MAILLAGE/ORIE_PEAU_2D by assigning the top group in GROUP_MA. We keep the same mesh name when using reuse.
Define the finite elements used: Command AFFE_MODELE to assign the MECANIQUE phenomenon and 2D plane stress modeling (C_PLAN) to all elements.
Read the traction curve provided in the file forma03b.21: Command LIRE_FONCTION/NOM_PARA = “EPSI”.
To**define the material**: Command DEFI_MATERIAU/ELAS and TRACTION (affect the traction curve).
Assign material to all elements: Command AFFE_MATERIAU.
Affect kinematic boundary conditions and loading: Command AFFE_CHAR_CINE/MECA_IMPO for symmetry on the quarter plate (the left and bottom groups).
Affect the load: Command AFFE_CHAR_ MECA/FORCE_CONTOUR for the force distributed over the top of the plate. The easiest way is to define a unit load (FY=1.0), which will then be multiplied by a ramp function over time.
Create a linear ramp function \(f=t\) to multiply the unit mechanical load: Command DEFI_FONCTION. For example, it varies between \((\mathrm{0.,}0.)\) and \((\mathrm{1000.,}1000.)\)
Create time discretion using commands DEFI_LIST_REEL and DEFI_LIST_INST. For example, we can determine 30 time steps up to 300s to match the maximum load to \(300\mathit{MPa}\). In DEFI_LIST_INST, activate the automatic cutting of the time step: ECHEC/EVENEMENT =” ERREUR “and ACTION =” DECOUPE”.
Calculate the evolution elasto-plastic: Command STAT_NON_LINE. We put COMPORTEMENT/RELATION = “VMIS_ISOT_TRAC”, the list of times defined previously in INCREMENT, the materials in CHAM_MATER, MODELE and also the boundary conditions and loading (CHARGE + FONC_MULT) in EXCIT.
5.3.2. Elastic calculation#
If we indicate INST_FIN =10s in the INCREMENT keyword of the STAT_NIN_LINE command, we will have applied a force of \(10\mathit{MPa}\), which is strictly equivalent to the elastic case.
Here are a few things to check:
Verify that we find the same results as in the previous modeling: displacements and constraints.
Check the plasticity indicator (VARI_ELGA). It must be zero everywhere, same for the cumulative plastic deformation (EPSP_ELGA).
Observe the convergence table: number of iterations of Newton.
Vary the temporal discretization (number of time steps).
Compare constraint \({\sigma }_{\theta \theta }\) on the edge of the hole and compare with the analytical solution.
5.3.3. Elasto-plastic calculation under load#
If we specify INST_FIN =230in the INCREMENT keyword of the STAT_NIN_LINE command, we will have applied a force of \(230\mathit{MPa}\) * . *
We must have plasticization of the structure because the resulting stress in a part of the structure is greater than the elastic limit (which is equal to \(\mathrm{200MPa}\)).
With the discretization in 30 steps of the loading time and the activation of the cutting of the time step, there will be no convergence. Newton’s algorithm is failing. To try to converge, you can play with several parameters:
Increase the discretization of the load (be careful not to go too far, less than 100 steps) so that the calculation is not too long!).
Increase Newton’s maximum number of iterations to 10 by default (STAT_NON_LINE/CONVERGENCE/ITER_GLOB_MAXI).
Activate linear search (STAT_NON_LINE/RECH_LINEAIRE).
Increase the number of possible subdivisions of the time step in managing non-convergence (DEFI_LIST_INST/ECHEC/SUBD_NIVEAU).
The combination of all the previous techniques.
Here is a combination that works in our case:
Discretization of loading in 50 steps;
STAT_NON_LINE/CONVERGENCE/ITER_GLOB_MAXI = 20;
Linear search activated;
Cut up to five sub-levels (DEFI_LIST_INST/ECHEC/SUBD_NIVEAU =5).
With this combination, convergence is obtained in 358 steps (instead of the initial 30, because of the cuts) and more than 4300 iterations.
Some interesting results to observe:
At the final moment of the calculation, for the maximum load, we can notice on the cumulative plastic deformation isovalues, the location of the plastic deformations (internal variable V1) in the vicinity of \(B\). Visualization at Gauss points can be used to visualize the cumulative equivalent plastic deformation at the locations where it is calculated.
For a loading under \(\mathrm{66,7}\mathrm{MPa}\), there is no lamination.
Until \(230\mathit{MPa}\), we are constantly in charge.
The maximum value of the Von Mises criterion at Gauss points is always less than or equal to \(270\mathrm{Mpa}\), which shows that the solution does indeed verify the law of behavior.
Let’s look at the convergence table at any step:
Calculation time: 2.297250000000e+02 - Cutting level: 2
We are at moment \(229.725\), we have split the initial list twice and we are converging into 8 Newton iterations. We notice:
The linear search was not very expensive: one iteration only (by default STAT_NON_LINE/RECH_LINEAIRE/ITER_LINE_MAXI =3). We also see that there is no linear search in prediction.
Convergence is tested on criterion RESI_GLOB_RELA. Without changing the default values of the STAT_NON_LINE/CONVERGENCE command, you must therefore have RESI_GLOB_RELA less than \(1.0\mathrm{\times }{10}^{\text{-}6}\) to achieve convergence.
The tangent matrix is only calculated as a prediction, over several time steps, we see that it is always calculated at the first iteration. This corresponds to the default setting for the STAT_NON_LINE/NEWTON command: REAC_INCR =1 and REAC_ITER = 0.
You can ask the STAT_NON_LINEd command to show more information: follow the value of a degree of freedom (keyword SUIVI_DDL), or ask to show where the convergence criterion is the worst (keyword AFFICHAGE/INFO_RESIDU =” OUI “). This last adjustment makes it possible, for example, to know which location in the structure controls convergence.
At the end of the transition statistical information is displayed:
Statistics on the whole transition.
Number of time steps: 358
Number of Newton iterations: 4353
Number of factorizations in the matrix: 358
Number of behavior integrations: 8715
Number of K.U=F resolutions: 4353
Number of linear search iterations: 4004
Time CPU consumed in the transition: 2 m 29 s
including time « lost » in cuts: 55.640 s -> the instant list is 62.8% effective
Matrix assembly time: 0.790 s
Construction time second member: 5.690 s
Total matrix factorization time: 1.390 s
Total integration behavior time: 2 m 4 s
Total time K.U=F resolution: 3.030 s
Time other operations: 14.120 s
We are not going to detail all the information but make a few remarks:
For this we can see that the most consuming item is the integration of the law of behavior, well ahead of factorization and system resolution. This is often the case in 2D, but it is mainly due to the fact that we use a non-complete version of Newton. The matrix is only factored once per time step (we can also see it on the number of factorizations: 358, like the number of time steps).
The initial time list divided into 50 time steps was not the most efficient: the time lost in the calculation due to convergence failures and therefore the redistribution of the time step is about a third of the total time.
To substantially improve convergence, all you need to do is activate the full Newton: STAT_NON_LINE/NEWTON/REAC_ITER =1. Over 50 steps of time, the following results are obtained:
Statistics on the whole transition.
Number of time steps: 50
Number of Newton iterations: 152
Number of factorizations in the matrix: 152
Number of behavior integrations: 202
Number of K.U=F resolutions: 152
Time CPU consumed in the transition: 5.330 s
Matrix assembly time: 0.220 s
Second member construction time: 0.520 s
Total matrix factorization time: 0.650 s
Total integration behavior time: 3.280 s
Total time resolution K.U=F: 0.110 s
Time other operations: 0.550 s
In addition to the increase in the convergence speed (30 times faster), there is no reduction in the time step. The list of moments can be divided even more roughly.
5.3.4. Elasto-plastic calculation under charge and then discharge#
We are now going to proceed with the discharge. For this, a new ramp in the shape of a hat is defined:
\(F\mathrm{=}0.\) for \(t\mathrm{=}0\).
\(F\mathrm{=}230.\) for \(t\mathrm{=}230.\).
\(F\mathrm{=}0.\) for \(t\mathrm{=}300.\).
Then you have to define a new list of related moments: Command DEFI_LIST_REEL (30 time steps until \(t\mathrm{=}230.\), then 10 time steps until \(t\mathrm{=}300.\)) and the command DEFI_LIST_INST by activating the automatic cutting of the time step with ECHEC/EVENEMENT =” ERREUR “and ACTION =” DECOUPE”).
Do a new calculation (a new STAT_NON_LINE command) by taking INST_FIN =300. in the INCREMENT keyword of the STAT_NON_LINE command, and using the new ramp and the new instant list.
If we resume the strategy for carrying out the calculation to the end in plasticity (i.e. up to \(p\mathrm{=}230\mathit{MPa}\) using STAT_NON_LINE/NEWTON/REAC_ITER =1), by minimizing the time division (previous exercise), we see that this strategy must be improved on the plastic dump part because it does not converge with these settings.
It should be noted that the discharge occurs elastically and creates an inelastic deformation when the load is zero. For the calculation to converge, it is therefore necessary to activate the elastic matrix in prediction (STAT_NON_LINE/NEWTON/PREDICTION =” ELASTIQUE “).
The aim is to carry out the elasto-plastic calculation with isotropic work hardening given by a tensile curve such that the uniaxial stress tends to a constant value (\(270\mathrm{MPa}\)).
There is therefore a limit load for this structure for which an upper bound is known: \({p}_{\text{lim}}<243\mathit{MPa}\). In this modeling, we will show how to carry out the calculation in clear of this limit load thanks to piloting.
5.4. Tested sizes and results#
We test the value of the stress components for the elastic calculation (loading of \(\mathrm{10MPa}\)), we must find the same thing as in the A modeling, i.e.:
Component |
Reference type |
Value |
Tolerance |
SIGM_NOEU — \(\mathit{SIYY}\) in \(B\) |
AUTRE_ASTER |
Identical to A |
|
SIGM_NOEU — \(\mathit{SIXX}\) in \(A\) |
AUTRE_ASTER |
Identical to A |
|
We test the value of the stress components and the internal variables for the elasto-plastic calculation (at the moment corresponds to the load of \(\mathrm{230MPa}\)):
Component |
Reference Type |
Tolerance |
SIGM_NOEU — \(\mathit{SIYY}\) in \(B\) |
|
1,00E -006% |
SIGM_NOEU — \(\mathit{SIXX}\) in \(A\) |
|
1,00E -006% |
VARI_NOEU — \(\mathit{V1}\) in \(B\) |
|
1,00E -006% |
VARI_NOEU — \(\mathit{V2}\) in \(B\) |
|
1,00E -006% |
VARI_NOEU — \(\mathit{V1}\) in \(A\) |
|
1,00E -006% |
VARI_NOEU — \(\mathit{V2}\) in \(A\) |
|
1,00E -006% |
We test the value of the stress components and the internal variables for the elasto-plastic calculation with discharge (at the moment corresponds to the final unloading of \(0\)):
Component |
Reference Type |
Tolerance |
SIGM_NOEU — \(\mathit{SIYY}\) in \(B\) |
|
1,00E -006% |
SIGM_NOEU — \(\mathit{SIXX}\) in \(A\) |
|
1,00E -006% |
VARI_NOEU — \(\mathit{V1}\) in \(B\) |
|
1,00E -006% |
VARI_NOEU — \(\mathit{V2}\) in \(B\) |
|
1,00E -006% |
VARI_NOEU — \(\mathit{V1}\) in \(A\) |
|
1,00E -006% |
VARI_NOEU — \(\mathit{V2}\) in \(A\) |
|
1,00E -006% |