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 poroelasticity 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 halfspace. Therefore, we use the usual momentum balance equation:
0 = grad(sigma) + F
and we use stressstrain 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 halfspace: 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:
 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.
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?
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?
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?
Which sort of solver would you recommend for this problem?
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?
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.1365246X.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)'%(x10.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""",}
participants (1)

Hanno Klemm