Stiffness coefficients in granular micromechanics determined with three-point bending experiments and digital image correlation for a bimodulus material

1. Introduction

• Tension/compression asymmetry shifts the neutral axis in bending, providing a macroscopic signature of bimodulus behavior.The shift occurs away from the section centroid when tensile and compressive stiffnesses differ.Source: §1 · p. 2
• Unlike continuum switching models, granular micromechanics represents asymmetry through orientation-dependent normal and tangent interactions.Each orientation is classified by its normal displacement, allowing tension and compression to coexist at one material point.Source: §1 · p. 2
• The study aims to identify simultaneously four independent parameters governing bimodulus interaction response from measured full-field strain data.The analysis is restricted to the classical first-gradient formulation.Source: §1 · p. 3
• Strain-gradient extensions are outside this identification procedure, although they are a possible direction for representing localized effects near the loading region.Such extensions can introduce additional kinematical descriptors and characteristic lengths.Source: §1 · p. 3
• The inverse analysis compares measured and simulated full-field normal strain over the ROI rather than relying only on derived quantities such as neutral-axis position.DIC measurements provide the experimental input.Source: §1 · p. 3
• The optimization couples MATLAB and COMSOL numerical evaluations, using pattern search for exploration followed by quasi-Newton local refinement.The two-stage strategy addresses the lack of directly available gradient information.Source: §1 · pp. 3–4

2.1. Kinematics and energy formulation

• The local strain-energy density aggregates normal and tangent interaction contributions over orientations, with displacements expressed using Green–Saint-Venant strain and a characteristic length.The macroscopic energy is represented quadratically through a fourth-order stiffness tensor.Source: §2.1 · pp. 4–5

2.2. Stiffness identification and tension/compression asymmetry

• The model assigns tensile or compressive normal and tangent stiffnesses according to each interaction’s normal displacement, rather than classifying the whole material point.The orientation integral is evaluated numerically with discrete orientations and updated interaction states [p60?].Source: §2.2 · pp. 5–6
• The first-gradient response identifies stiffnesses conditional on the prescribed characteristic length, because the macroscopic stiffness depends on L^2kη and L^2kτ.Rescaling L requires inverse-square rescaling of interaction stiffnesses; the ratios nη and nτ remain unchanged.Source: §2.2 · p. 6

2.3. Beam setting and neutral-axis description

• In three-point bending, the neutral axis is where G11 vanishes; bimodulus asymmetry shifts it from mid-depth toward the stiffer side.Figure 1 identifies the ROI, comparison cross-section Γ, its offset d from the loading point, and neutral-axis deviation YC.Source: §2.3 · p. 7, Figure 1 · p. 7

3.1. Specimen, three-point bending setup, and loading protocol

• Three-point bending tests on prismatic beams were designed to produce clear strain fields for inverse identification.The specimen is simply supported and loaded at midspan; the campaign used a load–unload–reload protocol with a 3500 N peak selected to keep ROI strains essentially elastic.Source: §3.1 · p. 8, §3.2.1 · p. 9

3.2. DIC principle and experimental setup

• DIC tracks a speckle pattern to provide full-field in-plane displacement and strain measurements, with this study focusing on G11.A single-camera, two-dimensional setup measures in-plane kinematics; Ncorr and PyDIC outputs are cross-checked through a common processing pipeline.Source: §3.2 · pp. 8–9, §3.2.2 · p. 10

3.2.1. Region of interest and loading protocol

• The ROI is centered on the specimen span, spans the beam height, and is used consistently for tests and numerical comparisons.The inverse comparison uses zero-load and peak-load configurations from the second loading phase; tests showed negligible strain-field variation across the investigated loading rates.Source: §3.2.1 · p. 9

3.2.2. DIC outputs and neutral-axis extraction

• At each ROI position, the local neutral axis is located where G11 changes sign, and the resulting coordinates are averaged to obtain YC relative to mid-depth.This compact experimental measure is used in the initial identification stage.Source: §3.2.2 · p. 10

3.3. No-load assessment of DIC accuracy and precision

• No-load accuracy and precision are treated as global indicators based on ROI-averaged G11, not as pointwise estimates of the noise in the identification objective.Four correlations were made against one reference image, so they are repeated correlations rather than fully independent measurements.Source: §3.3 · pp. 10–11

3.4. Material screening and final selection

• Candidate materials were screened for measurable elastic strain and stiffness contrast, but each alternative had a material-specific obstacle to identification.The screened set included gypsum/plaster, 3D-printed polymers, natural stones, cast iron, and NHL 3.5 lime mortar; preliminary checks indicated cracking in the mortar before a sufficiently extended measurable elastic field developed.Source: §3.4 · p. 11
• EG-92E graphite was selected because it provides DIC-suitable strains, remains elastic over the considered loading range, and has sufficient tension/compression contrast.Its nominal compressive modulus sets the material scale; the supplier did not report tensile Young’s modulus, so tensile interaction stiffnesses were initialized from the neutral-axis asymmetry ratio.Source: §3.4 · p. 12

• A physically informed initial parameter set guides the nonlinear regression toward mechanically meaningful values.The measured neutral-axis shift informs the tension/compression asymmetry estimate.Source: §4.1 · pp. 12–13
• The specimen geometry and adopted DIC region of interest are summarized in Table 2.The ROI comparison domain is centered on the specimen span and extends over its full height.Source: Table 2 · p. 13
• The final identification compares numerical and experimental G11 fields over the selected region of interest using a hybrid optimization procedure.Pattern search explores the parameter space before quasi-Newton refinement.Source: §4.1 · p. 13

4.2. Experimental strain field and initial parameter set

• The local neutral axis is extracted where G11 changes sign, and its average position summarizes the measured field.The vertical profile Γ is selected away from the loading point to reduce localized stress-transfer effects in the bending-dominated comparison.Source: §4.2 · p. 13, Figure 3 · p. 15
• Supplier-reported compressive modulus sets only the initial compressive stiffness scale, while the four identified quantities exclude the fixed length L.The prescribed values are L = 10 mm and ν2D = 0.16; the reported stiffnesses correspond to that length scale.Source: §4.2 · pp. 14–15
• The measured neutral-axis deviation estimates the initial bimodulus ratio; for this specimen, YC ≈ −4.20 mm gives n̄ = 1.92.The ratio is assigned equally to the normal and tangent interaction channels in the initial set.Source: §4.2 · p. 14
• Table 3 summarizes the initial and final physical parameter sets, q0 and qopt.The initial values are subsequently corrected through the full-field optimization.Source: Table 3 · p. 15, §4.2 · p. 15

4.2.1. No-load DIC accuracy and precision

• No-load normal strain remains on the order of 10^-5 and serves as a global reference for interpreting measured strain magnitudes.The assessment uses four no-load correlations from five unloaded images.Source: §4.2.1 · p. 16, Table 4 · p. 16

4.3.1. Numerical implementation

• For each trial set, the coupled MATLAB–COMSOL solver updates tension or compression stiffness according to the sign of normal interaction displacement uη.The resulting numerical G11 field is compared with DIC over the same ROI; J is the full-field NRMSE normalized by the observed strain level.Source: §4.3.1 · pp. 16–17
• The numerical bending model applies load over finite 1 mm boundary segments and increases load in five 700 N steps to 3500 N.Five increments were selected after convergence analysis showed negligible stiffness-tensor variation with five or more increments.Source: §4.3.1 · p. 17

4.3.2. Hybrid optimization procedure

• Pattern search explores four scaled parameters without derivatives, then quasi-Newton refinement starts from its best set and continues minimizing the same full-field J.The trial model is rerun and its result accepted or rejected according to whether J decreases.Source: Figure 4 · p. 18, §4.3.2 · pp. 18–19

4.3.3. Evolution of the objective function and identified parameters

• The objective function shows a marked early reduction followed by a smoother decrease as the search approaches a stable parameter region.The plot reports accepted updates; the full identification required 128 workflow evaluations.Source: §4.3.3 · p. 20, §4.3.2 · p. 20, Figure 5 · p. 21
• All four identified parameters change in a coordinated manner as J decreases, with tangent stiffnesses varying more relatively than normal stiffnesses.The plotted parameter histories show accepted states only.Source: §4.3.3 · pp. 20–21, Figure 6 · p. 21

4.4. Comparison of the strain fields before and after optimization

• At Γ, the optimized neutral-axis position is 29.65 mm versus 29.31 mm experimentally; the average values are 29.53 mm and 30.20 mm, respectively.The experimental average was used to initialize the bimodulus ratio, while optimization retained a close prediction over the region.Source: §4.4 · pp. 22–23
• 18.6% to 12.3%: optimization lowers objective J, while improving agreement in the through-depth profile and full-field strain distribution.The largest remaining differences occur near the load-application region.Source: §4.4 · pp. 22–23, Figure 7 · p. 22
• The remaining mismatch is concentrated near the loading region, where concentrated-load effects and local mechanisms may exceed the present first-gradient model’s representation.The study notes strain-gradient formulations as a possible way to represent higher-order kinematics and microstructural regularization.Source: §4.4 · p. 23
• The identified parameter set reproduces the measured bending response with good consistency, supporting full-field G11 data as a basis for inverse stiffness identification.The procedure combines DIC measurements, a granular-micromechanics model, and optimization.Source: §4.4 · p. 23

4.5. Parameter-resolution assessment

• The objective uses the complete two-dimensional G11 field across the ROI, combining strain profiles from multiple longitudinal positions rather than a single section.The analysis perturbs each parameter independently while holding the other three at their identified values.Source: §4.5 · p. 24
• Threshold crossings occur at −12.10% and +15.81% for ekc η, and −16.86% and +19.39% for nη; neither ekc τ nor nτ crosses within ±25%.The threshold is eJth = 1.31 × 10−4, from eJopt = 9.47 × 10−5 plus noise of 3.660 × 10−5.Source: §4.5 · pp. 24–25
• The normal-interaction parameters ekc η and nη are better resolved because the measured G11 field is less sensitive to the tangent-interaction parameters ekc τ and nτ.G22 was also measured, but its smaller magnitude and comparable measurement noise gave it insufficient signal-to-noise for this identification.Source: §4.5 · p. 25
• The four parameters should be interpreted collectively as a set reproducing the measured full-field bending response, not as four independently and equally resolved material constants.More accurate G22 measurements with improved signal-to-noise could provide complementary information about the tangent-interaction parameters.Source: §4.5 · pp. 25–26
• For hi = 0.01, finite-difference response changes for ekc η and nη exceed the 3.66 × 10−5 noise level, while those for ekc τ and nτ fall below it.The increment is retained because gradients are evaluated on the complete deterministic objective; it is not treated as an experimentally resolved parameter variation.Source: §4.5 · p. 26

5. Conclusions

• The identified parameter set reduces discrepancy with experimental strain data, improving both the through-depth profile along Γ and the spatial distribution of G11 over the ROI.The procedure combines DIC, forward numerical simulations, and hybrid optimization.Source: §5 · p. 26
• The four parameters are not equally resolved: normal-interaction parameters are more strongly constrained than tangent-interaction parameters by this bending experiment.Source: §5 · p. 26
• Only G11 was used for identification because the smaller G22 and G12 components were more sensitive to measurement noise and insufficiently robust.Source: §5 · pp. 26–27
• The calibration applies to the elastic regime: the experiment sought to remain elastic, and the forward model excludes plasticity and other inelastic effects.Source: §5 · p. 27
• Residual differences are concentrated near the load-application region, where localized load transfer and idealized contact are not fully captured by the first-gradient model.Source: §5 · p. 27
• Reported stiffness coefficients are conditional on prescribed L; changing L requires a 1/L^2 rescaling, while a strain-gradient model could identify L with the four bimodulus parameters.Source: §5 · p. 27
• Post-processing could relate the identified granular stiffnesses to equivalent tensile and compressive Young’s moduli and Poisson’s ratios under idealized uninormal loading.Source: §5 · p. 27
• The authors frame the procedure as a first step, with richer kinematics and improved local modeling among the proposed extensions.Source: §5 · p. 27