[FElupe v.4.0.0] Non-homogenous shear loadcase on a cube with a nearly-incompressible non-affine micro-sphere model

[adtzlr]

import felupe
import matadi

from pypardiso import spsolve

# create a hexahedron-region on a cube
mesh    = felupe.Cube(b=(20, 20, 10), n=(17, 17, 9))
region  = felupe.RegionHexahedron(mesh)

# add displacement, pressure and volume-ratio fields
fields = felupe.FieldsMixed(region)

# apply a non-homogenous shear on the cube
boundaries, dof0, dof1, ext0 = felupe.dof.shear(fields, move=5, top=20)

# define the constitutive material behavior
mat = matadi.MaterialHyperelastic(
    fun=matadi.models.miehe_goektepe_lulei,
    bulk=5000.0,
    mu=0.1475,
    N=3.273,
    p=9.31,
    U=9.94,
    q=0.567,
)

lab = matadi.Lab(mat)
data = lab.run(stretch_max=2.0)
lab.plot(data)

umat = matadi.ThreeFieldVariation(mat)

# newton-rhapson procedure
res = felupe.newtonrhapson(
    fields, 
    umat=umat, 
    dof1=dof1, 
    dof0=dof0, 
    ext0=ext0, 
    # solver=spsolve,
    kwargs={"parallel": True},
    tol=1e-4
)

# save result
felupe.save(region, res.x, filename="result.vtk")

[adtzlr]

And here is another version with incremental prescribed displacements on the boundary move:

import numpy as np
import matplotlib.pyplot as plt

import felupe
import matadi

# create a hexahedron-region on a cube
mesh    = felupe.Cube(b=(20, 20, 10), n=(17, 17, 9))
region  = felupe.RegionHexahedron(mesh)

# add displacement, pressure and volume-ratio fields
fields = felupe.FieldsMixed(region)

# define the constitutive material behavior
mat = matadi.MaterialHyperelastic(
    fun=matadi.models.miehe_goektepe_lulei,
    bulk=5000.0,
    mu=0.1475,
    N=3.273,
    p=9.31,
    U=9.94,
    q=0.567,
)

# plot material behavior on homogenous loadcases
lab = matadi.Lab(mat)
data = lab.run(stretch_max=2.0)
lab.plot(data)

umat = matadi.ThreeFieldVariation(mat)
body = felupe.SolidBody(umat, fields)

# apply a non-homogenous shear on the cube
boundaries, dof0, dof1, ext0 = felupe.dof.shear(fields, top=20)

Move = np.linspace(0, 15, 6)
Force = np.zeros((len(Move), 3))

for increment, move in enumerate(Move):
    
    print("\nStart of increment %d" % (increment + 1))
    print("=====================")
    
    print("\nprescribed values on boundary `move` = %1.4g \n" % move)
    
    boundaries["move"].value = move
    ext0 = felupe.dof.apply(fields, boundaries, dof0)
    
    # newton-rhapson procedure
    res = felupe.newtonrhapson(
        fields, 
        bodies=[body], 
        dof1=dof1, 
        dof0=dof0, 
        ext0=ext0, 
        kwargs={"parallel": True},
        tol=1e-3,
    )
    
    # update fields
    fields -= fields
    fields += res.x
    
    # get reaction force on boundary "move"
    Force[increment] = felupe.tools.force(fields, res.fun, boundaries["move"])
    
    # save result
    felupe.save(region, res.x, filename="result_%d.vtk" % (1 + increment))
    felupe.save(region, res.x, filename="result.vtk")

# plot force-displacement curve
plt.figure()
plt.plot(Move, Force[:,0], "C0o-")
plt.xlabel("Displacement X")
plt.ylabel("Reaction Force X")
plt.title("Non-homogenous shear loading")
plt.grid()

[adtzlr]

The output of the above script gives:

Start of increment 1
=====================

prescribed values on boundary `move` = 0 


Newton-Rhapson solver
=====================

| # |  norm(dx) |
|---|-----------|
| 1 | 1.904e-15 |

Solution converged in 1 iterations within 2.127 seconds.


Start of increment 2
=====================

prescribed values on boundary `move` = 3 


Newton-Rhapson solver
=====================

| # |  norm(dx) |
|---|-----------|
| 1 | 9.127e+01 |
| 2 | 1.203e+00 |
| 3 | 7.767e-03 |
| 4 | 9.237e-07 |

Solution converged in 4 iterations within 9.22 seconds.


Start of increment 3
=====================

prescribed values on boundary `move` = 6 


Newton-Rhapson solver
=====================

| # |  norm(dx) |
|---|-----------|
| 1 | 9.123e+01 |
| 2 | 1.187e+00 |
| 3 | 7.538e-03 |
| 4 | 8.367e-07 |

Solution converged in 4 iterations within 9.241 seconds.


Start of increment 4
=====================

prescribed values on boundary `move` = 9 


Newton-Rhapson solver
=====================

| # |  norm(dx) |
|---|-----------|
| 1 | 9.117e+01 |
| 2 | 1.157e+00 |
| 3 | 7.213e-03 |
| 4 | 7.217e-07 |

Solution converged in 4 iterations within 9.203 seconds.


Start of increment 5
=====================

prescribed values on boundary `move` = 12 


Newton-Rhapson solver
=====================

| # |  norm(dx) |
|---|-----------|
| 1 | 9.114e+01 |
| 2 | 1.109e+00 |
| 3 | 7.191e-03 |
| 4 | 1.435e-06 |

Solution converged in 4 iterations within 9.397 seconds.


Start of increment 6
=====================

prescribed values on boundary `move` = 15 


Newton-Rhapson solver
=====================

| # |  norm(dx) |
|---|-----------|
| 1 | 9.118e+01 |
| 2 | 1.120e+00 |
| 3 | 1.956e-02 |
| 4 | 5.737e-05 |

Solution converged in 4 iterations within 9.16 seconds.

[adtzlr]