PDF
Parameter Estimation of Viscoplastic Polymers
Introduction
The Bergstrom–Boyce viscoplastic model for engineering rubbers and soft biological tissues has been successful in capturing nonequilibrium material behavior such as strain-rate dependence and hysteresis under cyclic loading (Ref. 1). In order to accurately predict the time-dependent behavior of such materials using finite element models, the Bergstrom–Boyce model needs to be calibrated against a rich set of experimental data. In this example, we demonstrate how to set up a parameter estimation problem in which the material parameters of the original Bergstrom–Boyce model are calibrated against cyclic uniaxial tension and compression data at two different strain rates. Further, it is shown how additional experimental data obtained at elevated temperatures can be used to include temperature dependence in the model formulation. For further information on the mathematical formulation of the inverse problem and the Bergstrom–Boyce implementation in COMSOL Multiphysics, please refer to the models Parameter Estimation of Hyperelastic Materials and Chloroprene Rubber Compression Test in the Nonlinear Structural Mechanics Module Application Library, respectively.
Model Definition
The original Bergstrom–Boyce model consists of two networks: one hyperelastic Arruda–Boyce network that determines the equilibrium behavior in parallel with a second Arruda–Boyce network in series with a viscoplastic component that captures the nonequilibrium viscoplastic flow. Since the networks act in parallel, the total deformation gradient F is decomposed multiplicatively as , with , , and denoting the elastic deformation gradient in the equilibrium network, the elastic deformation gradient in the nonequilibrium network, and the viscoplastic deformation gradient, respectively.
The strain energy density of the nearly incompressible Arruda–Boyce model reads
(1)
Here, μ0 is the shear modulus of the network; N is the number of chain segments (related to the locking stretch ); cp are known constants derived from a five-term approximation of the inverse Langevin function; and Cel, , and Jel denote the right elastic Cauchy–Green deformation tensor of the network, its first isochoric invariant, and the elastic volume ratio, respectively. The shear moduli of the two networks are related by the energy factor βv, such that . In the nearly incompressible formulation, the volumetric strain energy density Wvol acts as a penalty term on the volume ratio Jel given a large bulk modulus K with respect to the shear modulus of the network.
The viscoplastic flow occurs in the direction of the equivalent stress in the nonequilibrium network, and its magnitude depends on the viscoplastic multiplier λ,
(2)
Herein, A is the viscoplastic rate coefficient (SI unit: 1/s), is the effective chain stretch in the viscoplastic element, ε = 0.001 is a numerical correction factor to avoid division by zero (Ref. 2), is a material parameter controlling the strain hardening, the flow resistance σres is a parameter introduced for dimensional consistency, and n is the stress hardening exponent. In summary, the six unknown material parameters that need to be estimated from experimental data are , N, βv, A, c, and n, see Table 1. It is worth noting that A and σres are dependent parameters; here, we choose to fix σres and only include A in the parameter estimation.
Table 1: Material parameters of the Bergstrom–Boyce model, together with initial values for parameter estimation.
Parameter
Name
Initial guess
Estimate?
Shear modulus, equilibrium network
mu0_eq
1[MPa]
yes
Number of chain segments
Nsegm
5
yes
Bulk modulus
K
1[GPa]
no
Energy factor
beta
1
yes
Viscoplastic rate coefficient
A
1[s^-1]
yes
Strain hardening exponent
c
-0.5
yes
Stress hardening exponent
n
5
yes
Temperature hardening exponent
m
10
yes
Flow resistance
sig_res
sqrt(3)[MPa]
no
In this example, it is assumed that cyclic uniaxial tension and compression data obtained at room temperature (T = Tref = 293.15 K) for two nominal strain rates, 0.001 1/s and 0.1 1/s, are provided. The maximum nominal strains applied in the tension and compression experiments are 0.6 and −0.3, respectively. The stress–strain curves are reported in Figure 1. Note that the data used here were generated using the Bergstrom–Boyce model with parameter values inspired by Ref. 1; we will therefore compare the final calibrated material parameters with those used to generate the data in order to verify the implementation.
Figure 1: Cyclic uniaxial tension and compression data inspired by Ref. 1.
In a second step, suppose that we in addition to the data in Figure 1 also have performed tension experiments at different temperatures. The temperature-dependent tension data is shown in Figure 2. Adding a temperature dependence to the Bergstrom–Boyce model of elastomeric materials can be done by accounting for (i) the entropic stiffening of the polymer network; (ii) the increase in viscoplastic flow at higher temperatures; and (iii) thermal expansion. Here, the latter is assumed negligible. The entropic stiffening is modeled by noting that, based on the statistical theory of rubber elasticity, the shear modulus is proportional to the absolute temperature, so that we can write
(3)
Thermal effects can be included in the expression for the viscoplastic multiplier by the temperature-dependent function g(T), so that Equation 2 becomes
(4)
In this example, a power-law formulation will be used,
(5)
with the temperature exponent m that needs to be estimated from experimental data.
Figure 2: Cyclic tension experiments at three different temperatures.
Results and Discussion
The model prediction for the initial guess of the parameter values in Table 1 is shown in Figure 3. After running the parameter estimation study, the results for the calibrated material model are shown in Figure 4. The fit to the experimental data is excellent, and the final material parameters agree well with those used to generate the data, see Table 2. The same holds for the temperature-dependent cyclic tension data, which is shown together with the final model prediction in Figure 5. Note that, for the values of the material parameters considered here, the viscoplastic flow shows the strongest temperature dependence.
Figure 3: Model prediction with the initial values of the material parameters.
Figure 4: Model prediction with the calibrated material parameters.
Figure 5: Model prediction of the temperature-dependent cyclic tension data.
Table 2: Material parameters of the calibrated Bergstrom–Boyce model.
Parameter
Name
Estimated value
Reference value
Shear modulus, equilibrium network
mu0_eq
0.60013[MPa]
0.6[MPa]
Number of chain segments
Nsegm
8.0152
8.0
Energy factor
beta
2.4981
2.5
Viscoplastic rate coefficient
A
5.0502[s^-1]
5.0[s^-1]
Strain hardening exponent
c
-0.99809
-1.0
Stress hardening exponent
n
4.0004
4.0
Temperature hardening exponent
m
24.971
25.0
Notes About the COMSOL Implementation
In parameter estimation problems, it is good practice to first set up and test the forward model before solving the inverse problem. This is particularly important when the model is highly nonlinear, for which a robust and efficient solver is required. When the experimental data consists of multiple load cases with different boundary conditions, it can be more efficient to solve the load cases in parallel than in series. This is demonstrated here by creating two unit cube elements, one for the compression test and one for the tension test.
The Parameter Estimation functionality is available in COMSOL Multiphysics in the context menu of a Component or under Optimization in the Physics toolbar, wherein each Global Least-Squares Objective node is used to link an experimental data file to the corresponding model variables. To solve the inverse problem, the model needs to be combined with a study containing a Parameter Estimation study step. When multiple objectives are selected in the study step, the total objective function that is minimized will be the sum of all objectives selected.
For most least-squares problems for nonlinear material models, the Levenberg–Marquardt algorithm with a finite difference approximation of the Jacobian is a robust and efficient choice of optimization solver. By default, the Levenberg–Marquardt solver is set to terminate if either the increment of the (scaled) parameters or the maximum angle between the error vector and the Jacobian is smaller than a given optimality tolerance (default 1e-3). In the settings of the Optimization Solver, you can optionally include an additional termination criterion based on the relative change of the objective function by selecting the Terminate also for defect reduction checkbox, which can be useful if the solver reaches a relatively flat local minimum in parameter space where improvements in the objective function are small. The default termination criteria are normally more robust, however, and these will be used here. To monitor the progress of the optimization visually, we will set up a plot while solving that compares the current model prediction with the experimental data.
References
1. J. S. Bergström and M. C. Boyce, “Constitutive modeling of the large strain time-dependent behavior of elastomers,” J. Mech. Phys. Solids, vol. 46, pp. 931–954, 1998.
2. J. S. Bergström and M. C. Boyce, “Constitutive modelling of the time-dependent and cyclic loading of elastomers and application to soft biological tissues,” Mech. Mater., vol. 33, pp. 523–530, 2001.
Application Library path: Nonlinear_Structural_Materials_Module/Viscoplasticity/parameter_estimation_polymer_viscoplasticity
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 > Time Dependent.
6
Click  Done.
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 four data files for the experimental data at 293 K.
Compression, 0.001 1/s, T=293 K
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 Compression, 0.001 1/s, T=293 K in the Label text field.
4
Locate the Data section. Click  Import.
5
Browse to the model’s Application Libraries folder and double-click the file parameter_estimation_polymer_viscoplasticity_compression_1e-3_T293K.txt.
Compression, 0.1 1/s, T=293 K
1
Right-click Compression, 0.001 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Table, type Compression, 0.1 1/s, T=293 K in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file parameter_estimation_polymer_viscoplasticity_compression_1e-1_T293K.txt.
Tension, 0.001 1/s, T=293 K
1
Right-click Compression, 0.1 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Table, type Tension, 0.001 1/s, T=293 K in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file parameter_estimation_polymer_viscoplasticity_tension_1e-3_T293K.txt.
Tension, 0.1 1/s, T=293 K
1
Right-click Tension, 0.001 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Table, type Tension, 0.1 1/s, T=293 K in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file parameter_estimation_polymer_viscoplasticity_tension_1e-1_T293K.txt.
For future use, import also the tension data at 310 K and 323 K.
Tension, 0.1 1/s, T=310 K
1
Right-click Tension, 0.1 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Table, type Tension, 0.1 1/s, T=310 K in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file parameter_estimation_polymer_viscoplasticity_tension_1e-1_T310K.txt.
Tension, 0.1 1/s, T=323 K
1
Right-click Tension, 0.1 1/s, T=310 K and choose Duplicate.
2
In the Settings window for Table, type Tension, 0.1 1/s, T=323 K in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file parameter_estimation_polymer_viscoplasticity_tension_1e-1_T323K.txt.
Stress–Strain Data
Plot the tension and compression data at 293 K.
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Stress-Strain 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 Logarithmic strain (1).
6
Select the y-axis label checkbox. In the associated text field, type Cauchy stress (MPa).
7
Locate the Legend section. From the Position list, choose Upper left.
Compression, 0.001 1/s, T=293 K
1
Right-click Stress–Strain Data and choose Table Graph.
2
In the Settings window for Table Graph, type Compression, 0.001 1/s, T=293 K in the Label text field.
3
Locate the Data section. From the x-axis data list, choose Logarithmic strain (-).
4
From the Plot columns list, choose Manual.
5
In the Columns list box, select Cauchy stress (MPa).
6
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose Dotted.
7
Find the Line markers subsection. From the Marker list, choose Point.
8
Click to expand the Legends section. Select the Show legends checkbox.
9
Find the Include subsection. Clear the Headers checkbox.
10
Select the Label checkbox.
Compression, 0.1 1/s, T=293 K
1
Right-click Compression, 0.001 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Table Graph, type Compression, 0.1 1/s, T=293 K in the Label text field.
3
Locate the Data section. From the Table list, choose Compression, 0.1 1/s, T=293 K.
Tension, 0.001 1/s, T=293 K
1
Right-click Compression, 0.1 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Table Graph, type Tension, 0.001 1/s, T=293 K in the Label text field.
3
Locate the Data section. From the Table list, choose Tension, 0.001 1/s, T=293 K.
Tension, 0.1 1/s, T=293 K
1
Right-click Tension, 0.001 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Table Graph, type Tension, 0.1 1/s, T=293 K in the Label text field.
3
Locate the Data section. From the Table list, choose Tension, 0.1 1/s, T=293 K.
4
In the Stress–Strain Data toolbar, click  Plot.
Global Definitions
Now, set up the forward model. Start by defining the material parameters.
Material Parameters
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, type Material Parameters in the Label text field.
3
Locate the Parameters section. In the table, enter the following settings:
Name
Expression
Value
Description
mu0_eq
1[MPa]
1 MPa
Shear modulus, equilibrium network
Nsegm
5
5
Number of chain segments
K
1[GPa]
1000 MPa
Bulk modulus
beta
1
1
Energy factor
A
1[s^-1]
1 1/s
Viscoplastic rate coefficient
c
-0.5
-0.5
Strain hardening exponent
n
5
5
Stress hardening exponent
m
10
10
Temperature hardening exponent
sig_res
sqrt(3)[MPa]
1.7321 MPa
Flow resistance
Load Parameters
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Load Parameters in the Label text field.
3
Locate the Parameters section. In the table, enter the following settings:
Name
Expression
Value
Description
strain_rate
1e-3[s^-1]
0.001 1/s
Nominal strain rate
t_end
1/strain_rate
1000 s
Loading duration
emax_ten
0.6
0.6
Nominal strain, tension load case
emax_comp
-0.3
-0.3
Nominal strain, compression load case
Tref
293.15[K]
293.15 K
Reference temperature
T
293.15[K]
293.15 K
Temperature
Load Cycle
Add a Triangle function for prescribing the cyclic loading. The function will be given in terms of dimensionless units so that it can be used for both the tension and compression load case, independent of the strain rate.
1
In the Home toolbar, click  Functions and choose Global > Triangle.
2
In the Settings window for Triangle, type Load Cycle in the Label text field.
3
Locate the Parameters section. In the Lower limit text field, type 0.
4
In the Upper limit text field, type 2.
5
Click to expand the Smoothing section. Clear the Size of transition zone checkbox.
Geometry 1
Create two unit cubes, one for the compression and one for the tension load case.
1
In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2
In the Settings window for Geometry, locate the Units section.
3
From the Length unit list, choose mm.
Block 1 (blk1)
In the Geometry toolbar, click  Block.
Array 1 (arr1)
1
In the Geometry toolbar, click  Transforms and choose Array.
2
Select the object blk1 only.
3
In the Settings window for Array, locate the Size section.
4
In the y size text field, type 2.
5
Locate the Displacement section. In the y text field, type 2.
6
In the Home toolbar, click  Build All.
7
Click the  Go to Default View button in the Graphics toolbar.
Materials
Both the Arruda–Boyce and Bergstrom–Boyce Viscoplasticity material property groups need to be added for defining all material parameters in the Bergstrom–Boyce constitutive model.
Bergstrom–Boyce Material
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
In the Settings window for Material, type Bergstrom-Boyce Material in the Label text field.
3
Click to expand the Material Properties section. In the Material properties tree, select Solid Mechanics > Hyperelastic Material > Arruda–Boyce.
4
Click  Add to Material.
5
In the Material properties tree, select Solid Mechanics > Viscoplastic Material > Bergstrom–Boyce Viscoplasticity.
6
Click  Add to Material.
7
In the Model Builder window, expand the Component 1 (comp1) > Materials > Bergstrom–Boyce Material (mat1) node, then click Arruda–Boyce (ArrudaBoyce).
8
In the Settings window for Arruda–Boyce, locate the Output Properties section.
9
In the table, enter the following settings:
Property
Variable
Expression
Unit
Size
Number of segments
Nseg
Nsegm
1
1x1
Macroscopic shear modulus
mu0
mu0_eq*(1+(T-Tref)/Tref)
MPa
1x1
10
In the Model Builder window, under Component 1 (comp1) > Materials > Bergstrom–Boyce Material (mat1) click Bergstrom–Boyce viscoplasticity (BergstromBoyce).
11
In the Settings window for Bergstrom–Boyce Viscoplasticity, locate the Output Properties section.
12
In the table, enter the following settings:
Property
Variable
Expression
Unit
Size
Viscoplastic rate coefficient
A_BB
A
1/s
1x1
Flow resistance
sigRes_BB
sig_res
MPa
1x1
Stress exponent
n_BB
n
1
1x1
Strain exponent
c_BB
c
1
1x1
Cutoff stress
sigmaco_BB
0[MPa]
MPa
1x1
Solid Mechanics (solid)
Continue to set up the Bergstrom–Boyce material model and the uniaxial tension and compression load cases under the Solid Mechanics interface. Because the states of stress and strain will be homogeneous, we can use linear shape functions and reduced integration to reduce the computational cost of the forward model.
1
In the Model Builder window, under Component 1 (comp1) click Solid Mechanics (solid).
2
In the Settings window for Solid Mechanics, locate the Structural Transient Behavior section.
3
From the list, choose Quasistatic.
4
Click to expand the Discretization section. From the Displacement field list, choose Linear.
Hyperelastic Material 1
1
In the Physics toolbar, click  Domains and choose Hyperelastic Material.
2
Click in the Graphics window and then press Ctrl+A to select both domains.
3
In the Settings window for Hyperelastic Material, locate the Hyperelastic Material section.
4
From the Material model list, choose Arruda–Boyce.
5
From the Volumetric strain energy list, choose Miehe.
6
In the κ text field, type K.
7
Locate the Quadrature Settings section. Select the Reduced integration checkbox.
Polymer Viscoplasticity 1
1
In the Physics toolbar, click  Attributes and choose Polymer Viscoplasticity.
2
In the Settings window for Polymer Viscoplasticity, locate the Viscoplasticity Model section.
3
Find the Hyperelastic element subsection. In the βv text field, type beta.
4
Locate the Thermal Effects section. From the g  (T) list, choose Power law.
5
In the Tref text field, type Tref.
6
In the m text field, type m.
7
Locate the Time Stepping section. From the Method list, choose Domain ODEs.
When the number of degrees of freedom is small, formulating the viscoplastic rate equations in terms of domain ODEs yields better performance than solving them locally.
8
Locate the Model Input section. From the T list, choose User defined. In the associated text field, type T.
Roller 1
Add ideal boundary conditions for the two load cases.
1
In the Physics toolbar, click  Boundaries and choose Roller.
2
Select Boundaries 1–3 and 6–8 only.
Uniaxial Compression
1
In the Physics toolbar, click  Boundaries and choose Prescribed Displacement.
2
In the Settings window for Prescribed Displacement, type Uniaxial Compression in the Label text field.
3
Select Boundary 11 only.
4
Locate the Prescribed Displacement section. From the Displacement in x direction list, choose Prescribed.
5
In the u 0 x text field, type emax_comp*tri1(t/t_end)*1[mm].
Uniaxial Tension
1
Right-click Uniaxial Compression and choose Duplicate.
2
In the Settings window for Prescribed Displacement, type Uniaxial Tension in the Label text field.
3
Select Boundary 12 only.
4
Locate the Prescribed Displacement section. In the u 0 x text field, type emax_ten*tri1(t/t_end)*1[mm].
Mesh 1
Mesh each domain with a single hexahedral element.
Mapped 1
1
In the Mesh toolbar, click  More Generators and choose Mapped.
2
Select Boundaries 1 and 6 only.
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.
Definitions
Before setting up the solver, add global variables for the volume-averaged Cauchy stress and logarithmic strain in the two load cases. These will be used for the comparison with experimental data.
Average 1 (aveop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Average.
2
Select Domain 1 only.
Average 2 (aveop2)
1
Right-click Average 1 (aveop1) and choose Duplicate.
2
Select Domain 2 only.
Global Stress and Strain Variables
1
In the Model Builder window, right-click Definitions and choose Variables.
2
In the Settings window for Variables, type Global Stress and Strain Variables in the Label text field.
3
Locate the Variables section. In the table, enter the following settings:
Name
Expression
Unit
Description
sig_comp
aveop1(solid.sxx)
MPa
Cauchy stress, compression
elog_comp
aveop1(solid.elogxx)
 
Logarithmic strain, compression
sig_ten
aveop2(solid.sxx)
MPa
Cauchy stress, tension
elog_ten
aveop2(solid.elogxx)
 
Logarithmic strain, tension
Forward Problem
Set up a first study to verify that the forward model is defined correctly. A parametric sweep is used to compute the response for the two strain rates.
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, type Forward Problem in the Label text field.
3
Locate the Study Settings section. Clear the Generate default plots checkbox.
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
Click  Add.
4
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
strain_rate (Nominal strain rate)
0.001 0.1
1/s
Step 1: Time Dependent
1
In the Model Builder window, click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,0.02,1)*2*t_end.
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
Add manual scalings to the dependent variables to improve the convergence rate.
2
In the Model Builder window, expand the Solution 1 (sol1) node.
3
In the Model Builder window, expand the Forward Problem > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 node, then click Viscoplastic Strain Tensor, Local Coordinate System (comp1.solid.hmm1.pvp1.evp).
4
In the Settings window for Field, locate the Scaling section.
5
In the Scale text field, type 1.
6
In the Model Builder window, under Forward Problem > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Equivalent Viscoplastic Strain (comp1.solid.hmm1.pvp1.evpe).
7
In the Settings window for Field, locate the Scaling section.
8
In the Scale text field, type 1.
9
In the Model Builder window, under Forward Problem > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Auxiliary Pressure (comp1.solid.hmm1.pw).
10
In the Settings window for Field, locate the Scaling section.
11
In the Scale text field, type 1[MPa].
12
In the Model Builder window, under Forward Problem > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Displacement Field (comp1.u).
13
In the Settings window for Field, locate the Scaling section.
14
In the Scale text field, type emax_ten.
15
In the Model Builder window, under Forward Problem > Solver Configurations > Solution 1 (sol1) click Time-Dependent Solver 1.
16
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
17
From the Steps taken by solver list, choose Strict.
When solving inelastic problems, strict time stepping is often required to avoid interpolation errors.
18
Find the Algebraic variable settings subsection. From the Consistent initialization list, choose Off.
The boundary and initial conditions in the model are consistent, so Consistent initialization can be disabled.
19
In the Model Builder window, under Forward Problem > Solver Configurations > Solution 1 (sol1) > Time-Dependent Solver 1 click Fully Coupled 1.
20
In the Settings window for Fully Coupled, click to expand the Method and Termination section.
21
From the Nonlinear method list, choose Constant (Newton).
The constant Newton solver is a fast and robust choice for nonlinear problems.
22
From the Jacobian update list, choose Once per time step.
23
In the Maximum number of iterations text field, type 25.
24
In the Study toolbar, click  Compute.
Results
Compare the initial model prediction with the experimental data.
Initial Model Prediction
1
In the Model Builder window, right-click Stress–Strain 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/Parametric Solutions 1 (sol2).
Initial Prediction, Compression
1
Right-click Initial Model Prediction and choose Global.
2
In the Settings window for Global, type Initial Prediction, Compression in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
sig_comp
MPa
Cauchy stress, compression
4
Locate the x-Axis Data section. From the Parameter list, choose Expression.
5
In the Expression text field, type elog_comp.
6
Click to expand the Coloring and Style section. From the Color list, choose Cycle (reset).
7
Click to expand the Legends section. Clear the Show legends checkbox.
Initial Prediction, Tension
1
Right-click Initial Prediction, Compression and choose Duplicate.
2
In the Settings window for Global, type Initial Prediction, Tension in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
sig_ten
MPa
Cauchy stress, tension
4
Locate the x-Axis Data section. In the Expression text field, type elog_ten.
5
Locate the Coloring and Style section. From the Color list, choose Cycle.
6
In the Initial Model Prediction toolbar, click  Plot.
Component 1 (comp1)
Now, move on to the inverse problem. For this, we need to define four Least-Squares Objectives, one for each experiment at 293 K.
Compression, 0.001 1/s, T=293 K
1
In the Physics toolbar, click  Optimization and choose Parameter Estimation.
2
In the Settings window for Least-Squares Objective, type Compression, 0.001 1/s, T=293 K in the Label text field.
3
Locate the Experimental Data section. From the Data source list, choose Result table.
In the Column Settings section, we need to provide definitions for each column in the data file. Here, the first column contains the time, the second contains the logarithmic strain, and the third contains Cauchy stress data. Since the tests are displacement controlled, we will only use the stress data for the evaluation of the objective function and assume that the strain is measured with negligible error. The Logarithmic strain column is therefore ignored.
4
Locate the Data Column Settings section. In the table, enter the following settings:
Columns
Type
Settings
Logarithmic strain (-)
Ignored column
 
5
In the table, click to select the cell at row number 3 and column number 1.
In the Model expression field, enter the name of the global variable we defined earlier for the Cauchy stress in the compression load case.
6
In the Model expression text field, type comp1.sig_comp.
7
In the Column name text field, type comp_slow.
The Variable name is used to refer to variables for the data and the model expression in postprocessing.
8
In the Unit text field, type MPa.
In the Experimental Conditions section, specify the value of the model parameters that are constant for this particular dataset, but varies between datasets. In this example, we want to add the parameters strain_rate and T and set them to the values corresponding to the experiment. Since the temperature in these four experiments is constant and equal to the value given under Global Definitions, we do not need to specify it explicitly here.
9
Locate the Experimental Conditions section. Click  Add.
10
In the table, enter the following settings:
Name
Expression
strain_rate (Nominal strain rate)
0.001[s^-1]
Compression, 0.1 1/s, T=293 K
Continue in a similar fashion with the remaining three datasets.
1
Right-click Compression, 0.001 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Least-Squares Objective, type Compression, 0.1 1/s, T=293 K in the Label text field.
3
Locate the Experimental Data section. From the Result table list, choose Compression, 0.1 1/s, T=293 K.
4
Locate the Data Column Settings section. In the table, click to select the cell at row number 3 and column number 1.
5
In the Column name text field, type comp_fast.
6
Locate the Experimental Conditions section. In the table, enter the following settings:
Name
Expression
strain_rate (Nominal strain rate)
0.1[s^-1]
Tension, 0.001 1/s, T=293 K
1
Right-click Compression, 0.1 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Least-Squares Objective, type Tension, 0.001 1/s, T=293 K in the Label text field.
3
Locate the Experimental Data section. From the Result table list, choose Tension, 0.001 1/s, T=293 K.
4
Locate the Data Column Settings section. In the table, click to select the cell at row number 3 and column number 1.
5
In the Model expression text field, type comp1.sig_ten.
6
In the Column name text field, type ten_slow.
7
Locate the Experimental Conditions section. In the table, enter the following settings:
Name
Expression
strain_rate (Nominal strain rate)
0.001[s^-1]
Tension, 0.1 1/s, T=293 K
1
Right-click Tension, 0.001 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Least-Squares Objective, type Tension, 0.1 1/s, T=293 K in the Label text field.
3
Locate the Experimental Data section. From the Result table list, choose Tension, 0.1 1/s, T=293 K.
4
Locate the Data Column Settings section. In the table, click to select the cell at row number 3 and column number 1.
5
In the Column name text field, type ten_fast.
6
Locate the Experimental Conditions section. In the table, enter the following settings:
Name
Expression
strain_rate (Nominal strain rate)
0.1[s^-1]
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 > Time Dependent.
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
1
In the Settings window for Study, type Parameter Estimation 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 Selected Least-Squares objectives.
4
Locate the Objective Function section. In the table, select the Active checkboxes for Compression, 0.001 1/s, T=293 K, Compression, 0.1 1/s, T=293 K, Tension, 0.001 1/s, T=293 K, and Tension, 0.1 1/s, T=293 K.
5
Locate the Estimated Parameters section. Click  Add six times.
6
In the table, enter the following settings:
Parameter
Initial value
Scale
Lower bound
Upper bound
Unit
mu0_eq (Shear modulus, equilibrium network)
1[MPa]
1
0
 
MPa
Nsegm (Number of chain segments)
5
1
2
 
 
beta (Energy factor)
1
1
0
 
 
A (Viscoplastic rate coefficient)
1[s^-1]
1
0
 
1/s
c (Strain hardening exponent)
-0.5
1
-1
0
 
n (Stress hardening exponent)
5
1
1.01
20
 
7
Locate the Parameter Estimation Method section. From the Method list, choose Levenberg–Marquardt.
8
From the Least-squares time/parameter list method list, choose Use only least-squares data points.
Solution 5 (sol5)
Generate the solver sequence and edit the default solver settings as we did for the first study.
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 5 (sol5) node.
3
In the Model Builder window, expand the Parameter Estimation > Solver Configurations > Solution 5 (sol5) > Dependent Variables 1 node, then click Viscoplastic Strain Tensor, Local Coordinate System (comp1.solid.hmm1.pvp1.evp).
4
In the Settings window for Field, locate the Scaling section.
5
In the Scale text field, type 1.
6
In the Model Builder window, under Parameter Estimation > Solver Configurations > Solution 5 (sol5) > Dependent Variables 1 click Equivalent Viscoplastic Strain (comp1.solid.hmm1.pvp1.evpe).
7
In the Settings window for Field, locate the Scaling section.
8
In the Scale text field, type 1.
9
In the Model Builder window, under Parameter Estimation > Solver Configurations > Solution 5 (sol5) > Dependent Variables 1 click Auxiliary Pressure (comp1.solid.hmm1.pw).
10
In the Settings window for Field, locate the Scaling section.
11
In the Scale text field, type 1.
12
In the Model Builder window, under Parameter Estimation > Solver Configurations > Solution 5 (sol5) > Dependent Variables 1 click Displacement Field (comp1.u).
13
In the Settings window for Field, locate the Scaling section.
14
In the Scale text field, type emax_ten.
15
In the Model Builder window, expand the Parameter Estimation > Solver Configurations > Solution 5 (sol5) > Optimization Solver 1 node, then click Time-Dependent Solver 1.
16
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
17
From the Steps taken by solver list, choose Strict.
18
Find the Algebraic variable settings subsection. From the Consistent initialization list, choose Off.
19
In the Model Builder window, expand the Parameter Estimation > Solver Configurations > Solution 5 (sol5) > Optimization Solver 1 > Time-Dependent Solver 1 node, then click Fully Coupled 1.
20
In the Settings window for Fully Coupled, locate the Method and Termination section.
21
From the Nonlinear method list, choose Constant (Newton).
22
From the Jacobian update list, choose Once per time step.
23
In the Maximum number of iterations text field, type 25.
Results
Before solving the study, create a plot comparing the model prediction with the experimental data to monitor the progress of the optimization solver visually. Start by copying the plot group Stress–Strain Data that already contains all experimental data curves.
Parameter Estimation: Bergstrom–Boyce Viscoplasticity
1
In the Model Builder window, right-click Stress–Strain Data and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Parameter Estimation: Bergstrom-Boyce Viscoplasticity in the Label text field.
3
Locate the Data section. From the Dataset list, choose Parameter Estimation/Solution 5 (sol5).
Model Prediction, Compression 0.001 1/s
1
Right-click Parameter Estimation: Bergstrom–Boyce Viscoplasticity and choose Global.
2
In the Settings window for Global, type Model Prediction, Compression 0.001 1/s in the Label text field.
Each Least-Squares Objective feature defines a variable lso.variable_name.model, lso being the tag of the feature and variable_name the string entered in the Variable name field of the Value column, that contains the corresponding model expression.
3
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
lso1.comp_slow.model
MPa
Least-squares model value
4
Locate the x-Axis Data section. From the Parameter list, choose Expression.
5
In the Expression text field, type elog_comp.
6
Locate the Coloring and Style section. From the Color list, choose Cycle (reset).
7
Locate the Legends section. Clear the Show legends checkbox.
Model Prediction, Compression 0.1 1/s
1
Right-click Model Prediction, Compression 0.001 1/s and choose Duplicate.
2
In the Settings window for Global, type Model Prediction, Compression 0.1 1/s in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
lso2.comp_fast.model
MPa
Least-squares model value
4
Locate the Coloring and Style section. From the Color list, choose Cycle.
Model Prediction, Tension 0.001 1/s
1
Right-click Model Prediction, Compression 0.1 1/s and choose Duplicate.
2
In the Settings window for Global, type Model Prediction, Tension 0.001 1/s in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
lso3.ten_slow.model
MPa
Least-squares model value
4
Locate the x-Axis Data section. In the Expression text field, type elog_ten.
Model Prediction, Tension 0.1 1/s
1
Right-click Model Prediction, Tension 0.001 1/s and choose Duplicate.
2
In the Settings window for Global, type Model Prediction, Tension 0.1 1/s in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
lso4.ten_fast.model
MPa
Least-squares model value
Parameter Estimation
Parameter Estimation
1
In the Model Builder window, under Parameter Estimation click Parameter Estimation.
2
In the Settings window for Parameter Estimation, click to expand the Output section.
3
Select the Plot checkbox.
4
In the table, enter the following settings:
Plot group
Plot window
Parameter Estimation: Bergstrom-Boyce Viscoplasticity
Graphics
5
Select the Show individual objective values checkbox.
6
Select the Table graph checkbox.
7
In the Study toolbar, click  Compute.
Solid Mechanics (solid)
The calibrated material parameters of the Bergstrom–Boyce model can now be copied from the Objective Probe Table, or accessed within the model using the withsol operator. We will use the latter to create a new Hyperelastic Material with the calibrated material properties. This is needed for the next step where we will estimate the temperature dependence.
Hyperelastic Material 2
1
In the Model Builder window, under Component 1 (comp1) > Solid Mechanics (solid) right-click Hyperelastic Material 1 and choose Duplicate.
2
In the Settings window for Hyperelastic Material, locate the Hyperelastic Material section.
3
From the μ0 list, choose User defined. In the associated text field, type withsol('sol5', mu0_eq)*(1+(T-Tref)/Tref).
4
From the N list, choose User defined. In the associated text field, type withsol('sol5', Nsegm).
Polymer Viscoplasticity 1
1
In the Model Builder window, expand the Hyperelastic Material 2 node, then click Polymer Viscoplasticity 1.
2
In the Settings window for Polymer Viscoplasticity, locate the Viscoplasticity Model section.
3
Find the Hyperelastic element subsection. In the βv text field, type withsol('sol5', beta).
4
Find the Inelastic element subsection. From the A list, choose User defined. In the associated text field, type withsol('sol5', A).
5
From the σres list, choose User defined. In the associated text field, type sig_res.
6
From the n list, choose User defined. In the associated text field, type withsol('sol5', n).
7
From the σco list, choose User defined. Find the Isotropic hardening model subsection. From the c list, choose User defined. In the associated text field, type withsol('sol5', c).
Results
Plot the temperature-dependent tension data.
Stress–Strain Data: Temperature Dependence
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Stress-Strain Data: Temperature Dependence in the Label text field.
3
Locate the Data section. From the Dataset list, choose None.
4
Locate the Title section. From the Title type list, choose Label.
5
Locate the Legend section. From the Position list, choose Upper left.
Tension, 0.1 1/s, T=293 K
1
Right-click Stress–Strain Data: Temperature Dependence and choose Table Graph.
2
In the Settings window for Table Graph, locate the Data section.
3
From the Table list, choose Tension, 0.1 1/s, T=293 K.
4
From the x-axis data list, choose Logarithmic strain (-).
5
From the Plot columns list, choose Manual.
6
In the Columns list box, select Cauchy stress (MPa).
7
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose Dotted.
8
Find the Line markers subsection. From the Marker list, choose Point.
9
Locate the Legends section. Select the Show legends checkbox.
10
Find the Include subsection. Select the Label checkbox.
11
Clear the Headers checkbox.
12
In the Label text field, type Tension, 0.1 1/s, T=293 K.
Tension, 0.1 1/s, T=310 K
1
Right-click Tension, 0.1 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Table Graph, type Tension, 0.1 1/s, T=310 K in the Label text field.
3
Locate the Data section. From the Table list, choose Tension, 0.1 1/s, T=310 K.
Tension, 0.1 1/s, T=323 K
1
Right-click Tension, 0.1 1/s, T=310 K and choose Duplicate.
2
In the Settings window for Table Graph, type Tension, 0.1 1/s, T=323 K in the Label text field.
3
Locate the Data section. From the Table list, choose Tension, 0.1 1/s, T=323 K.
4
In the Stress–Strain Data: Temperature Dependence toolbar, click  Plot.
Parameter Estimation
Add new Least-Squares Objectives and a Parameter Estimation study to fit the temperature dependence.
Tension, 0.1 1/s, T=310 K
1
In the Model Builder window, right-click Tension, 0.1 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Least-Squares Objective, type Tension, 0.1 1/s, T=310 K in the Label text field.
3
Locate the Experimental Data section. From the Result table list, choose Tension, 0.1 1/s, T=310 K.
4
Locate the Data Column Settings section. In the table, click to select the cell at row number 3 and column number 1.
5
In the Column name text field, type ten_fast_310K.
6
Locate the Experimental Conditions section. Click  Add.
7
In the table, enter the following settings:
Name
Expression
T (Temperature)
310.15[K]
Tension, 0.1 1/s, T=323 K
1
Right-click Tension, 0.1 1/s, T=310 K and choose Duplicate.
2
In the Settings window for Least-Squares Objective, type Tension, 0.1 1/s, T=323 K in the Label text field.
3
Locate the Experimental Data section. From the Result table list, choose Tension, 0.1 1/s, T=323 K.
4
Locate the Data Column Settings section. In the table, click to select the cell at row number 3 and column number 1.
5
In the Column name text field, type ten_fast_323K.
6
Locate the Experimental Conditions section. In the table, enter the following settings:
Name
Expression
T (Temperature)
323.15[K]
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 > Time Dependent.
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: Temperature Dependence
1
In the Settings window for Study, type Parameter Estimation: Temperature Dependence 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.
The only remaining unknown material parameter is the temperature hardening exponent m. To estimate its value, it is sufficient to use the last three objectives.
2
In the Settings window for Parameter Estimation, locate the Experimental Data section.
3
From the Data source list, choose Selected Least-Squares objectives.
4
Locate the Objective Function section. In the table, enter the following settings:
Objective functions from components
Active
Compression, 0.001 1/s, T=293 K
 
Compression, 0.1 1/s, T=293 K
 
Tension, 0.001 1/s, T=293 K
 
Tension, 0.1 1/s, T=293 K
√
Tension, 0.1 1/s, T=310 K
√
Tension, 0.1 1/s, T=323 K
√
5
Locate the Estimated Parameters section. Click  Add.
6
In the table, enter the following settings:
Parameter
Initial value
Scale
Lower bound
Upper bound
Unit
m (Temperature hardening exponent)
10
10
 
 
 
7
Locate the Parameter Estimation Method section. From the Method list, choose Levenberg–Marquardt.
8
From the Least-squares time/parameter list method list, choose Use only least-squares data points.
Step 1: Time Dependent
1
In the Model Builder window, click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Physics and Variables Selection section.
3
Select the Modify model configuration for study step checkbox.
4
In the tree, select Component 1 (comp1) > Solid Mechanics (solid), Controls spatial frame > Hyperelastic Material 1.
5
Click  Disable.
Solution 6 (sol6)
Generate the solver sequence and edit the default solver settings as we did for the first study, but modify the finite difference interval for the Jacobian computation in the Optimization Solver to 1e-3. Using a larger finite difference step can be necessary if the sensitivity to the estimated parameters is low. Refer to the previous study for the exact step-by-step instructions.
1
In the Study toolbar, click  Show Default Solver.
2
Click  Compute.
Results
Plot the model prediction of the temperature-dependent data.
Parameter Estimation: Temperature Dependence
1
In the Model Builder window, right-click Stress–Strain Data: Temperature Dependence and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Parameter Estimation: Temperature Dependence in the Label text field.
3
Locate the Data section. From the Dataset list, choose Parameter Estimation: Temperature Dependence/Solution 6 (sol6).
Model Prediction, Tension 0.1 1/s, T=293 K
1
Right-click Parameter Estimation: Temperature Dependence and choose Global.
2
In the Settings window for Global, type Model Prediction, Tension 0.1 1/s, T=293 K in the Label text field.
3
Locate the Data section. From the Dataset list, choose Parameter Estimation: Temperature Dependence/Solution 6 (sol6).
4
From the Parameter selection (T) list, choose From list.
5
In the Parameter values (T) list box, select 293.15.
6
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
lso4.ten_fast.model
MPa
Least-squares model value
7
Locate the x-Axis Data section. From the Parameter list, choose Expression.
8
In the Expression text field, type elog_ten.
9
Locate the Coloring and Style section. From the Color list, choose Cycle (reset).
10
Locate the Legends section. Clear the Show legends checkbox.
Model Prediction, Tension 0.1 1/s, T=310 K
1
Right-click Model Prediction, Tension 0.1 1/s, T=293 K and choose Duplicate.
2
In the Settings window for Global, type Model Prediction, Tension 0.1 1/s, T=310 K in the Label text field.
3
Locate the Data section. In the Parameter values (T) list box, select 310.15.
4
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
lso5.ten_fast_310K.model
MPa
Least-squares model value
5
Locate the Coloring and Style section. From the Color list, choose Cycle.
Model Prediction, Tension 0.1 1/s, T=323 K
1
Right-click Model Prediction, Tension 0.1 1/s, T=310 K and choose Duplicate.
2
In the Settings window for Global, type Model Prediction, Tension 0.1 1/s, T=323 K in the Label text field.
3
Locate the Data section. In the Parameter values (T) list box, select 323.15.
4
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
lso6.ten_fast_323K.model
MPa
Least-squares model value
5
In the Parameter Estimation: Temperature Dependence toolbar, click  Plot.
Result Templates
Finish up by collecting the complete set of calibrated material parameters in an Evaluation Group.
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/Solution 5 (sol5) > Solid Mechanics > Estimated Parameters (std2).
4
Click the Add Result Template button in the window toolbar.
5
In the Results toolbar, click  Result Templates to close the Result Templates window.
Results
Calibrated Material Parameters
In the Settings window for Evaluation Group, type Calibrated Material Parameters in the Label text field.
Global Evaluation 2
1
Right-click Calibrated Material Parameters and choose Global Evaluation.
2
In the Settings window for Global Evaluation, locate the Data section.
3
From the Dataset list, choose Parameter Estimation: Temperature Dependence/Solution 6 (sol6).
4
From the Parameter selection (T) list, choose Last.
5
From the Parameter selection (strain_rate) list, choose Last.
6
From the Time selection list, choose Last.
7
Locate the Expressions section. In the table, enter the following settings:
Expression
Unit
Description
m
 
Temperature hardening exponent
8
In the Calibrated Material Parameters toolbar, click  Evaluate.