PDF

Parameter Estimation of Elastoplastic Materials
Introduction
This tutorial model demonstrates how to estimate the material parameters of an elastoplastic constitutive model based on cyclic shear data. The mathematical definition of the inverse problem is formulated in Parameter Estimation of Hyperelastic Materials, wherein a similar inverse problem is solved for the case of hyperelasticity.
In this example, we assume that the material is isotropic and that the Young’s modulus E and Poisson’s ratio ν are known a priori (for example, extracted directly from the initial part of the stress–strain curve), and we will focus on how to estimate the plastic material behavior. The problem setup is inspired by Ref. 1.
Model Definition
The load case considered is a cyclic shear test, consisting of strain-controlled cyclic loading up to a maximum engineering shear strain γ = 2.5%. The stress–strain curve is shown in Figure 1. The elastic material parameters E = 200 GPa and ν = 0.3 are considered known.
Figure 1: Cyclic simple shear data.
Before setting up the model, it is a good idea to conduct a preliminary analysis of the stress–strain curve to choose an appropriate elastoplastic material model. First, note that initial yielding starts at a shear stress τ of about 200 MPa. We can thus estimate the initial yield stress σy0 as
(1)
by using the definition of the von Mises stress ; s being the deviatoric stress.
Next, we analyze the hardening behavior. Note that the maximum stress in the first loading is about 300 MPa, and the onset of plastic flow during unloading starts around τ = 150 MPa. This indicates a significant Bauschinger effect, which means that the material model needs to include both isotropic and kinematic hardening. Assuming von Mises-associated plasticity, the yield function F can be expressed as
(2)
where sb denotes the back stress, σy is the isotropic hardening function, and εpe is the equivalent plastic strain. For the isotropic hardening part, we select a nonlinear Voce model
(3)
Herein, the two material parameters to estimate are the saturation stress σsat and the saturation exponent β. In particular, this model can capture both cyclic hardening and softening depending on the sign of σsat.
To model kinematic hardening, an evolution equation for the back stress is needed. The simplest form is a linear relation where the back stress is proportional to the plastic strain tensor,
(4)
Here, Ck is the kinematic hardening modulus. In COMSOL Multiphysics, the user input controlling the kinematic hardening is the tangent modulus Ek, which is related to Ck according to
(5)
Since the initial yield stress is given by Equation 1, the inverse problem needs to be solved for three parameters: σsat, β, and Ck.
The linear kinematic hardening model will be compared with the Armstrong–Frederick kinematic hardening model, for which the back stress is computed by integrating the evolution equation
(6)
The second term, which depends on the dimensionless kinematic hardening parameter γk, models the rate of decrease of the hardening modulus upon accumulation of plastic strain, which leads to a stress–strain response that stabilizes upon repeated cyclic loading.
The material parameters for both models along with an initial guess of their values are provided in Table 1.
Results and Discussion
The prediction for the linear kinematic hardening model with the initial parameter values in Table 1 is shown in Figure 2. Although this initial prediction is rather poor, in particular with regard to the cyclic hardening, it is suitable as a starting point for the parameter estimation study.
After having estimated the parameters σsat, β, and Ck of the linear kinematic hardening model, the results with the calibrated material parameters are shown in Figure 3. The final values of the parameters are reported in Table 2.
When switching to the Armstrong–Frederick kinematic hardening model, all four parameters σsat, β, Ck, and γk are estimated. This leads to a significant improvement in the overall prediction, see Figure 4.
Figure 2: Model prediction with the initial guess of the material parameters.
Figure 3: Model prediction after estimating the parameters of the linear kinematic hardening model.
Figure 4: Model prediction after estimating the parameters of the Armstrong–Frederick kinematic hardening model.
Notes About the COMSOL Implementation
The Parameter Estimation functionality is available in COMSOL Multiphysics in the context menu of a Component or under Optimization in the Physics toolbar. To solve an inverse problem, add one Least-Squares Objective subnode for each experiment together with a study containing a Parameter Estimation study step. For most least-squares problems with nonlinear material models, the Levenberg–Marquardt algorithm with a finite difference approximation of the Jacobian is a robust and efficient choice of optimization solver.
References
1. J. Fu, F. Barlat, J.-H. Kim, and F. Pierron, “Identification of nonlinear kinematic hardening constitutive model parameters using the virtual fields method for advanced high strength steels,” Int. J. Solids Struct., vol. 102–103, pp. 30–43, 2016.
Application Library path: Nonlinear_Structural_Materials_Module/Plasticity/parameter_estimation_plasticity
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Structural Mechanics > Solid Mechanics (solid).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies > Stationary.
6
Root
1
In the Model Builder window, click the root node.
2
In the root node’s Settings window, locate the Unit System section.
3
From the Unit system list, choose MPa. The MPa base unit system is often convenient to use when working with structural mechanics problems.
Results
Start by importing the cyclic shear data.
Cyclic Shear Data
1
In the Model Builder window, expand the Results node.
2
Right-click Results > Tables and choose Table.
3
In the Settings window for Table, type Cyclic Shear Data in the Label text field.
4
Locate the Data section. Click  Import.
5
Cyclic Shear Data
1
Go to the Cyclic Shear Data window.
2
Click the Table Graph button in the window toolbar.
Results
Table Graph 1
1
In the Settings window for Table Graph, locate the Data section.
2
From the x-axis data list, choose Engineering shear strain (-).
3
From the Plot columns list, choose Manual.
4
In the Columns list box, select Shear stress (MPa).
5
Locate the Coloring and Style section. Find the Line markers subsection. From the Marker list, choose Point.
Cyclic Shear Data
1
In the Model Builder window, under Results click 1D Plot Group 1.
2
In the Settings window for 1D Plot Group, type Cyclic Shear Data in the Label text field.
3
Click to expand the Title section. From the Title type list, choose Label.
4
Locate the Plot Settings section.
5
Select the x-axis label checkbox. In the associated text field, type Engineering shear strain (1).
6
Select the y-axis label checkbox. In the associated text field, type Shear stress (MPa).
7
In the Cyclic Shear Data toolbar, click  Plot.
Global Definitions
Now, set up the forward model. Start by creating an interpolation function for the applied shear strain as a function of time.
Shear Strain
1
In the Home toolbar, click  Functions and choose Global > Interpolation.
2
In the Settings window for Interpolation, type Shear Strain in the Label text field.
3
Locate the Definition section. From the Data source list, choose Result table.
4
Locate the Data Column Settings section. In the table, click to select the cell at row number 1 and column number 1.
5
In the Unit text field, type s.
6
7
In the Name text field, type shear_strain.
8
In the Unit text field, type 1.
9
Parameters 1
1
In the Model Builder window, click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Geometry 1
Block 1 (blk1)
1
In the Geometry toolbar, click  Block.
2
In the Settings window for Block, locate the Size and Shape section.
3
In the Width text field, type L.
4
In the Depth text field, type L.
5
In the Height text field, type L.
Materials
Material 1 (mat1)
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
In the Settings window for Material, locate the Material Contents section.
3
4
Click to expand the Material Properties section. In the Material properties tree, select Solid Mechanics > Elastoplastic Material > Elastoplastic Material Model.
5
Click  Add to Material.
6
Locate the Material Contents section. In the table, enter the following settings:
7
Locate the Material Properties section. In the Material properties tree, select Solid Mechanics > Elastoplastic Material > Voce.
8
Click  Add to Material.
9
Locate the Material Contents section. In the table, enter the following settings:
10
Locate the Material Properties section. In the Material properties tree, select Solid Mechanics > Elastoplastic Material > Armstrong–Frederick.
11
Click  Add to Material.
12
Locate the Material Contents section. In the table, enter the following settings:
Solid Mechanics (solid)
When solving inverse problems, the forward model will be solved multiple times for each iteration of the optimization solver. If the material tests are designed such that the distributions of stress and strain are homogeneous, it is computationally advantageous to set up a single element model of the experiment by using idealized boundary conditions, linear shape functions, and reduced integration.
1
In the Model Builder window, expand the Material 1 (mat1) node, then click Component 1 (comp1) > Solid Mechanics (solid).
2
In the Settings window for Solid Mechanics, locate the Structural Transient Behavior section.
3
4
Click to expand the Discretization section. From the Displacement field list, choose Linear.
Linear Elastic Material 1
1
In the Model Builder window, under Component 1 (comp1) > Solid Mechanics (solid) click Linear Elastic Material 1.
2
In the Settings window for Linear Elastic Material, locate the Geometric Nonlinearity section.
3
From the Formulation list, choose Geometrically linear.
4
Locate the Quadrature Settings section. Select the Reduced integration checkbox.
Linear Kinematic Hardening
1
In the Physics toolbar, click  Attributes and choose Plasticity.
First, we define an elastoplastic model with Voce isotropic hardening and linear kinematic hardening.
2
In the Settings window for Plasticity, type Linear Kinematic Hardening in the Label text field.
3
Locate the Plasticity Model section. Find the Isotropic hardening model subsection. From the list, choose Voce.
4
Find the Kinematic hardening model subsection. From the list, choose Linear.
Fixed Constraint 1
Define the idealized boundary conditions for simple shear in the xz-plane. Since we use linear shape functions, we do not need to apply any constraints to the x boundaries for them to remain flat during deformation.
1
In the Physics toolbar, click  Boundaries and choose Fixed Constraint.
2
Roller 1
1
In the Physics toolbar, click  Boundaries and choose Roller.
2
Prescribed Displacement 1
1
In the Physics toolbar, click  Boundaries and choose Prescribed Displacement.
2
3
In the Settings window for Prescribed Displacement, locate the Prescribed Displacement section.
4
From the Displacement in x direction list, choose Prescribed.
5
In the u0x text field, type shear_strain(t)*L.
6
From the Displacement in y direction list, choose Prescribed.
7
From the Displacement in z direction list, choose Prescribed.
Definitions
Add a variable for the volume-averaged shear stress. This global variable will be used later on in the Least-Squares Objective.
Average 1 (aveop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Average.
2
Variables 1
1
In the Model Builder window, right-click Definitions and choose Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Mesh 1
Mesh the unit cube with a single hexahedral element.
Mapped 1
1
In the Mesh toolbar, click  More Generators and choose Mapped.
2
Distribution 1
1
Right-click Mapped 1 and choose Distribution.
2
In the Settings window for Distribution, locate the Edge Selection section.
3
From the Selection list, choose All edges.
4
Locate the Distribution section. In the Number of elements text field, type 1.
Swept 1
In the Mesh toolbar, click  Swept.
Distribution 1
1
Right-click Swept 1 and choose Distribution.
2
In the Settings window for Distribution, locate the Distribution section.
3
In the Number of elements text field, type 1.
4
In the Model Builder window, right-click Mesh 1 and choose Build All.
Forward Problem
Solve the forward model once to check that everything is set up correctly.
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, locate the Study Settings section.
3
Clear the Generate default plots checkbox.
4
In the Label text field, type Forward Problem.
Step 1: Stationary
1
In the Model Builder window, under Forward Problem click Step 1: Stationary.
2
In the Settings window for Stationary, click to expand the Study Extensions section.
3
Select the Auxiliary sweep checkbox.
4
5
6
In the Study toolbar, click  Compute.
Results
Compare the initial model prediction with the shear data.
Initial Model Prediction
1
In the Model Builder window, right-click Cyclic Shear Data and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Initial Model Prediction in the Label text field.
3
Locate the Data section. From the Dataset list, choose Forward Problem/Solution 1 (sol1).
4
Locate the Legend section. From the Position list, choose Lower right.
Table Graph 1
1
In the Model Builder window, expand the Initial Model Prediction node, then click Table Graph 1.
2
In the Settings window for Table Graph, click to expand the Legends section.
3
Select the Show legends checkbox.
4
From the Legends list, choose Manual.
5
Global 1
1
In the Model Builder window, right-click Initial Model Prediction and choose Global.
2
In the Settings window for Global, locate the y-Axis Data section.
3
4
Locate the x-Axis Data section. From the Parameter list, choose Expression.
5
In the Expression text field, type shear_strain(t).
6
Click to expand the Legends section. Find the Include subsection. Clear the Solution checkbox.
7
In the Initial Model Prediction toolbar, click  Plot.
You can compute a scalar metric of the quality of the model fit with respect to the data by adding a Comparison subnode. In this example, the root-mean-square (RMS) value is chosen for the comparison.
Comparison 1
1
Right-click Global 1 and choose Comparison.
2
In the Settings window for Comparison, locate the Comparison section.
3
From the Metric list, choose RMS.
4
In the Initial Model Prediction toolbar, click  Plot. The RMS for the initial model prediction should be about 91.5.
Component 1 (comp1)
Now, set up the inverse problem to improve the model parameters.
Least-Squares Objective 1
1
In the Physics toolbar, click  Optimization and choose Parameter Estimation.
2
In the Settings window for Least-Squares Objective, locate the Experimental Data section.
3
From the Data source list, choose Result table.
The first column contains the time. Since the plasticity model is rate independent, we can solve it either using the Time-Dependent or the Stationary continuation solver. Here, we will use the latter. Therefore, set the column type to Parameter.
4
Locate the Data Column Settings section. In the table, enter the following settings:
5
From the Name list, choose t (Time parameter).
6
In the Unit text field, type s.
The second column containing the shear strain data can be ignored in the optimization problem, since we use these data to prescribe the boundary conditions.
7
The third column contains the shear stress data, which will be used in the evaluation of the objective function. In the Model expression field, enter the expression for the global variable of the shear stress τ.
8
9
In the Model expression text field, type comp1.tau.
10
In the Column name text field, type shear_stress.
11
In the Unit text field, type MPa.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies > Stationary.
4
Click the Add Study button in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Parameter Estimation: Linear Kinematic Hardening
1
In the Settings window for Study, type Parameter Estimation: Linear Kinematic Hardening in the Label text field.
2
Locate the Study Settings section. Clear the Generate default plots checkbox.
Parameter Estimation
1
In the Study toolbar, click  Optimization and choose Parameter Estimation.
2
In the Settings window for Parameter Estimation, locate the Experimental Data section.
3
From the Data source list, choose All Least-Squares objectives.
4
Locate the Estimated Parameters section. Click  Add three times.
5
6
Locate the Parameter Estimation Method section. From the Method list, choose Levenberg–Marquardt.
7
From the Least-squares time/parameter list method list, choose Use only least-squares data points.
Step 1: Stationary
Activate the continuation solver for the time parameter t.
1
In the Model Builder window, click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Study Extensions section.
3
Select the Auxiliary sweep checkbox.
4
From the Least-squares continuation parameter list, choose t (Time parameter).
Before computing the study, set up a plot that can be used while solving to monitor the comparison between the data and the model prediction. Start by generating the solution dataset by getting the initial values.
5
In the Study toolbar, click  Get Initial Value.
Results
Parameter Estimation: Linear Kinematic Hardening
1
In the Model Builder window, right-click Initial Model Prediction and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Parameter Estimation: Linear Kinematic Hardening in the Label text field.
3
Locate the Data section. From the Dataset list, choose Parameter Estimation: Linear Kinematic Hardening/Solution 2 (sol2).
Parameter Estimation: Linear Kinematic Hardening
Parameter Estimation
1
In the Model Builder window, under Parameter Estimation: Linear Kinematic Hardening click Parameter Estimation.
2
In the Settings window for Parameter Estimation, click to expand the Output section.
3
Select the Plot checkbox.
4
5
In the Study toolbar, click  Compute.
The RMS is now down to around 19.4 and the model is improved, but the shape and cyclic evolution of the hardening are not well captured. Next, see if the prediction can be improved with a nonlinear kinematic hardening model.
Solid Mechanics (solid)
Armstrong-Frederick Kinematic Hardening
1
In the Model Builder window, right-click Linear Kinematic Hardening and choose Duplicate.
2
In the Settings window for Plasticity, type Armstrong-Frederick Kinematic Hardening in the Label text field.
3
Locate the Plasticity Model section. Find the Kinematic hardening model subsection. From the list, choose Armstrong–Frederick.
Add Study
1
In the Study toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies > Stationary.
4
Click the Add Study button in the window toolbar.
5
In the Study toolbar, click  Add Study to close the Add Study window.
Parameter Estimation: Armstrong-Frederick Kinematic Hardening
1
In the Settings window for Study, type Parameter Estimation: Armstrong-Frederick Kinematic Hardening in the Label text field.
2
Locate the Study Settings section. Clear the Generate default plots checkbox.
Parameter Estimation
1
In the Study toolbar, click  Optimization and choose Parameter Estimation.
2
In the Settings window for Parameter Estimation, locate the Experimental Data section.
3
From the Data source list, choose All Least-Squares objectives.
4
Locate the Estimated Parameters section. Click  Add four times.
5
6
Locate the Parameter Estimation Method section. From the Method list, choose Levenberg–Marquardt.
7
From the Least-squares time/parameter list method list, choose Use only least-squares data points.
Step 1: Stationary
1
In the Model Builder window, click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Study Extensions section.
3
Select the Auxiliary sweep checkbox.
4
From the Least-squares continuation parameter list, choose t (Time parameter).
Create a plot while solving also for this study.
5
In the Study toolbar, click  Get Initial Value.
Results
Parameter Estimation: Armstrong-Frederick Kinematic Hardening
1
In the Model Builder window, right-click Parameter Estimation: Linear Kinematic Hardening and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Parameter Estimation: Armstrong-Frederick Kinematic Hardening in the Label text field.
3
Locate the Data section. From the Dataset list, choose Parameter Estimation: Armstrong-Frederick Kinematic Hardening/Solution 3 (sol3).
Parameter Estimation: Armstrong-Frederick Kinematic Hardening
Parameter Estimation
1
In the Model Builder window, under Parameter Estimation: Armstrong-Frederick Kinematic Hardening click Parameter Estimation.
2
In the Settings window for Parameter Estimation, locate the Output section.
3
Select the Plot checkbox.
4
5
In the Study toolbar, click  Compute.
The final prediction is now in excellent agreement with the data with an RMS of approximately 2.9.
Results
Collect the final values of the material parameters by creating Evaluation Groups from the Result Templates window.
Result Templates
1
In the Results toolbar, click  Result Templates to open the Result Templates window.
2
Go to the Result Templates window.
3
In the tree, select Parameter Estimation: Linear Kinematic Hardening/Solution 2 (sol2) > Solid Mechanics > Estimated Parameters (std2).
4
Click the Add Result Template button in the window toolbar.
Results
Estimated Parameters: Linear Kinematic Hardening
In the Settings window for Evaluation Group, type Estimated Parameters: Linear Kinematic Hardening in the Label text field.
Result Templates
1
Go to the Result Templates window.
2
In the tree, select Parameter Estimation: Armstrong-Frederick Kinematic Hardening/Solution 3 (sol3) > Solid Mechanics > Estimated Parameters (std3).
3
Click the Add Result Template button in the window toolbar.
4
In the Results toolbar, click  Result Templates to close the Result Templates window.
Results
Estimated Parameters: Armstrong-Frederick Kinematic Hardening
1
In the Settings window for Evaluation Group, type Estimated Parameters: Armstrong-Frederick Kinematic Hardening in the Label text field.
The values of the material parameters can now be copied to a new Material node, which can be used for further analysis.