Does the Neumann boundary condition work in the right way?

PyLith
1

Hi PyLith developers and other users,

I am testing a simple 2D case in which the left side of a rectangle is pressed by a force, the right side is fixed in the horizontal direction, and the bottom side is fixed in the vertical direction. So, the Neumann boundary condition is applied to the left side. However, the output always gives me the solution for the case with an extension applied to the left side, even if compression is exerted on the left side.

In the spatialdb, “traction-shear” is 0, and “traction-normal” is non-zero (which is used to realize either compression or extension).

Did other users encounter a similar issue when using PyLith-2.2.1 and PyLith-2.2.2?

But the 3D case works (7.9.5.3 Step02.cfg Dirichlet and Neumann Boundary Conditions): by changing the sign of normal stress in the spatialdb, we can see the impacts from extension and compression.

Best,
SZ

PS:

The rectangle is 500 km long (x-axis) and 100 km wide (Y-axis). The left boundary is located at x-position of -250.0 km.

Here is the format of spatialdb of Neumann boundary condition

#SPATIAL.ascii 1
SimpleDB {
num-values = 2
value-names = traction-shear traction-normal
value-units = MPa MPa
num-locs = 2 // number of locations
data-dim = 1 // locations on a line
space-dim = 2
cs-data = cartesian {
to-meters = 1.0e+3 // specify coordinates in km
space-dim = 2
}
}
// Columns are
// (1) x coordinate (km)
// (2) y coordinate (km)
// (3) T-shear (MPa)
// (4) T-normal (MPa)
-250.0 -20.0 0.0 20.0
-250.0 -70.0 0.0 20.0

Here is part of the set in *.cfg file.

[pylithapp.timedependent.bc]
boundary_LC_xneg = pylith.bc.Neumann
[pylithapp.timedependent.bc.boundary_LC_xneg]
label = boundary_LC_xneg
db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC on -x
db_initial.iohandler.filename = spatialdb/tractions_pressure_LC.spatialdb
db_initial.query_type = nearest
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
quadrature.cell.quad_order = 2

3

Screen Shot 2023-03-01 at 11.42.09 AM
Screen Shot 2023-03-01 at 11.42.09 AM1254×826 24.1 KB

Let consider compression on the right-side boundary. The rectangle is horizontally 1000-km long and vertically 100-km wide. The elastic properties are similar to the PREM model. We consider viscoelasticity in the lower crust and upper mantle. In case the compression is not large enough (e.g., <2000 MPa), gravitation will make the entire system extensional. In my case, the x-directed displacement on the right side at the first few steps can be as large as several kilometers, greater than the size of the mesh. That x-directed displacement will conitnue to grow. Following the method of adding gravity (run the simulation and get the updated stress state as the input to a new simulation), it’s impossible for me to get an updated stress state that can give me negligible deformation in the new simulation.

4

The stresses due to gravity should be equal to the integral of density*g. For a uniform half space, the stress will increase linearly with depth. For a more complex function of density with depth, the stress will have a higher order variation, either piecewise linear with different slopes in the case of multiple layers with uniform properties, or a quadratic variation for a linear variation in density with depth.

Are you computing your initial stresses to account for the variation in density with depth? Likewise, does your initial normal traction on the boundary exactly match your axial stress from gravity on the boundary?