Antiplane and plane strain kernels for slip and force

[mallickrishg]

Tasks for @mallickrishg and @brendanjmeade

• [ ] Derive quadratic slip greens functions in terms of the 3 quadratics we already use, instead of $a_0 + a_1x+a_2x^2$
• [x] Think about how to absorb all the antiplane calculations into existing bemcs functions like these. The options are
• rewrite the existing bemcs functions for antiplane and rename everything with suffix 'ap' (antiplane) or 'ps' (plane strain)
• add specific flags/variables within els to specify what the geometry is, and then modify existing functions. The latter requires more planning but fewer lines of new code.
• @brendanjmeade your advice is much needed and appreciated!
• [ ] Incorporate all the new kernels for linearly varying force on line segments into bemcs.
• Note that there is subtle difference in the kernels we want compared to the ones we have.
• In antiplane, we have $uz,\epsilon{zx},\epsilon_{zy}$ in response to a force $f_z$ where the body force comes from $\nabla^2 u_z = f_z$.
• This is different from the plane strain kernels $\left(u_x,uy,\sigma{xx},\sigma{yy},\sigma{xy}\right)$ where the body force is $\left(f_x/\mu,f_y/\mu\right)$.
• Also, for consistency I think we should derive $\epsilon{xx},\epsilon{yy},\epsilon{xy}$ kernels - because strains are the actual terms needed for the heterogeneous material calculations, not stress.

[brendanjmeade]

• rewrite the existing bemcs functions for antiplane and rename everything with suffix 'ap' (antiplane) or 'ps' (plane strain)
• add specific flags/variables within els to specify what the geometry is, and then modify existing functions. The latter requires more planning but fewer lines of new code.

Hard to say which is right before building, but I'd lean towards option 2 if the arguments passed and variables returned are the same. That would make it one-stop shopping.

• Also, for consistency I think we should derive ϵxx,ϵyy,ϵxy kernels - because strains are the actual terms needed for the heterogeneous material calculations, not stress.

The way I've previously done this is by converting stress to strain, with the lines starting at:

https://github.com/brendanjmeade/bemcs/blob/7ec2b95615cff1023985a3640f4218c797969f81/bemcs/bemcs.py#L2134

[mallickrishg]

• rewrite the existing bemcs functions for antiplane and rename everything with suffix 'ap' (antiplane) or 'ps' (plane strain)
• add specific flags/variables within els to specify what the geometry is, and then modify existing functions. The latter requires more planning but fewer lines of new code.

Hard to say which is right before building, but I'd lean towards option 2 if the arguments passed and variables returned are the same. That would make it one-stop shopping.

That would be perfect. Anyway this is something I will keep in mind as I add the functions

The way I've previously done this is by converting stress to strain, with the lines starting at:

https://github.com/brendanjmeade/bemcs/blob/7ec2b95615cff1023985a3640f4218c797969f81/bemcs/bemcs.py#L2134

aah nice! that is pretty convenient! So i guess we can still write all the functions in terms of $u,\sigma$, and then just use the matrix operator to convert from stress to strain.

[mallickrishg]

tasks for @mallickrishg

• [x] Derive quadratic slip kernels for antiplane using same basis functions as in plane strain case, and push code
• [x] Push code for linear force kernels for antiplane
• [x] Write kernel calculation code for slip, displacement and stress Greens functions, with appropriate rotations for antiplane
• [x] Write example notebooks showing the use of slip and force kernels (antiplane)
• [ ] Push code for linear force kernels for plane strain (low priority for now)

tasks for @brendanjmeade

• [ ] Create function that carries out numerical integration of constant body force kernels over arbitrary triangles: scipy.dblquad() integration for interior point evaluation, quadpy for external evaluation. This will require some arbitrary rotation for dblquad(), but not necessary for quadpy()
• [ ] Write example notebook showing how to carry out the kernel computation

[mallickrishg]

@brendanjmeade i realized something really odd about antiplane. Dealing with singularities due to slip on non-planar faults in antiplane is more severe than in plane strain! This is purely visual, but when I try similar experiments for non-planar faults as we tried for plane strain, even adding in the smoothness conditions to quadratic slip does not make singularities go away (see attached image for comparison of antiplane vs planestrain for a 3 segment fault with a bend).

Antiplane stress components show a pronounced singularity at the bend

Plane stress components do not show a singularity at the bend!

[mallickrishg]

One of the big differences between antiplane and plane strain is of course slip is a scalar (antiplane) vs a vector (plane strain).

So, for non-planar faults the gradient continuity condition in antiplane strain ends up being independent of geometry i.e., $\frac{ds}{d\zeta}$ is only a function of fault segment length where $\zeta$ is the along-fault direction - it has no information about its orientation in space.

In contrast, the plane strain gradient continuity conditions requires continuity of $\frac{ds_x}{d\zeta},\frac{ds_y}{d\zeta}$, which is not perfect because we actually need continuity of the Jacobian $\nabla s$ but for small changes in fault angles, it seems to work just fine.

The presence of these singularities in antiplane has sort of stumped me. Thoughts welcome (and needed).

[mallickrishg]

To be clear this issue does not exist for planar faults

[brendanjmeade]

@mallickrishg Thanks for sharing this. Certainly didn't expect this. Just to check a basic point, the Burgers vector varies smoothly for these anti-plane strain experiments? I don't have an immediate technical response, other than to say that this seems somewhat in contrast to the Hess (1973) result that we sort of used as support for the plane strain case.

[brendanjmeade]

Also, are you sure there aren't singularities in the planar case too? Visually, there seems to be smoothing near +- 0.5. That might be the black dots denoting element boundaries. Does this still look smooth when visualizing the log of the computed stresses?

[mallickrishg]

@mallickrishg Thanks for sharing this. Certainly didn't expect this. Just to check a basic point, the Burgers vector varies smoothly for these anti-plane strain experiments? I don't have an immediate technical response, other than to say that this seems somewhat in contrast to the Hess (1973) result that we sort of used as support for the plane strain case.

Yes, the burgers vector (scalar) various smoothly. The point I want to highlight though is that: the actual values of slip along the fault remains the same irrespective of the relative angles of the 2 fault segments, as long as length of fault 1 = fault 2. This is because the slip continuity and slip gradient continuity equations contain no information about the relative angles of the 2 fault segments.

So, it makes sense that we have a singularity with nonplanar faults since the slip distribution perfectly cancels singularities only for an equivalent planar fault case. I can't think of a way to actually incorporate fault angles into these calculations.

[mallickrishg]

Also, are you sure there aren't singularities in the planar case too? Visually, there seems to be smoothing near +- 0.5. That might be the black dots denoting element boundaries. Does this still look smooth when visualizing the log of the computed stresses?

The black spot is the element boundary. When I plot stress components evaluated exactly on the fault, I see it is continuous. Haven't plotted log.

[brendanjmeade]

I can't think of a way to actually incorporate fault angles into these calculations.

I think I'm understanding this better now. I don't have any quick hit ideas. I would ask this question: Is there an application for the anti-plane strain case that really needs this? Most everything we've recently discussed is for the plate strain case.

[mallickrishg]

I would ask this question: Is there an application for the anti-plane strain case that really needs this? Most everything we've recently discussed is for the plate strain case.

Thankfully, I can't think of any serious applications that need this. Just thought to highlight this issue as a real surprise!

[brendanjmeade]

I think you've found something deeply interesting that I certainly had not anticipated!

[mallickrishg]

Happy Friday! I have a working example for vertically stratified elastic structure. I currently implemented it for constant force kernels only, I will work on the linear force kernels next.

[brendanjmeade]

This notebook makes for more than a happy Friday!

Amazing notebook that's clear and super easy to modify too. For fun, I did an example with a 100x increase in shear modulus at the material properties at the 3rd (non-topo) layer (see attached). Honestly, I'm pretty amazed that this hasn't been done decades ago in geophysics. This just opens the door to so many fun and meaningful calculations!

I'm done with classes in a few weeks and am so looking forward to getting triangle integration done for the antiplane strain case.

[mallickrishg]

More improvements in the notebook (please pull the most recent repo). I have now provided proper documentation describing the method, and how we implement it. Added surface displacement plotting and comparison with a homogeneous half-space.

Now I need to think about how to deal with the linear force kernels - implementing the spatial awareness part (especially for triple junctions) is turning out to be a lot more tricky than the quadratic slip kernels.

[mallickrishg]

@brendanjmeade - i have been thinking about how to deal with triple junctions using linear basis functions for forces along line segments and I have come across a surprising problem. For a system with $N_t$ triple junctions, I find that we fall short of the number of equations needed to make the system matrix square by $\frac{N_t}{2}$ (specifically for meshes that have no open nodes)!

The boundary conditions for nodes that I apply are

• continuity of force for overlapping elements
• continuity of force at triple junctions

What surprises me about this problem is that the deficit in equations is $\frac{N_t}{2}$ and not $N_t$! Let me know if you have time this week to zoom and chat about this.

[brendanjmeade]

I'm sure you've done more algebra here than I have. I know I've been answering like this a lot lately, but can we simply punt on doing this problem because we're probably going to use triangles for most real-world problems?

[mallickrishg]

I think you are right, there may not be a great case for triple junctions with forces on line segments. Anyway we already have the traditional constant force on a line BEM working, the triangle forces are certainly the next big step.