Linear Elasticity in 3D


This example provides a demonstration of using PyMKS to compute the linear strain field for a two-phase composite material in 3D, and presents a comparison of the computational efficiency of MKS, when compared with the finite element method. The example first provides information on the boundary conditions, used in MKS. Next, delta microstructures are used to calibrate the first-order influence coefficients. The influence coefficients are then used to compute the strain field for a random microstructure. Lastly, the calibrated influence coefficients are scaled up and are used to compute the strain field for a larger microstructure and compared with results computed using finite element analysis.

Elastostatics Equations and Boundary Conditions

A review of the governing field equations for elastostatics can be found in the Linear Elasticity in 2D example. The same equations are used in the example with the exception that the second lame parameter (shear modulus) \(\mu\) is defined differently in 3D.

\[\mu = \frac{E}{2(1+\nu)}\]

In general, generating the calibration data for the MKS requires boundary conditions that are both periodic and displaced, which are quite unusual boundary conditions. The ideal boundary conditions are given by:

\[u(L, y, z) = u(0, y, z) + L\bar{\varepsilon}_{xx}\]
\[u(0, L, L) = u(0, 0, L) = u(0, L, 0) = u(0, 0, 0) = 0\]
\[u(x, 0, z) = u(x, L, z)\]
\[u(x, y, 0) = u(x, y, L)\]
import numpy as np
from sklearn.pipeline import Pipeline
import dask.array as da

from pymks import (

%matplotlib inline
%load_ext autoreload
%autoreload 2

Modeling with MKS

Calibration Data and Delta Microstructures

The first-order MKS influence coefficients are all that is needed to compute a strain field of a random microstructure, as long as the ratio between the elastic moduli (also known as the contrast) is less than 1.5. If this condition is met, we can expect a mean absolute error of 2% or less, when comparing the MKS results with those computed using finite element methods [1].

Because we are using distinct phases and the contrast is low enough to only need the first order coefficients, delta microstructures and their strain fields are all that we need to calibrate the first-order influence coefficients [2].

The generate_delta function can be used to create the two delta microstructures needed to calibrate the first-order influence coefficients for a two phase microstructure. This function uses the Python module SfePy to compute the strain fields using finite element methods.

x_delta = generate_delta(n_phases=2, shape=(9, 9, 9)).persist()

    x_delta[0, x_delta.shape[1] // 2],
    x_delta[1, x_delta.shape[1] // 2],
    titles=['[0]', '[1]'],

Using delta microstructures for the calibration of the first-order influence coefficients is essentially the same as using a unit impulse response to find the kernel of a system in signal processing. Delta microstructures are composed of only two phases. One phase is located only at the center cell of the microstructure, and the rest made up of the other phase.

Generating Calibration Data

This example models a two-phase microstructure with elastic moduli values of 80 and 120 and Poisson’s ratio values of 0.3 and 0.3, respectively. The macroscopic imposed strain is set to 0.02. All of these parameters used in the simulation are used in the solve_fe function to calculate the elastic strain.

strain_xx = lambda x: solve_fe(
        elastic_modulus=(80, 120),
        poissons_ratio=(0.3, 0.3),

y_delta = strain_xx(x_delta).persist()

Observe the strain field.

    y_delta[0, x_delta.shape[1] // 2, :, :],

Calibrating First Order Influence Coefficients

Calibrate the influence coefficients by creating a model pipeline using the PrimitiveTransformer and the LocalizationRegressor.

model = Pipeline(steps=[
    ('discretize', PrimitiveTransformer(n_state=2, min_=0.0, max_=1.0)),
    ('regressor', LocalizationRegressor())
[13]:, y_delta);

Observe the influence coefficients.

to_real = lambda x: coeff_to_real(x.steps[1][1].coeff).real

coeff = to_real(model)
    coeff[x_delta.shape[1] // 2, :, :, 0],
    coeff[x_delta.shape[1] // 2, :, :, 1],
    titles=['Influence coeff [0]', 'Influence coeff [1]']

The influence coefficients have a Gaussian-like shape.

Predict of the Strain Field for a Random Microstructure

Use the calibrated model to compute the strain field for a random two phase microstructure and compare it with the results from a finite element simulation. The strain_xx helper function is used to generate the strain field.


x_data = da.random.randint(2, size=(1,) + x_delta.shape[1:]).persist()
%time y_data = strain_xx(x_data).persist()
CPU times: user 2min 26s, sys: 38.2 s, total: 3min 4s
Wall time: 27.7 s
    x_data[0, x_delta.shape[1] // 2, :, :],
    y_data[0, x_delta.shape[1] // 2, :, :],

Note that the calibrated influence coefficients can only be used to reproduce the simulation with the same boundary conditions that they were calibrated with.

Now to get the strain field from the model, pass the same microstructure to the predict method.


%time y_predict = model.predict(x_data).persist()
CPU times: user 20.9 ms, sys: 0 ns, total: 20.9 ms
Wall time: 19 ms

Finally, compare the results from finite element simulation and the MKS model.

    y_data[0, x_delta.shape[1] // 2, :, :],
    y_predict[0, x_delta.shape[1] // 2, :, :],
        r'$\mathbf{\varepsilon_{xx}}$ - FE',
        r'$\mathbf{\varepsilon_{xx}}$ - MKS'

Observe the difference between the two plots.

    (y_data - y_predict)[0, x_delta.shape[1] // 2, :, :],
    titles=['FE - MKS']

The MKS model is able to capture the strain field for the random microstructure after being calibrated with delta microstructures.

Resizing the Coefficeints to use on Larger Microstructures

The influence coefficients that were calibrated on a smaller microstructure can be used to predict the strain field on a larger microstructure though spectral interpolation [3], but accuracy of the MKS model drops slightly. To demonstrate how this is done, let’s generate a new larger \(m\) by \(m\) random microstructure and its strain field.

new_shape = tuple(np.array(x_delta.shape[1:]) * 3)
x_large = da.random.randint(2, size=(1,) + new_shape).persist()

The influence coefficients that have already been calibrated need to be resized to match the shape of the new larger microstructure that we want to compute the strain field for. This can be done by passing the shape of the new larger microstructure into the coeff_resize method.


Use the resize coefficients to calculate large strain field. The coefficients can not now be used for the smaller microstructures.


%time y_large = model.predict(x_large).persist()
CPU times: user 35.8 ms, sys: 8.96 ms, total: 44.8 ms
Wall time: 36.7 ms
plot_microstructures(y_large[0, x_delta.shape[1] // 2], titles=[r'$\mathbf{\varepsilon_{xx}}$'])


[1] Binci M., Fullwood D., Kalidindi S.R., A new spectral framework for establishing localization relationships for elastic behav ior of composites and their calibration to finite-element models. Acta Materialia, 2008. 56 (10): p. 2272-2282 doi:10.1016/j.actamat.2008.01.017.

[2] Landi, G., S.R. Niezgoda, S.R. Kalidindi, Multi-scale modeling of elastic response of three-dimensional voxel-based microstructure datasets using novel DFT-based knowledge systems. Acta Materialia, 2009. 58 (7): p. 2716-2725 doi:10.1016/j.actamat.2010.01.007.

[3] Marko, K., Kalidindi S.R., Fullwood D., Computationally efficient database and spectral interpolation for fully plastic Taylor-type crystal plasticity calculations of face-centered cubic polycrystals. International Journal of Plasticity 24 (2008) 1264–1276 doi;10.1016/j.ijplas.2007.12.002.

[ ]: