Describing smooth slip or eigen strains in antiplane geometry

[mallickrishg]

I was thinking about this today for a variety of reasons, but the most tantalising one for me is to attack the poroelasticity problem. I know you have mentioned working on the faulting part of the problem for glaciology - which also sounds really interesting!

We already have smooth slip on faults in an elastic medium, and so we can evaluate the full stress tensor anywhere in the medium. The poroelastic (and similarly the thermoelastic problem) requires a scalar source within the volume whose spatial gradients drive flow, and these pressure gradients also cause elastic deformation of the bulk. Since the pressure field is a scalar, I think it's greens functions can be taken from solutions to Laplace's equation, same as for faulting in anti-plane. In Banerjee 1994 (updated BEM textbook), he has a section where he discusses the use of particular integrals to solve this problem as a boundary value problem. Unfortunately even after reading that chapter, I don't actually understand how it works.

However, we have to integrate these sources over triangles - which means we need to think about how to deal with quadratic basis functions that vary in $x,y$. I don't have any answers right now, just thought I will open this issue so that we can think about it in the coming days.

[mallickrishg]

@brendanjmeade - some half baked thoughts for the weekend

I was looking into what the appropriate boundary conditions are for dealing with non-localised deformation sources, with the big-picture goal being to figure out how to incorporate body forces in the volume.

I started with something small: a rectangular box source in antiplane strain. Valere & Sylvain Barbot derived these solutions previously, but they have an ad-hoc way of dealing with stresses. Basically they integrate the point source (dislocation type) solution to Laplace's equation over a rectangle to get a spatially continuous displacement field. But then the stress field becomes $\sigma_{xz} = \frac{\mu}{2}\frac{\partial uz}{\partial x} - \mu\epsilon{xz}^i\Pi(x,z)$ where the second term is a constant source term defined ONLY within the rectangle domain. This makes $\sigma$ discontinuous! My interpretation is that it seems like they were basically trying to mimic integrated force/traction solutions, which naturally produce discontinuous stress fields but continuous displacement fields in contrast to slip-type solutions!

So i tried to use antiplane constant slip elements to reconstruct the Lambert&Barbot solution. The boundary conditions I chose were:

• left edge of the rectangle: $u_z = -1$
• right edge of rectangle: $u_z = 1$
• top & bottom: traction-free

The BEM solution is similar to the Lambert&Barbot solution, but not identical (attached figure). And you can see that the BEM solve predicts opposite sign displacements in the interior of the rectangle than what one would intuitively think is going on inside a 'shear zone' that accommodates some strain.

All this to say that it might be a good idea to attempt both antiplane & planestrain volumetric sources using the smooth-slip formulation to construct elementary geometries, like rectangles or triangles, instead of pursuing the task of numerical integration of polynomial functions. Thoughts?

For quadratics over something like a rectangle, we could just impose the describe 3 quadratic basis functions as antisymmetric spatially varying displacement boundary conditions on the lateral edges, and set the top and bottom boundaries to traction-free. OR we need to formulate this in terms of traction elements?

I need to think a lot more about how to deal with the next big step i.e., how to describe smoothness over a common edge shared by two volumetric elements.......

[brendanjmeade]

This is really interesting. I definitely don't understand it well at this point.

Let me ask a bit of a zoom-out question. If we want to include body forces, we can return to the Kelvin kernel. If we can integrate this with some shape functions over a single right triangle, shouldn't that give us exactly what we need? That still seems like a meaningful path to me. While we haven't done the integration yet the singularity is relatively weak and adding a shape function might even reduce the order of the singularity? Just sharing some quick thoughts. Don't stop with the good energy you have going on this!

[mallickrishg]

If we want to include body forces, we can return to the Kelvin kernel. If we can integrate this with some shape functions over a single right triangle, shouldn't that give us exactly what we need? That still seems like a meaningful path to me.

Absolutely! That would be the way to do it. What I am suggesting here is a hack that I think (and hope) gives us the same answer with our existing tools, without having to deal with Mathematica (potentially).

[brendanjmeade]

I get that and this could be a really interesting approach. I think I would like to get back to working directly with the Kelvin kernels but I want to get the first paper out the door first.

[mallickrishg]

Sound strategy, can't argue with that. It's just that I am trying to understand how the volumetric source eigen strain solutions by Barbot et al., 2017 are related to the Kelvin source solutions. I'm sure we will eventually figure all this out.

[mallickrishg]

After this, I will get back to contributing to the manuscript.

I just realised something about the point force integration over a volume that might be really interesting! I think the singularities totally go away when we move from a line to a rectangle source. So perhaps there is no need to even seek out quadratic basis functions? (some early celebration 👯 )

The singularity in stress components, at end points, appears because of $\left(\log\left(x^2 + (a-y)^2\right)-\log\left(x^2 + (a+y)^2\right)\right) \to \infty$ as $y\to \pm a$, type of terms when the source is a line. However, when the source is rectangular, I think these terms cancel out leaving the stress components bounded. Attached are rectangle sources in antiplane & planestrain numerically integrated using integral2 in MATLAB. I have verified that the $\log$ singularities do indeed cancel out along the edges and at corners for the antiplane problem using Mathematica.

Attached below are 2 figures - the first one is for antiplane while the second one is for plane strain (force components $(f_x,f_y) = (1,0)$).

[brendanjmeade]

Whoa!!! I don't understand it but I like it. Do you think something similar would happen over triangles as well?

[mallickrishg]

I'm glad you asked! It works just the same for triangles too, at least numerically. Below is $\sigma_{xx}$ when $(f_x,f_y) = (1,0)$ for a right triangle.

[brendanjmeade]

Now this is incredibly good news. I need to understand this better and I also wonder why we didn't get this before. Is it because integral2 is numerical that we're getting these solutions now whereas before we/I was still trying to do this with symbolic integration tools?

[mallickrishg]

@brendanjmeade That's my feeling too. Plotting things and looking at them helps solve problems???? My new/old/new moto to live by.

The only reason I realized this now was because I started playing with the anti-plane problem, which I was doing by hand and realized that the singularities seem to go away for finite dimensional volume sources. And it seems like that might also be the case for the Kelvin problem!

[brendanjmeade]

I definitely see that logic. Anticline has a singularity that's pretty similar to the Kelvin kernels so the same integration technique should work for both. And log is not a super bad singularity so I'm willing to work with a numerical solution here. So nice!

[mallickrishg]

If you're interested, you can play around with the integration with this Mathematica notebook.

https://www.dropbox.com/scl/fi/wpt5rw0zzoc1vedg7kp4o/kelvinproblem_stress.nb?rlkey=6j9vp9ghyfbekb3w0aqc4d9rf&dl=0

[brendanjmeade]

Thanks for sharing the link. I definitely want to play with this but I have to have beginning of the semester discipline right now!

[mallickrishg]

@brendanjmeade Good news - Looks like any standard quadrature technique with sufficiently high precision can deal with the 2-d Kelvin-type integrals quite well. So we may not really need MATLAB, or we don't pay a serious penalty for not using integral2.

Here is an example of integrating a $\log(r)$ kernel (for a body force in antiplane) and plotting the stress components using integral2 vs Gauss-Legendre (61x61 points) and tanh-sinh quadrature (400x400 points). The volume source is present in the domain $\left[-1\leq x \leq 1, -1\leq y \leq 1\right]$

[brendanjmeade]

Thanks Rishav. My visual inspection of these suggests that integral2 is way better. No? Both of the others are noisy and that's not good. I've requested a pdf of the original integral2 paper from my library. Hoping there is code in there.

[mallickrishg]

you are right, that the standard quadrature methods are noisier. But I assume there is no real penalty for using them since the integrals itself are non-singular, and the coincident evaluation for the stress components is just 0. But if you can get integral2 from your library - that's even better.

[mallickrishg]

I've been thinking about how to relate these body force solutions to either poroelastic or viscoelastic strains.

These equivalent body forces, for the antiplane problem at least, can be related to an application of an arbitrary viscous strain source $\epsilon^v{1x},\epsilon^v{1y}$ as $f = -\nabla.(C:\epsilon)$, where $C$ is the constitutive relation tensor and the $1$ or $z$ is the antiplane displacement direction (taken from Segall 2010, Barbot&Fialko, 2010). For an elastic response, the above expression reduces to: $f= -2\mu\left(\frac{\partial\epsilon^v{1x}}{\partial x} + \frac{\partial\epsilon^v{1y}}{\partial y}\right)$

That means the solutions plotted in the previous post correspond to an equivalent body force resulting from linearly varying viscous strains i.e., $\epsilon^v_{1x} = a_0 + a1x$ OR $\epsilon^v{1y} = b_0 + b_1y$. This is sort of where I am stuck conceptually - how can a spatially varying strain source not produce singular stresses? I get that the resulting force is essentially a boxcar function in space, but the underlying viscous strain is linearly varying. Say at one edge of the rectangle, strain is 0 and so is continuous with the rest of the medium. Then at the other end it is non-zero but then immediately outside the rectangle domain, viscous strain is 0. I need to think about this more, a lot more.

I hope that the next time we schedule a chat to discuss the body force problem, I would have figured out these details.

[mallickrishg]

YESSS!!! I think I figured out how to deal with viscous/pore-pressure strain sources as equivalent body forces! We need to schedule a time, when you are properly free, for me to show you figures and convince you of how the Banerjee & Butterfield formulation can be used for viscoelastic or poroelastic problems with arbitrary spatial variations in viscosity or permeability respectively.

Below I lay out the framework for how we can use the fictitious force method, using line sources and volume sources, to do the relevant problem. This is more for me, so that in a couple months or so when we work on this problem I will have notes to work off.

As i mentioned in my previous post, in antiplane the equivalent body force can be written as: $$f = -\nabla . \left(C:\epsilon^v\right)$$ for a constitutive relation tensor $C$, and eigen strain tensor $\epsilon^v$. Note that the total strain tensor is $$\epsilon = \epsilon^v + \epsilon^e$$ where the second term is the elastic response. Since we only care about the elastic response to this imposed eigen strain tensor, the body force reduces to: $$fz = -2\mu\left( \frac{\partial \epsilon^v{zx}}{\partial x} + \frac{\partial \epsilon^v{zy}}{\partial y} \right)$$ where $\epsilon^v{zx},\epsilon^v{zy}$ are the only non-zero components of strain in this geometry. For the following discussion, let us only focus on $\epsilon^v{zx}$ as a source.

Let the eigen strain source be of constant strength, $\alpha_0$, defined over a rectangular domain $-w_x\leq x \leq w_x,-w_y \leq y \leq wy$. We can write this source as: $$\epsilon^v{zx} = \alpha_0 \Pi\left(\frac{x}{w_x}\right)\Pi\left(\frac{y}{w_y}\right)$$ where $\Pi$ refers to box-car functions. The equivalent body force is: $$f_z = -2\mu\alpha_0 \Pi\left(\frac{y}{w_y}\right) \left[ \delta\left(\frac{x}{w_x}-1\right) - \delta\left(\frac{x}{w_x}+1\right) \right]$$ The $\delta$ functions are applied exactly at the edges of the eigen source domain (in the $x$-direction). This means that we can simply use two equal and opposite fictitious force sources defined over a line of length $2w_y$ placed at $x=\pm w_x$ to model the equivalent body force. In fact, that is exactly what I did in the attached figure. I compared them to solutions from Lambert & Barbot 2016 who integrated a displacement greens function $u_z = \frac{\cos\theta}{r}$ over a 2-d domain.

Of course, the problem with singularities still exists because our equivalent body forces are $\delta$ functions! However, we can remedy this the same way we did with bemcs. More in next post.

[mallickrishg]

Say we have a spatially varying eigen strain source of the form, $$\epsilon^v_{zx} = \left(\alpha_0 + \alpha_1 x\right) \Pi\left(\frac{x}{w_x}\right)\Pi\left(\frac{y}{wy}\right)$$ the equivalent body force requires differentiating $\epsilon^v{zx}$ w.r.t $x$, $$f_z = -2\mu\Pi\left(\frac{y}{w_y}\right)\left(\left(\alpha_0+\alpha_1 x\right) \left[ \delta\left(\frac{x}{w_x}-1\right) - \delta\left(\frac{x}{w_x}+1\right) \right] + \alpha_1\Pi\left(\frac{x}{w_x}\right)\right)$$

Evaluating the $\delta$ functions at $x = \pm w_x$ in $f_z$ gives, $$f_z = -2\mu\Pi\left(\frac{y}{w_y}\right)\left(\left(\alpha_0 \pm \alpha_1 w_x\right) \left[ \delta\left(\frac{x}{w_x}-1\right) - \delta\left(\frac{x}{w_x}+1\right) \right] + \alpha_1\Pi\left(\frac{x}{w_x}\right)\right)$$ The first term is the same as before - line forces with strengths $-\left(\alpha_0 - \alpha_1w_x\right)$ at $x=-w_x$ and $\alpha_0 + \alpha_1w_x$ at $x=w_x$, while the second term comes from integrating a point force of strength $\alpha_1$ uniformly over the rectangular domain. The stress singularities still exist because of the line forces, but DO NOT exist for the volume integral!

So if we want to now deal with a meshed volume/bulk (with rectangles), we simply need to ensure that the strengths of the line sources $\pm \left(\alpha_0 \pm \alpha_1 w_x\right)$ at the common edge shared by two rectangles are equal. That cancels the singularities and makes stresses continuous everywhere!

PROBLEM SOLVED!

[brendanjmeade]

Whoa! This is seriously cool, and I'm very grateful for the details explain. That's helping me a lot and allows me to dwell on the core ideas.

I have a lot to understand here, so let me ask a few basic questions:

1. Does "eigenstrain" just mean deformation associated with a deformation source interior to some boundary?
2. The idea that the singularity effect decreases with integration dimension makes sense.
3. I can't yet tell how similar this is compared with Butterfield and Banerjee (at least what I recall). The core idea there is that all non-elastic terms (linear or non-linear) can be lumped into effective body force terms. In the formulation you've presented above, there might(?) an assumption that there is linear relationship between strain and force? Not saying that's bad just trying to understand how these approaches might be similar or different.
4. Can this be done with triangles too? I'm asking because I'm interested in things that are conformal with complex geometries.

I'm looking forward to learning more and pushing this forward

Can you Zoom next Thursday (September 28th)? This is super exciting!

[mallickrishg]

@brendanjmeade the answer to all your questions is yes.

3. I can't yet tell how similar this is compared with Butterfield and Banerjee (at least what I recall). The core idea there is that all non-elastic terms (linear or non-linear) can be lumped into effective body force terms. In the formulation you've presented above, there might(?) an assumption that there is linear relationship between strain and force? Not saying that's bad just trying to understand how these approaches might be similar or different.

It is exactly the same idea! eigen strain is just a way to parameterize a unit strength source term and getting a linear elastic response. To solve a poro-elastic or viscoelastic system, with whatever non-linearities or arbitrary spatial heterogeneities in non-elastic material properties (permeability for P.E and viscosity for V.E), we will need to provide a second set of equations for potential flow or viscous flow to couple with the elastic problem.

4. Can this be done with triangles too? I'm asking because I'm interested in things that are conformal with complex geometries.

Actually I don't know the answer to question 4(if the approach works for triangles), I need to create a simple example and test it.

I have been exploring a lot of this stuff here: https://github.com/mallickrishg/testing_greensfunctions Unfortunately, much of the work is in antiplane, and the code is in MATLAB. However, the fundamental idea of constructing greens functions using line sources + volume sources is extendible to plane strain readily, and we will of course need to write it in python.

Thursday at 12.30pm Boston time?

[mallickrishg]

Attached is a figure for displacement field and then stress components for 3 rectangular strain sources chosen such that they are spatially continuous eigen strain values. The resulting deformation fields are all smooth!

The displacement transect is taken at $y=0$

[brendanjmeade]

It is exactly the same idea! eigen strain is just a way to parameterize a unit strength source term and getting a linear elastic response. To solve a poro-elastic or viscoelastic system, with whatever non-linearities or arbitrary spatial heterogeneities in non-elastic material properties (permeability for P.E and viscosity for V.E), we will need to provide a second set of equations for potential flow or viscous flow to couple with the elastic problem.

I'm slowly getting my head around this. With the parameterization that you described, it still looks like there's the assumption that there's a linear relationship between strain and force. That's the $f=\nabla (C:\epsilon)$ term. To me, that seems like a difference from the idea of simply letting the force be anything and then using the Kelvin kernels. Thoughts?

[mallickrishg]

You are spot on about the linear relationship between strain and force! This is what took me the longest time to wrap my head around.

However, that relationship is between the applied non-elastic component of the strain tensor and the equivalent body force. The assumption is that the total strain tensor is a linear sum of the elastic component and the non-elastic component. So, as long as the total strain tensor has commutative properties, we lump all the non-linearities into the non-elastic term and the Kelvin kernels allow us to calculate the elastic response from the non-elastic terms.

[mallickrishg]

Just adding the derivation/discussion from yesterday here, and extending it for triangles:

Starting with an anti-plane geometry, the equilibrium equation is

$$ \nabla.\sigma = 0$$

where $\sigma = \left[ \sigma_x,\sigma_y \right]$ are the only two non-zero components of stress.

Assuming linear elasticity, we know $\sigma = 2\mu\epsilon_e$ where the elastic part of the strain tensor can be decomposed as $\epsilon_e = \epsilon - \epsilon^v$. The total strain tensor $\epsilon = \left[\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}\right]$. Substituting this into the equilibrium equation, we get:

$$\nabla^2 u - \left(\frac{\partial\epsilon^v_x}{\partial x} +\frac{\partial\epsilon^v_y}{\partial y}\right) = 0$$

where the equivalent body force $f = \left(\frac{\partial\epsilon^v_x}{\partial x} +\frac{\partial\epsilon^v_y}{\partial y}\right)$. Typically we start with solutions for $u(x,y)$ in terms of a point source $\delta(0,0)$ which gives us Greens functions of the form

$$G(x,y,0,0) = \frac{1}{4\pi}\log\left(x^2 + y^2\right)$$

The elastic response to any spatially variable body force source can be calculated by integrating $G(x,y,\zeta,\xi)$ over the domain $\zeta,\xi$ with the appropriate spatially variable strength.

For the case of rectangles that are of constant strength, these sources have been called eigen strains OR constant strain sources - I think we agreed that those are incorrectly named and have been the source of our confusion for a long time now, since these sources do not result in constant strain or stress, and it also does not represent a constant strength body force (as I show below). They do however involve a constant strength for the non-elastic component of the strain tensor over a given domain. For now I just call it 'constant strength' as opposed to linearly varying strength, but we need a more correct nomenclature.

For the constant strength case, let us consider $\epsilon^v = \left[\alpha_0,0\right]$ which can then be defined as:

$$\epsilon^v_x = \alpha_0 \Pi\left(\frac{x}{w_x}\right)\Pi\left(\frac{y}{w_y}\right), \epsilon^v_y = 0$$

The integration of the Greens function needs to then be done for a body force

$$f = \frac{\partial\epsilon^v_x}{\partial x} = \alpha_0 \Pi\left(\frac{y}{w_y}\right) \left[ \delta\left(\frac{x}{w_x}-1\right) - \delta\left(\frac{x}{w_x}+1\right) \right]$$

which reduces to applying line forces of equal and opposite strengths, which represent the $\delta$ force functions, at $x = w_x$ and $x = -w_x$ with length $w_y$! There is no need for a volume integral - only the boundaries matter.

Similarly for a triangle, if we have 3 edges of lengths $a_1,a_2,a_3$ with unit normal vectors $\hat{n}_1,\hat{n}_2,\hat{n}_2$, we can perfectly capture a 'constant strength' source by applying $\delta$ force functions with exact same lengths and orientations as the three edges with strengths given by

• the dot product $\hat{n}.[1,0]$$ when we consider $\frac{\partial}{\partial x}$ terms and
• the dot product $\hat{n}.[0,1]$ when we consider $\frac{\partial}{\partial y}$ terms.

[mallickrishg]

Below is a figure for the still incorrectly names 'constant strength' source within a triangle, showing the solutions from Barbot 2018 is exactly the same as putting line forces along the triangle edges with appropriate amplitudes. What I need to figure out is how to deal with spatially variable strength within the triangle? The nice thing is that we already have solutions for the line forces and volume integrals in the triangle, we just need to figure out how to combine them into a correct equivalnet body force that represents linearly varying non-elastic strain within a triangle.

[brendanjmeade]

This helps so much and I love documenting this. I think I'm stuck on the first step: Going from the divergence of the stress tensor to the Poisson equation. Here is an attempt to do this algebraically. I'm showing this as an image because I"m having trouble getting some of the longer equations to render in GitHub (must be the rain!)

This is close, but it's not quite the Laplacian. The second term is off because it would need to be $\partial_y \partial_y u_y$ instead of $\partial_x \partial_y u_y$. This is what I'm stuck on.

[mallickrishg]

You won't get the Laplacian if you are in plane strain, the solutions I laid out were for anti-plane only. Where $\nabla.\sigma = \frac{\partial\sigma{zx}}{\partial x} + \frac{\partial\sigma{zy}}{\partial y} = 0$.

And since $\sigma{zx} = 2\mu\epsilon{zx}^e$ and $\sigma{zy} = 2\mu\epsilon{zy}^e$, we get

$$\frac{\partial\epsilon^e{zx}}{\partial x} + \frac{\partial\epsilon^e{zy}}{\partial y} = 0$$

which then becomes the Laplacian for displacements, with an equivalent body force on the RHS.

For plane strain, I imagine it should reduce to a bi-harmonic on LHS and equivalent body forces on the RHS. I haven't done this on paper since I was working in anti-plane this entire time, but that's on the top of my TODO list.

[brendanjmeade]

Ahh, this explains it! Happy to work through the plane strain and 3D cases when you are. I'm betting all of this is either in Malvern or standard seismology textbook too.

[mallickrishg]

Below is a derivation of the biharmonic for plane strain, and in that case the equivalent body sources appear on the RHS and as long as we have solutions for point forces I.e., the Kelvin kernels, we can construct any arbitrary shaped force functions.

[brendanjmeade]

Nice! I'll have to think a bit about what bit body force terms look like.