PDF

Parameter Estimation of Hyperelastic Materials
Introduction
Rubber-like components in tires, seals, insulators, soft sensors, and actuators are often modeled as hyperelastic materials (Ref. 1). In order to accurately predict the behavior of such components using finite element models, the material model chosen needs to be calibrated and validated against experimental data. In contrast to isotropic linear elastic materials, for which the Young’s modulus and the Poisson’s ratio can be estimated directly from the stress and lateral contraction measured in a uniaxial tension test, the calibration of hyperelastic material models typically requires consideration of multiple load cases over the full range of deformation expected in the final application. This tutorial model demonstrates how to set up this so-called inverse problem in order to estimate the material parameters of a hyperelastic model from experimental data. The data are reproduced from Ref. 2, wherein large deformation uniaxial tension, pure shear, and equibiaxial tension tests were performed on a soft elastomeric material employed in a tactile sensor. The procedure is validated by comparing the material parameters of an Ogden hyperelastic model against the results reported in Ref. 2.
Model Definition
Parameter estimation problems consist of three components: (i) experimental data; (ii) a forward model that represents the physics of the experiments; and (iii) an optimization algorithm that compares the two and updates the model parameters to minimize the difference. This can be formulated mathematically as a nonlinear least-squares minimization problem,
(1)
with
(2)
Herein, q is the vector of material parameters that we want to estimate, N is the number of experiments, Mn is the number of data points per experiment, is the mth data point of experiment n, and Pnm, q) denotes the corresponding model prediction given the experimental parameter λm.
In this example, we consider N = 3 experiments (uniaxial tension, pure shear, and equibiaxial tension), for which the measured quantity Pn is the first Piola–Kirchhoff stress and λm is the applied stretch in the loading direction. A schematic of the three experiments is shown in Figure 1. The data from Ref. 2 are reproduced in Figure 2.
Figure 1: Illustration of the three homogeneous load cases considered.
Figure 2: Experimental data from Ref. 2.
In Ref. 2, these data were used to fit the parameters of an incompressible second-order Ogden model. This form of the Ogden strain energy density function reads
(3)
which consists of four unknown material parameters that we will estimate, that is, q = (μ1α1μ2α2). Note that for the model to yield physically admissible predictions, the parameters need to satisfy the constraints μpαp > 0 for all terms p. The parameters along with an initial guess of their values are provided in Table 1.
It is important to note that experimental data from standardized material tests are usually given in terms of stress–strain curves for a homogeneous state of stress and deformation. In this case, the forward model can be set up with a single element, reduced integration, and idealized boundary conditions. This reduces the computational cost significantly compared to solving the full physical problem for every model evaluation within the optimization solver. However, if the assumption of homogeneity does not hold for the experimental data at hand, you may have to perform the inverse analysis by simulating the full geometry and comparing the global force–displacement curve instead of stress–strain data.
Results and Discussion
The model prediction for the initial guess of the parameter values in Table 1 is shown in Figure 3. Note that the uniaxial and pure shear behavior is of the correct order of magnitude, but the equibiaxial response is off by an order of magnitude for large stretches.
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 with those reported in Ref. 2, see Table 2.
Figure 3: Model prediction with the initial values of the material parameters.
Figure 4: Model prediction with the calibrated material parameters. Note that the equibiaxial curves are displayed on the second y-axis to better visualize all datasets in a single plot.
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. When the experimental data consists of multiple load cases with different boundary conditions, it can be more efficient to solve them in parallel than in series. This is demonstrated here by creating three unit cube elements, one for each experiment.
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 Least-Squares Objective node adds an objective corresponding to Equation 2 to the model. To solve the inverse problem, these need 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, see Equation 1. 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.
References
1. P. Steinmann, M. Hossain, and G. Possart, “Hyperelastic models for rubber-like materials: consistent tangent operators and suitability for Treloar’s data,” Arch. Appl. Mech., vol. 82, pp. 1183–1217, 2012.
2. C. Sferrazza, A. Wahlsten, C. Trueeb, and R. d’Andrea, “Ground Truth Force Distribution for Learning-Based Tactile Sensing: A Finite Element Approach,” IEEE Access, vol. 7, pp. 173438–173449, 2019.
Application Library path: Nonlinear_Structural_Materials_Module/Hyperelasticity/parameter_estimation_hyperelasticity
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 experimental data files to result tables.
Uniaxial 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 Uniaxial Data in the Label text field.
4
Locate the Data section. Click  Import.
5
Pure Shear Data
1
In the Results toolbar, click  Table.
2
In the Settings window for Table, type Pure Shear Data in the Label text field.
3
Locate the Data section. Click  Import.
4
Equibiaxial Data
1
In the Results toolbar, click  Table.
2
In the Settings window for Table, type Equibiaxial Data 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_hyperelasticity_equibiaxial.txt.
Experimental Data
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Experimental 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 Applied stretch (1).
6
Select the y-axis label checkbox. In the associated text field, type Nominal stress (MPa).
7
Locate the Legend section. From the Position list, choose Upper left.
Uniaxial Data
1
Right-click Experimental Data and choose Table Graph.
2
In the Settings window for Table Graph, type Uniaxial Data in the Label text field.
3
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose None.
4
From the Color list, choose Cycle (reset).
5
Find the Line markers subsection. From the Marker list, choose Point.
6
Click to expand the Legends section. Select the Show legends checkbox.
7
From the Legends list, choose Manual.
8
Pure Shear Data
1
Right-click Uniaxial Data and choose Duplicate.
2
In the Settings window for Table Graph, type Pure Shear Data in the Label text field.
3
Locate the Data section. From the Table list, choose Pure Shear Data.
4
Locate the Coloring and Style section. From the Color list, choose Cycle.
5
Locate the Legends section. In the table, enter the following settings:
Equibiaxial Data
1
Right-click Pure Shear Data and choose Duplicate.
2
In the Settings window for Table Graph, type Equibiaxial Data in the Label text field.
3
Locate the Data section. From the Table list, choose Equibiaxial Data.
4
Locate the Legends section. In the table, enter the following settings:
5
In the Experimental Data toolbar, click  Plot.
Global Definitions
Continue with defining the material parameters and setting up the forward problem.
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Geometry 1
Create three unit cubes, one for each 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
3
In the Settings window for Array, locate the Size section.
4
In the y size text field, type 3.
5
Locate the Displacement section. In the y text field, type 3.
6
In the Geometry toolbar, click  Build All.
7
Click the  Zoom Extents button in the Graphics toolbar.
Solid Mechanics (solid)
Set up the three load cases in the Solid Mechanics interface. Since they result in homogeneous stresses and deformations, use linear shape functions with reduced integration to reduce the computation cost.
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
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 all domains.
3
In the Settings window for Hyperelastic Material, locate the Hyperelastic Material section.
4
From the Material model list, choose Ogden.
5
From the Compressibility list, choose Incompressible.
6
7
In the Ogden parameters table, enter the following settings:
8
Locate the Quadrature Settings section. Select the Reduced integration checkbox.
Roller 1
Add symmetry boundary conditions with the Roller node to suppress rigid body motions.
1
In the Physics toolbar, click  Boundaries and choose Roller.
2
Prescribed Displacement 1
The displacement in the x direction is identical for all load cases.
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 stretch-1.
Prescribed Displacement 2
Add a lateral constraint for the pure shear load case.
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 y direction list, choose Prescribed.
Prescribed Displacement 3
Finally, add another Prescribed Displacement node to prescribe the y displacement in the equibiaxial load case.
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 y direction list, choose Prescribed.
5
In the u0y text field, type stretch-1.
Mesh 1
Mesh each object 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.
Definitions
Define global variables for the nominal stress in the three load cases. These can be expressed as volume averages over each domain.
Average 1 (aveop1)
1
In the Model Builder window, expand the Component 1 (comp1) > Definitions node.
2
Right-click Definitions and choose Nonlocal Couplings > Average.
3
4
In the Settings window for Average, locate the Advanced section.
5
From the Frame list, choose Material  (X, Y, Z).
Average 2 (aveop2)
1
Right-click Average 1 (aveop1) and choose Duplicate.
2
Average 3 (aveop3)
1
Right-click Average 2 (aveop2) and choose Duplicate.
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
Forward Problem
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.
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
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
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 Auxiliary Pressure (comp1.solid.hmm1.pw).
4
In the Settings window for Field, locate the Scaling section.
5
From the Method list, choose Manual.
6
In the Scale text field, type 10[MPa].
7
In the Model Builder window, under Forward Problem > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Displacement Field (comp1.u).
8
In the Settings window for Field, locate the Scaling section.
9
From the Method list, choose Manual.
10
In the Model Builder window, expand the Forward Problem > Solver Configurations > Solution 1 (sol1) > Stationary Solver 1 node, then click Fully Coupled 1.
11
In the Settings window for Fully Coupled, click to expand the Method and Termination section.
12
From the Nonlinear method list, choose Constant (Newton).
13
In the Study toolbar, click  Compute.
Results
Compare the model prediction for the initial guess of the material parameters with the experimental stress–stretch curves.
Initial Model Prediction
1
In the Model Builder window, right-click Experimental 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).
Initial Model Prediction
1
Right-click Initial Model Prediction and choose Global.
2
In the Settings window for Global, type Initial Model Prediction in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
4
Click to expand the Coloring and Style section. From the Color list, choose Cycle (reset).
5
Click to expand the Legends section. Select the Show legends checkbox.
6
From the Legends list, choose Manual.
7
8
In the Initial Model Prediction toolbar, click  Plot.
The initial model parameters yield a prediction of the uniaxial and pure shear data that is off, but of the correct order of magnitude. However, they fail to capture the strain-stiffening behavior in equibiaxial tension.
Component 1 (comp1)
Now, we will set up the parameter estimation problem. Three Least-Squares Objective nodes will be added, one for each load case.
Uniaxial Tension Test
1
In the Physics toolbar, click  Optimization and choose Parameter Estimation.
2
In the Settings window for Least-Squares Objective, type Uniaxial Tension Test in the Label text field.
3
Locate the Experimental Data section. From the Data source list, choose Result table.
The first data column contains the applied stretch, which is the parameter for which the solution needs to be computed.
4
Locate the Data Column Settings section. In the table, enter the following settings:
5
From the Name list, choose stretch (Applied stretch).
6
In the Unit text field, type 1.
The second data column contains the nominal stress measured. These are the values that will be used to evaluate the objective in Equation 2.
7
The Model expression field is used to set the global variable that should be compared with the experimental data. Use the volume-averaged stress variables that we defined when setting up the forward model.
8
In the Model expression text field, type comp1.P_ua.
9
In the Column name text field, type UA.
10
In the Unit text field, type MPa.
Pure Shear Test
Proceed in a similar fashion for the two remaining objectives.
1
In the Parameter Estimation toolbar, click  Least-Squares Objective.
2
In the Settings window for Least-Squares Objective, type Pure Shear Test in the Label text field.
3
Locate the Experimental Data section. From the Data source list, choose Result table.
4
From the Result table list, choose Pure Shear Data.
5
Locate the Data Column Settings section. In the table, enter the following settings:
6
From the Name list, choose stretch (Applied stretch).
7
In the Unit text field, type 1.
8
9
In the Model expression text field, type comp1.P_ps.
10
In the Column name text field, type PS.
11
In the Unit text field, type MPa.
Equibiaxial Tension Test
1
In the Parameter Estimation toolbar, click  Least-Squares Objective.
2
In the Settings window for Least-Squares Objective, type Equibiaxial Tension Test in the Label text field.
3
Locate the Experimental Data section. From the Data source list, choose Result table.
4
From the Result table list, choose Equibiaxial Data.
5
Locate the Data Column Settings section. In the table, enter the following settings:
6
From the Name list, choose stretch (Applied stretch).
7
In the Unit text field, type 1.
8
9
In the Model expression text field, type comp1.P_eb.
10
In the Column name text field, type EB.
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
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 All Least-Squares objectives.
4
Locate the Estimated Parameters section. Click  Add four times.
5
Note that the inequality constraint on the Ogden material parameters can be enforced by setting a lower/upper bound of 0 on each μ,α-pair. With the settings entered here, we will force the first pair of parameters to be positive and the second pair to be negative.
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
Because the forward problem is nonlinear, it is important to activate the continuation solver for the stretch parameter by adding an Auxiliary sweep.
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 stretch (Applied stretch).
5
In the Initial value text field, type 1.
Solution 2 (sol2)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 2 (sol2) node.
3
In the Model Builder window, expand the Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Dependent Variables 1 node, then click Auxiliary Pressure (comp1.solid.hmm1.pw).
4
In the Settings window for Field, locate the Scaling section.
5
From the Method list, choose Manual.
6
In the Scale text field, type 10[MPa].
7
In the Model Builder window, under Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Dependent Variables 1 click Displacement Field (comp1.u).
8
In the Settings window for Field, locate the Scaling section.
9
From the Method list, choose Manual.
10
In the Model Builder window, expand the Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Optimization Solver 1 > Stationary Solver 1 node, then click Fully Coupled 1.
11
In the Settings window for Fully Coupled, locate the Method and Termination section.
12
From the Nonlinear method list, choose Constant (Newton).
Results
Before computing the study, it is often useful to set up a plot that monitors the optimization by comparing the model expression with the experimental data. This type of plot is useful to get visual feedback of the quality of the fit, which can be helpful to detect if the forward model or the solver settings need to be improved.
Parameter Estimation: Ogden Hyperelasticity
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Parameter Estimation: Ogden Hyperelasticity in the Label text field.
3
Locate the Data section. From the Dataset list, choose Parameter Estimation/Solution 2 (sol2).
4
Locate the Title section. From the Title type list, choose Label.
5
Locate the Plot Settings section.
6
Select the x-axis label checkbox. In the associated text field, type Applied stretch (1).
7
Select the y-axis label checkbox. In the associated text field, type Nominal stress (MPa).
Add a second y-axis to be able to visualize all datasets in a single plot. The second axis will be used to plot the much stiffer equibiaxial response.
8
Select the Two y-axes checkbox.
9
Select the Secondary y-axis label checkbox. In the associated text field, type Nominal stress (MPa).
10
Locate the Legend section. From the Position list, choose Upper left.
Initial Model Prediction
The Table Graphs displaying the experimental data can be copied from the plot group created to show the initial model prediction.
In the Model Builder window, expand the Results > Initial Model Prediction node.
Equibiaxial Data, Pure Shear Data, Uniaxial Data
1
In the Model Builder window, under Results > Initial Model Prediction, Ctrl-click to select Uniaxial Data, Pure Shear Data, and Equibiaxial Data.
2
Parameter Estimation: Ogden Hyperelasticity
In the Model Builder window, under Results right-click Parameter Estimation: Ogden Hyperelasticity and choose Paste Multiple Items.
Equibiaxial Data
1
In the Settings window for Table Graph, locate the y-Axis section.
2
Select the Plot on secondary y-axis checkbox.
3
In the Parameter Estimation: Ogden Hyperelasticity toolbar, click  Plot.
Uniaxial Model
1
Right-click Parameter Estimation: Ogden Hyperelasticity and choose Global.
2
In the Settings window for Global, type Uniaxial Model in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
4
Click to expand the Coloring and Style section. From the Color list, choose Cycle (reset).
5
Click to expand the Legends section. From the Legends list, choose Manual.
6
Pure Shear Model
1
Right-click Uniaxial Model and choose Duplicate.
2
In the Settings window for Global, type Pure Shear Model in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
4
Click to expand the Coloring and Style section. From the Color list, choose Cycle.
5
Locate the Legends section. In the table, enter the following settings:
Equibiaxial Model
1
Right-click Pure Shear Model and choose Duplicate.
2
In the Settings window for Global, type Equibiaxial Model in the Label text field.
3
Locate the y-Axis section. Select the Plot on secondary y-axis checkbox.
4
Locate the y-Axis Data section. In the table, enter the following settings:
5
Locate the Legends section. In the table, enter the following settings:
Uniaxial Model
A Filter subnode can be added to hide the part of the model prediction that extends beyond the data range of the experiments.
Filter 1
1
In the Model Builder window, right-click Uniaxial Model and choose Filter.
2
In the Settings window for Filter, locate the Point Selection section.
3
In the Logical expression for inclusion text field, type stretch<=2.
Filter 1
1
In the Model Builder window, right-click Pure Shear Model and choose Filter.
2
In the Settings window for Filter, locate the Point Selection section.
3
In the Logical expression for inclusion text field, type stretch<=1.75.
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
5
Select the Show individual objective values checkbox.
6
Select the Table graph checkbox.
7
In the Study toolbar, click  Compute.
Results
Finally, collect the estimated material parameters in an Evaluation Group. This is available 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/Solution 2 (sol2) > 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
Estimated Parameters
In the Settings window for Evaluation Group, type Estimated Parameters in the Label text field.
The material model is now calibrated and the final values can be used when simulating the behavior of a real-life component.