[sfepy] Questions regarding a linear elasticity problem

Hi,

thank you for SfePy. I would like to use it for checking an analytic implementation of a linear elastic poro-elasticity problem. However, I am a total beginner if it comes to using FEM methods in general, so this email is a request for help, partially probably also on very basic questions. If there's interest I would be happy to write up my journey as a tutorial for absolute beginners.

The problem that I want to solve is the following:

I want to look at a cuboidal inclusion in a stratified elastic half-space. Therefore, we use the usual momentum balance equation:

0 = grad(sigma) + F

and we use stress-strain relations of the kind

sigma_ij = C_ijkl(e_kl)

We can decompose the strain tensor in two parts: e = e_eff - e_I where e_eff is the effective part and e_I is the part stemming from the inclusion.

the stress tensor due to the inclusion is given by

sigma_ij,I = C_ijkl e_kl,I

We want to consider hydrostatic inclusions, therefore the inclusion stress tensor is diagonal: sigma_ij,I = p * delta_ij.

This strain inclusion produces the same deformation as the external force F = -grad p.

We assume that the inclusion is embedded in a stratified half-space: we have a free surface, and sections parallel to the surface with differing linear elastic isotropic material constants.

Now, with regards to solving this problem with SfePy, I have the following questions:

1. I would like to solve the problem on a hexahedral mesh of rectangular cubes. Creating the vertices, and connections is quite straightforward without the need for a mesh generator. However, it is not obvious to me, how to hand this mesh to SfePy:

I assume that I need to use Mesh.from_data(), furthermore, I assume that my vertices go under coors, and my connectivity matrix goes under conns, and that descs should be set to '3_8', but I am not sure about the other parameters.

2. In connection to that, is there a preferred "Dimensionality” for SfePy? I would usually define my grid in the thousands of m, my elastic constants in Pa (which makes for numbers in the order of 1e9), and which results in deformations in the order of cm. Does this cause numerical issues?

3. I had a look through the predfined equations and linear elasticity examples. It appears, as if I should use the

dw_lin_elastic

term for the force balance, but how do I encode the strain inclusion in the weak formalism?

One option would be to use that my force term is given by F=-grad p. However, if I assume a constant pressure change in my inclusion, this gives a pretty exciting delta function like derivative. Is there a better way of doing this?

4. In principle, I would like the top boundary of the box as a free surface, and the other ones should essentially "slide": so not being fixed to u=0. Is there a possibility to accomplish this?

5. Which sort of solver would you recommend for this problem?

6. If I would also like to get the stresses and strains, would you recommend to calculate them afterwards from the displacment or is there a better way of getting them?

7. Obviously, the problem is symmetric around the x = 0 and the y = 0 axis. Is it beneficial to exploit this?

Thanks in advance, Kind regards, Hanno

P.S. in case you are interested, I am trying to check my impementation of https://doi.org/10.1111/j.1365-246X.2007.03487.x : B. Kvshinov, “Reflectivity method for geomechanical equilibria".

Below is what I have so far (mostly cobbled together from various turorials):

import sys import numpy as np import matplotlib.pyplot as plt from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson from sfepy.discrete.fem import Mesh, FEDomain from sfepy.discrete import FieldVariable, Material, Integral, Function,
Equation, Equations, Problem from sfepy.discrete.conditions import Conditions, EssentialBC from sfepy.base.base import IndexedStruct from sfepy.solvers.ls import ScipyDirect from sfepy.solvers.nls import Newton

def create_mesh(x0=-4000., y0=-4000., z0=0., x1=4000, y1=4000, z1=2000, dx=100, dy=100, dz=20): # create the mesh

z_spacing = dz
x_spacing = dx
y_spacing = dy

the_grid = np.mgrid[x0:x1+1:x_spacing,y0:y1+1:y_spacing,z0:z1+1:z_spacing]

xr = the_grid[0].ravel()
yr = the_grid[1].ravel()
zr = the_grid[2].ravel()

vertices = np.vstack((xr, yr, zr)).T
nx, ny, nz = the_grid.shape[1:]

tetrahedra = []
dxs = np.array([ x_spacing, 0., -x_spacing,0.,  x_spacing,0., -x_spacing])
dys = np.array([0., y_spacing, 0, -y_spacing, 0, y_spacing, 0])
dzs = np.array([0, 0,0,z_spacing,0,0,0])
for i, v in enumerate(vertices):
    tet = [i,   i+nz*ny,   i+nz*ny+nz,   i+nz, 
           i+1, i+nz*ny+1, i+nz*ny+nz+1, i+nz+1]
    try:
        cube = vertices[tet]
        dc = cube[1:] - cube[:-1]
        if (np.array_equal(dc[:,0],dxs) & np.array_equal(dc[:,1],dys) & 
            np.array_equal(dc[:,2],dzs)):
            tetrahedra.append(tet)
    except IndexError:
        continue

tetrahedra = np.asarray(tetrahedra)
return vertices, tetrahedra

descs3d = ['3_8']

x0 = -4000. x1 = 4000. dx = 100. y0 = -4000. y1 = 4000. dy = 100. z0 = 0. z1 = 2000. dz = 20.

vertices, tets = create_mesh(x0, y0, z0, x1, y1, z1, dx, dy, dz)

# How do I do this? # mesh3d = Mesh.from_data('mesh', coors=vertices, # ngroups=1, conns=tets,

# Define the regions domain = FEDomain('domain', mesh3d)

# inclusion extent xr0 = -501. xr1 = 501. yr0 = -251. yr1 = 251. zr0 = 1000.-1. zr1 = 1101. rexp = f'vertices in ((x>{xr0})&(x<{xr1})&(y>{yr0})&(y<{yr1})&(z>{zr0})&(z<{zr1}))'

regions = {'Omega' : 'all', 'free_surface' : ('vertices in (z < 0.01)', 'facet'), 'base_boundary': ('vertices in (z > %f.2)'%(z1 - 0.5), 'facet'), 'Gamma_left': ('vertices in (x < %f.1)'%(x0 + 0.5), 'facet'), 'Gamma_right': ('vertices in (x > %f.1)'%(x1 - 0.5), 'facet'), 'Gamma_front': ('vertices in (y < %f.1)'%(y0+0.5), 'facet'), 'Gamma_back': ('vertices in (y > %f.1)'%(x1-0.5), 'facet'), 'reservoir_layer': (f'vertices in ((z > {zr0}) and (z<{zr1}','cell'), 'reservoir' : (rexp, 'cell'), 'overburden' : ('vertices in (z<{zr0})', 'cell'), 'underburden' : ('vertices in (z>{zr1})', 'cell')}

approx_order = 1

fields = {'displacement' : ('real', 3, 'Omega', approx_order)}

# What do I do about the inclusion here?

materials={'overburden' : ({'D':stiffness_from_youngpoisson(3, 1e9, 0.4)}), 'res_layer' : ({'D':stiffness_from_youngpoisson(3, 2e9, 0.25)}), 'underburden' : ({'D':stiffness_from_youngpoisson(3, 5e9, 0.3)})}

variables = { 'u' : ('unknown field', 'displacement', 0), 'v' : ('test field', 'displacement', 'u')}

# Not sure what to do here. bcs = {}

integrals = { 'i': 2, }

## # Force Balance :

equations = { 'elasticity' : """dw_lin_elastic.i.Omega( solid.D, v, u ) = 0""",}