A Stress-strain; strain-displacement in 2-D; uniqueness

[pyroknife]

I am working on a 2-D planar problem in the x-y direction, dealing with stresses, strain, displacements. Under the linear elastic relation and after substitution I can write the following:

$\begin{bmatrix} \sigma_{xx} & \sigma_{xy} \\ \sigma_{xy} & \sigma_{yy} \end{bmatrix} = \mu \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \end{bmatrix} + \mu \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial v}{\partial x}\\ \frac{\partial u}{\partial y} & \frac{\partial v}{\partial y} \end{bmatrix} + \lambda\begin{bmatrix} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} & 0\\ 0 & \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \end{bmatrix}$

I need to obtain the unknowns, which are the 4 displacement gradients
$\frac{du}{dx}; \frac{du}{dy}; \frac{dv}{dx}; \frac{dv}{dy}$

However, since the stress matrix is symmetrical, I only have 3 knowns. In structural problems, where does this 4th condition come from?

 

[Chestermiller]

In structural problems, you don't solve for the displacement gradients. You solve for the displacements.

 

[pyroknife]

In structural problems, you don't solve for the displacement gradients. You solve for the displacements.

This is for a finite difference numerical solution. I need the displacement gradient in order to compute the displacement. The stress is being specified at a boundary face (and each side of the face is a cell). SO I need to compute the gradient at that face and then knowing the gradient and one of the cell values, I can compute the other cell's value.

But either way, doesn't the solving require one more condition to make the system unique?

 

[Chestermiller]

This is for a finite difference numerical solution. I need the displacement gradient in order to compute the displacement. The stress is being specified at a boundary face (and each side of the face is a cell). SO I need to compute the gradient at that face and then knowing the gradient and one of the cell values, I can compute the other cell's value.

But either way, doesn't the solving require one more condition to make the system unique?

Maybe you could tell us the specific problem you're solving?

 

[pyroknife]

Maybe you could tell us the specific problem you're solving?

I'm trying to simulate linear elastic deformation on a square flat plate using the finite difference method. The plate is of size 1x1x0.1. So it has 2 square faces and 4 long rectangular faces. The 2 square sides are oriented in the x-y direction. Assume the 0.1 dimension is in the z direction. 2 of the square and 2 of the rectangular faces are fixed with 0 displacement boundary condition. Then let's say there's a compressive load applied to the plate on 2 of the rectangular faces; 1 in the -x direction and one in the -y direction.

This compressive load acts as a boundary condition on those 2 rectangular sides. The square plate is discretized into 10 cells of size 0.1x0.1x0.1. I am solving for displacements.

Since this compressive load is applied to the face as a stress BC; I need solve the system in the OP to obtain the displacement gradients at the face, which I will use to compute the displacements in the ghost cells, and the ghost cell displacements are in turn used to compute the cell centroid gradient of each cell of the real domain.

 

[pyroknife]

So for example, let's consider the load applied in the -x direction's face.
$\begin{bmatrix} -1 & 0 \\ 0 & \ 0 \end{bmatrix} = \mu \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \end{bmatrix} + \mu \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial v}{\partial x}\\ \frac{\partial u}{\partial y} & \frac{\partial v}{\partial y} \end{bmatrix} + \lambda\begin{bmatrix} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} & 0\\ 0 & \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \end{bmatrix}$

Due to the symmetry of the stress tensor, this is NOT a linear independent system of equations. We have 3 knowns (3 known stresses) and 4 unknowns (4 unknown displacements), which would yield a non-unique solution for the displacement gradients. So I was thinking that I need to impose an addition condition such that a unique solution can be obtained.

 

[Chestermiller]

I'm trying to simulate linear elastic deformation on a square flat plate using the finite difference method. The plate is of size 1x1x0.1. So it has 2 square faces and 4 long rectangular faces. The 2 square sides are oriented in the x-y direction. Assume the 0.1 dimension is in the z direction. 2 of the square and 2 of the rectangular faces are fixed with 0 displacement boundary condition. Then let's say there's a compressive load applied to the plate on 2 of the rectangular faces; 1 in the -x direction and one in the -y direction.

This compressive load acts as a boundary condition on those 2 rectangular sides.

So the problem is this:
:

with the arrows distributed uniformly over each of their rectangular faces, correct? What are the boundary conditions on the other two rectangular faces?

 

[pyroknife]

So the problem is this:
:View attachment 115027
with the arrows distributed uniformly over each of their rectangular faces, correct? What are the boundary conditions on the other two rectangular faces?

Yes! I actually just drew one out on paper. The boundary conditions on the other 2 rectangular faces AND the square faces are fixed with 0 displacement in all directions.

It's not a very realistic problem.

 

[pyroknife]

hmm so I just stumbled across this:
http://bluebox.ippt.pan.pl/~tzielins/doc/ICMM_TGZielinski_Elasticity.paper.pdf
Page 5, section 3.2. Not positive, but maybe this compatibility equation is the additional condition?

 

[Chestermiller]

Yes! I actually just drew one out on paper. The boundary conditions on the other 2 rectangular faces AND the square faces are fixed with 0 displacement in all directions.

It's not a very realistic problem.

OK. Let's assign coordinates. The square faces are at z = +h/2 and z = - h/2, where h = 0.1 m. The rectangular faces are at x = 0, x = L, y = 0, and y = L where L = 1 meter. So your are saying that,

at x = 0, u = v = 0,
at x = L, v = 0,
at y = 0, u = v = 0,
at y = L, u = 0,
at z = + h/2, u = v = 0
at z = - h/2, u = v = 0

Do you realize that the last two conditions require that u and v must be functions of z?

 

[Chestermiller]

hmm so I just stumbled across this:
http://bluebox.ippt.pan.pl/~tzielins/doc/ICMM_TGZielinski_Elasticity.paper.pdf
Page 5, section 3.2. Not positive, but maybe this compatibility equation is the additional condition?

Whoa. Hold your horses. Let's first make sure you are trying to solve the problem that you wish to solve.

 

[pyroknife]

OK. Let's assign coordinates. The square faces are at z = +h/2 and z = - h/2, where h = 0.1 m. The rectangular faces are at x = 0, x = L, y = 0, and y = L where L = 1 meter. So your are saying that,

at x = 0, u = v = 0,
at x = L, v = 0,
at y = 0, u = v = 0,
at y = L, u = 0,
at z = + h/2, u = v = 0
at z = - h/2, u = v = 0

Do you realize that the last two conditions require that u and v must be functions of z?

Regarding the last question: hmmm, so I'm really only trying to simulate a 2-D problem for this case. If u and v must be functions of z, then everything in the OP turns into a 3-D and all the matrices become 3x3. Maybe I can make the assumption that the thickness in the z-direction is infinite, and not 0.1m, so that I can constrain this problem to 2-D physics only.As for the BCs, 4/6 of those are what I'm imposing. Which are
at x = 0, u = v = 0,
at y = 0, u = v = 0,
at z = + h/2, u = v = 0
at z = - h/2, u = v = 0

However, at x=L and y = L, I plan to impose only the stress boundary conditions, and solve for all the displacements. So u&v are unknowns for both of these faces.

 

[Chestermiller]

Regarding the last question: hmmm, so I'm really only trying to simulate a 2-D problem for this case. If u and v must be functions of z, then everything in the OP turns into a 3-D and all the matrices become 3x3. Maybe I can make the assumption that the thickness in the z-direction is infinite, and not 0.1m, so that I can constrain this problem to 2-D physics only.As for the BCs, 4/6 of those are what I'm imposing. Which are
at x = 0, u = v = 0,
at y = 0, u = v = 0,
at z = + h/2, u = v = 0
at z = - h/2, u = v = 0

However, at x=L and y = L, I plan to impose only the stress boundary conditions, and solve for all the displacements. So u&v are unknowns for both of these faces.

OK. I think what you want is the following:

at z = +h/2, $\sigma_{xz} = \sigma_{yz} = 0$, w = 0
at z = - h/2, $\sigma_{xz} = \sigma_{yz} = 0$, w = 0

These will be equivalent to the thickness in the z direction being infinite

Also,

at x = 0, u = v = 0,
at y = 0, u = v = 0,
at x = L, $\sigma_{xx}=-\sigma_1,\ \ \sigma_{xy}=0$
at y = L, $\sigma_{yy}=-\sigma_2, \ \ \sigma_{xy}=0$

Is this what you had in mind?

 

[pyroknife]

OK. I think what you want is the following:

at z = +h/2, $\sigma_{xz} = \sigma_{yz} = 0$, w = 0
at z = - h/2, $\sigma_{xz} = \sigma_{yz} = 0$, w = 0

These will be equivalent to the thickness in the z direction being infinite

Also,

at x = 0, u = v = 0,
at y = 0, u = v = 0,
at x = L, $\sigma_{xx}=-\sigma_1,\ \ \sigma_{xy}=0$
at y = L, $\sigma_{yy}=-\sigma_2, \ \ \sigma_{xy}=0$

Is this what you had in mind?

Yes! Exactly.

 

[Chestermiller]

Yes! Exactly.

Gotta go now. Be back later.

By the way, at x = 0 and y = 0, are you sure you wouldn't rather have

at x = 0, u = 0, $\sigma_{xy}=0$
at y = 0, v = 0, $\sigma_{xy}=0$

 

[pyroknife]

Gotta go now. Be back later.

By the way, at x = 0 and y = 0, are you sure you wouldn't rather have

at x = 0, u = 0, $\sigma_{xy}=0$
at y = 0, v = 0, $\sigma_{xy}=0$

hmmm. Trying to figure out what this means. So instead of "fixing" both displacements to be 0, we fix one displacement and have the shear stress be 0.
Is this equivalent to constraining both u&v AND more? Because, e.g., v displacement on the x=0 face can only be caused by shear, and without shear, I expect v=0.

I'm going to be back in a bit too. Thanks for all the help.

 

[Chestermiller]

hmmm. Trying to figure out what this means. So instead of "fixing" both displacements to be 0, we fix one displacement and have the shear stress be 0.
Is this equivalent to constraining both u&v AND more? Because, e.g., v displacement on the x=0 face can only be caused by shear, and without shear, I expect v=0.

No. Zero shear stress is less constraining than zero tangential displacement. The points on the boundary are free to slide along the boundary in order to relieve any shear stress. It's like the boundary being a shower curtain rod, and the rings of the shower curtain being allowed to slide freely along the rod. At DuPont, we used to call this a "shower curtain boundary condition." So, if the shear stress is zero along the boundary, v is not equal to zero.

 

[pyroknife]

No. Zero shear stress is less constraining than zero tangential displacement. The points on the boundary are free to slide along the boundary in order to relieve any shear stress. It's like the boundary being a shower curtain rod, and the rings of the shower curtain being allowed to slide freely along the rod. At DuPont, we used to call this a "shower curtain boundary condition." So, if the shear stress is zero along the boundary, v is not equal to zero.

Ah I see that is interesting. So how would you obtain v along the boundary if you specify the shear stress? I guess this is the same issue as the other 2 compression boundary conditions.

 

[Chestermiller]

Ah I see that is interesting. So how would you obtain v along the boundary if you specify the shear stress? I guess this is the same issue as the other 2 compression boundary conditions.

We'll get to all that. But first, which of the two boundary conditions do you prefer using.

 

[pyroknife]

We'll get to all that. But first, which of the two boundary conditions do you prefer using.

Hmm Let me stick with my original one with u=v=0 for now. Setting shear stress=0 as opposed to both displacements=0 makes the computation harder. After figuring out the compressive BC computations, I would think that I would be able to solve the shear=0 BC.

 

[Chestermiller]

Hmm Let me stick with my original one with u=v=0 for now. Setting shear stress=0 as opposed to both displacements=0 makes the computation harder. After figuring out the compressive BC computations, I would think that I would be able to solve the shear=0 BC.

OK. So the stresses are:
$$\sigma_{xx}=(2\mu +\lambda)\frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial y}$$
$$\sigma_{yy}=(2\mu +\lambda)\frac{\partial v}{\partial y}+\lambda \frac{\partial u}{\partial x}$$
$$\sigma_{xy}=\mu\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)$$

The stress equilibrium equations are given by:
$$\frac{\partial \sigma_{xx}}{\partial x}+\frac{\partial \sigma_{xy}}{\partial y}=0$$
$$\frac{\partial \sigma_{xy}}{\partial x}+\frac{\partial \sigma_{yy}}{\partial y}=0$$

Now, what do you get when you substitute the equations for the stresses into the stress equilibrium equations?

 

[pyroknife]

OK. So the stresses are:
$$\sigma_{xx}=(2\mu +\lambda)\frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial y}$$
$$\sigma_{yy}=(2\mu +\lambda)\frac{\partial v}{\partial y}+\lambda \frac{\partial u}{\partial x}$$
$$\sigma_{xy}=\mu\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)$$

The stress equilibrium equations are given by:
$$\frac{\partial \sigma_{xx}}{\partial x}+\frac{\partial \sigma_{xy}}{\partial y}=0$$
$$\frac{\partial \sigma_{xy}}{\partial x}+\frac{\partial \sigma_{yy}}{\partial y}=0$$

Now, what do you get when you substitute the equations for the stresses into the stress equilibrium equations?

After substitution, I get the following:
$$ (2\mu+\lambda)\frac{\partial^2 u}{\partial x^2} +(\mu+\lambda) \frac{\partial^2 v}{\partial x\partial y} + \mu \frac{\partial^2 u}{\partial y^2} = 0$$
$$ (2\mu+\lambda)\frac{\partial^2 v}{\partial y^2} +(\mu+\lambda) \frac{\partial^2 u}{\partial x\partial y} + \mu \frac{\partial^2 v}{\partial x^2} = 0$$

 

[Chestermiller]

After substitution, I get the following:
$$ (2\mu+\lambda)\frac{\partial^2 u}{\partial x^2} +(\mu+\lambda) \frac{\partial^2 v}{\partial x\partial y} + \mu \frac{\partial^2 u}{\partial y^2} = 0$$
$$ (2\mu+\lambda)\frac{\partial^2 v}{\partial y^2} +(\mu+\lambda) \frac{\partial^2 u}{\partial x\partial y} + \mu \frac{\partial^2 v}{\partial x^2} = 0$$

Excellent. So now you have two partial differential equations in the 2 unknown displacements, u and v. These are the equations that you can solve numerically. Next, work out the boundary conditions in terms of u, v, and their derivatives. Let's see what you come up with.

Incidentally, for the case in which the shear stress is zero on all the boundaries, you get a very simple analytic solution. See if you can derive it. It might be a good test case for your numerical method.

 

[pyroknife]

Excellent. So now you have two partial differential equations in the 2 unknown displacements, u and v. These are the equations that you can solve numerically. Next, work out the boundary conditions in terms of u, v, and their derivatives. Let's see what you come up with.

Incidentally, for the case in which the shear stress is zero on all the boundaries, you get a very simple analytic solution. See if you can derive it. It might be a good test case for your numerical method.

So for now, let me just consider the x=L boundary condition. Aren't the equations that I need on post #6 on the previous page? However that post is written in terms of displacement gradients.

 

[pyroknife]

The x=L BC gives the following:
$$
\begin{bmatrix}
-1 & 0 \\
0 & \ 0
\end{bmatrix} = \mu
\begin{bmatrix}
\frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\
\frac{\partial v}{\partial x} & \frac{\partial v}{\partial y}
\end{bmatrix}
+ \mu
\begin{bmatrix}
\frac{\partial u}{\partial x} & \frac{\partial v}{\partial x}\\
\frac{\partial u}{\partial y} & \frac{\partial v}{\partial y}
\end{bmatrix}
+ \lambda\begin{bmatrix}
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} & 0\\
0 & \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}
\end{bmatrix}
$$

Or
$$
(2\mu+\lambda)\frac{\partial u}{\partial x} + \lambda \frac{\partial u}{\partial y} = -1 \\
\mu (\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} = 0 \\
(2\mu+\lambda) \frac{\partial v}{\partial y} + \lambda \frac{\partial u}{\partial x} = 0 \\
$$

But these are only 3 equations with 4 unknowns.

Using a central difference, I get the following:
$$
\frac{\partial u}{\partial x} = \frac{u_{i+1,j} - u_{i-1,j}}{2\Delta x}
$$

where i&j are the indices for the x&y directions, respectively.

 

[Chestermiller]

The x=L BC gives the following:
$$
\begin{bmatrix}
-1 & 0 \\
0 & \ 0
\end{bmatrix} = \mu
\begin{bmatrix}
\frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\
\frac{\partial v}{\partial x} & \frac{\partial v}{\partial y}
\end{bmatrix}
+ \mu
\begin{bmatrix}
\frac{\partial u}{\partial x} & \frac{\partial v}{\partial x}\\
\frac{\partial u}{\partial y} & \frac{\partial v}{\partial y}
\end{bmatrix}
+ \lambda\begin{bmatrix}
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} & 0\\
0 & \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}
\end{bmatrix}
$$

Or
$$
(2\mu+\lambda)\frac{\partial u}{\partial x} + \lambda \frac{\partial u}{\partial y} = -1 \\
\mu (\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} = 0 \\
(2\mu+\lambda) \frac{\partial v}{\partial y} + \lambda \frac{\partial u}{\partial x} = 0 \\
$$

But these are only 3 equations with 4 unknowns.

No. At x = L, the BCs are:

$v=0$

and

$\sigma_{xx}=(2\mu + \lambda)\frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial y}=-\sigma_1$
But, the 2nd term is equal to zero, so:
$$\frac{\partial u}{\partial x}=-\frac{\sigma_1}{(2\mu+\lambda)}$$

So you know the value of v and the normal derivative of u.

 

[pyroknife]

No. At x = L, the BCs are:

$v=0$

and

$\sigma_{xx}=(2\mu + \lambda)\frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial y}=-\sigma_1$
But, the 2nd term is equal to zero, so:
$$\frac{\partial u}{\partial x}=-\frac{\sigma_1}{(2\mu+\lambda)}$$

So you know the value of v and the normal derivative of u.

Wait, hold on. I was only going to define
$$
\sigma_{xx} = -\sigma_1
$$
at the x=L BC, and not constrain v = 0. Can I do this or is this insufficiently defined?
Or in other words, I'd be solving for v and u, given only the stress BC

 

[pyroknife]

Wait, hold on. I was only going to define
$$
\sigma_{xx} = -\sigma_1
$$
at the x=L BC, and not constrain v = 0. Can I do this or is this insufficiently defined?
Or in other words, I'd be solving for v and u, given only the stress BC

Sorry, let me correct myself. I mean that I would define
$$
\sigma_{xx} = -\sigma_1 \\
\sigma_{xy} = 0 \\
\sigma_{yy} = 0
$$

 

[Chestermiller]

Wait, hold on. I was only going to define
$$
\sigma_{xx} = -\sigma_1
$$
at the x=L BC, and not constrain v = 0. Can I do this or is this insufficiently defined?
Or in other words, I'd be solving for v and u, given only the stress BC

Oops. I meant for the shear stress to be zero. Sorry. So,
$$\frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}=0$$and$$\sigma_{xx}=(2\mu+\lambda)\frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial y}=-\sigma_1$$So, the normal derivatives of the displacements are:
$$\frac{\partial v}{\partial x}=-\frac{\partial u}{\partial y}$$ and
$$\frac{\partial u}{\partial x}=-\frac{\sigma_1}{(2\mu+\lambda)}-\frac{\lambda}{(2\mu+\lambda)}\frac{\partial v}{\partial y}$$

 

[Chestermiller]

Sorry, let me correct myself. I mean that I would define
$$
\sigma_{xx} = -\sigma_1 \\
\sigma_{xy} = 0 \\
\sigma_{yy} = 0
$$

$\sigma_{yy}$ can't be specified at x = L. The most you can get is the normal stress and the shear stress on a surface. These are the components of the stress vector.

 

[Chestermiller]

Yes, don't we have 4 unknowns in $$\frac{du}{dx} \\ \frac{du}{dy} \\ \frac{dv}{dx} \\ \frac{dv}{dy} $$? And not enough knowns?

There are only two unknowns: the two displacements u and v.

 

[pyroknife]

There are only two unknowns: the two displacements u and v.

Oops I deleted that post after reading your post prior to the quoted one.
I didn't realize that you can't specify $\sigma_{yy}$ on that face. So yes, I agree that u&v are the unknowns, but to obtain u&v, wouldn't I need to obtain the displacement gradients, and then apply the finite difference formula to compute u&v?

 

[Chestermiller]

Oops I deleted that post after reading your post prior to the quoted one.
I didn't realize that you can't specify $\sigma_{yy}$ on that face. So yes, I agree that u&v are the unknowns, but to obtain u&v, wouldn't I need to obtain the displacement gradients, and then apply the finite difference formula to compute u&v?

No. You are going to be solving the two differential equations for u and v in post #22, subject to the boundary conditions on u and v.

 

[pyroknife]

No. You are going to be solving the two differential equations for u and v in post #22, subject to the boundary conditions on u and v.

Yes, but the issue is when I discretize the domain into cells, there are issues when I get to cells that lie on the boundary. I plan to use a central difference for the second derivatives, so I would need to know the u&v values in ghost cells that lie outside the real domain. So before I can compute anything, I would need to solve for the u&v values in all the ghost cells given the BCs. I am applying the stress BCs on the face between the boundary cell and the ghost cell, and the displacement BCs are applied directly in the ghost cells. So, for the stress BCs, I do not know the u&v values in the ghost cell. Wouldn't I need to compute the u&v values in the ghost cells first using the stress BC?

 

[Chestermiller]

Yes, but the issue is when I discretize the domain into cells, there are issues when I get to cells that lie on the boundary. I plan to use a central difference for the second derivatives, so I would need to know the u&v values in ghost cells that lie outside the real domain. So before I can compute anything, I would need to solve for the u&v values in all the ghost cells given the BCs. I am applying the stress BCs on the face between the boundary cell and the ghost cell, and the displacement BCs are applied directly in the ghost cells. So, for the stress BCs, I do not know the u&v values in the ghost cell. Wouldn't I need to compute the u&v values in the ghost cells first using the stress BC?

No. Now that we have everything on the boundaries expressed in terms of the displacements and their derivatives, we no longer have to be working in terms of the stresses. This is where having some experience with solving PDEs numerically comes in. I can help with this. The only problem is properly handling the BCs in conjunction with the finite difference approximations to the differential equations at x = L and y = L.

It's late now, so I'll be back tomorrow. Meanwhile, please write out finite difference approximations to the differential equations at the interior grid points.

 

[pyroknife]

No. Now that we have everything on the boundaries expressed in terms of the displacements and their derivatives, we no longer have to be working in terms of the stresses. This is where having some experience with solving PDEs numerically comes in. I can help with this. The only problem is properly handling the BCs in conjunction with the finite difference approximations to the differential equations at x = L and y = L.

It's late now, so I'll be back tomorrow. Meanwhile, please write out finite difference approximations to the differential equations at the interior grid points.

Thanks for all the help. I will type it out in the next post here in a little bit.

 

[pyroknife]

The central difference formulas are:
\begin{equation}
\frac{\partial ^2 u}{\partial x^2} \approx \frac{u_{i+1,j}-2u_{i,j} + u_{i-1,j}}{\Delta x^2}
\end{equation}
\begin{equation}
\frac{\partial ^2 u}{\partial x \partial y} \approx \frac{u_{i+1,j+1}-u_{i+1,j-1} - u_{i-1,j+1} + u_{i-1,j-1}}{4\Delta x \Delta y}
\end{equation}

Plugging these into the PDEs yield:

\begin{equation}
(2\mu + \lambda) \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{\Delta x^2} + (\mu+\lambda)\frac{v_{i+1,j+1} - v_{i+1,j-1} - v_{i-1,j+1} + v_{i-1,j-1}}{4\Delta x\Delta y} + \mu \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{\Delta y^2} = 0
\end{equation}
\begin{equation}
(2\mu + \lambda) \frac{v_{i,j+1} - 2v_{i,j} + u_{i,j-1}}{\Delta y^2} + (\mu+\lambda)\frac{u_{i+1,j+1} - u_{i+1,j-1} - u_{i-1,j+1} + u_{i-1,j-1}}{4\Delta x\Delta y} + \mu \frac{v_{i+1,j} - 2v_{i,j} + v_{i-1,j}}{\Delta x^2} = 0
\end{equation}

these equations are valid from i,j=1,N. $\Delta x$ and $\Delta y$ are the spatial step size between grid points and are constant.
Assuming that if i=1 or j=1 or i=N, j=N is the boundary node. e.g., i = 0=>x=0, j=0=>y=0,i=N=>x=L,j=N=>y=L. Assuming that the displacement boundary conditions are applied at i = 0 and j = 0, which lie off the real domain, i.e., $x=-\Delta x$, $y=-\Delta y$. The stress BCs are applied half way in between i,j = N+1 and i,j=N. i.e., $x = L + \Delta x/2$, $y=L + \Delta y/2$

In the above finite difference implementation issues would arise at i or j = N, since we do not know the u or v values at i or j = N+1.

 

[Chestermiller]

The central difference formulas are:
\begin{equation}
\frac{\partial ^2 u}{\partial x^2} \approx \frac{u_{i+1,j}-2u_{i,j} + u_{i-1,j}}{\Delta x^2}
\end{equation}
\begin{equation}
\frac{\partial ^2 u}{\partial x \partial y} \approx \frac{u_{i+1,j+1}-u_{i+1,j-1} - u_{i-1,j+1} + u_{i-1,j-1}}{4\Delta x \Delta y}
\end{equation}

Plugging these into the PDEs yield:

\begin{equation}
(2\mu + \lambda) \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{\Delta x^2} + (\mu+\lambda)\frac{v_{i+1,j+1} - v_{i+1,j-1} - v_{i-1,j+1} + v_{i-1,j-1}}{4\Delta x\Delta y} + \mu \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{\Delta y^2} = 0
\end{equation}
\begin{equation}
(2\mu + \lambda) \frac{v_{i,j+1} - 2v_{i,j} + u_{i,j-1}}{\Delta y^2} + (\mu+\lambda)\frac{u_{i+1,j+1} - u_{i+1,j-1} - u_{i-1,j+1} + u_{i-1,j-1}}{4\Delta x\Delta y} + \mu \frac{v_{i+1,j} - 2v_{i,j} + v_{i-1,j}}{\Delta x^2} = 0
\end{equation}

these equations are valid from i,j=1,N. $\Delta x$ and $\Delta y$ are the spatial step size between grid points and are constant.
Assuming that if i=1 or j=1 or i=N, j=N is the boundary node. e.g., i = 0=>x=0, j=0=>y=0,i=N=>x=L,j=N=>y=L. Assuming that the displacement boundary conditions are applied at i = 0 and j = 0, which lie off the real domain, i.e., $x=-\Delta x$, $y=-\Delta y$. The stress BCs are applied half way in between i,j = N+1 and i,j=N. i.e., $x = L + \Delta x/2$, $y=L + \Delta y/2$

In the above finite difference implementation issues would arise at i or j = N, since we do not know the u or v values at i or j = N+1.

Nice job.

HANDLING THE BOUNDARY CONDITION AT x= L
We set the grid up so that there are grid points along the boundaries x = L and y = L. We need to derive appropriate finite difference relationships that apply at these boundaries so that we can determine the displacements on these boundaries. At x = L, we showed that the boundary conditions are as follows:
$$\frac{\partial v}{\partial x}=-\frac{\partial u}{\partial y}\tag{1}$$ and
$$\frac{\partial u}{\partial x}=-\frac{\sigma_1}{(2\mu+\lambda)}-\frac{\lambda}{(2\mu+\lambda)}\frac{\partial v}{\partial y}\tag{2}$$
With no loss of generality, we can take the partial derivatives of these equations in the direction tangent to the boundary (i.e., the y direction) to obtain:
$$\frac{\partial^2 v}{\partial x\partial y}=-\frac{\partial^2 u}{\partial y^2}\tag{3}$$ and
$$\frac{\partial^2 u}{\partial x\partial y}=-\frac{\lambda}{(2\mu+\lambda)}\frac{\partial^2 v}{\partial y^2}\tag{4}$$
If we substitute these relationships into the stress equilibrium equations, we obtain (specifically for the boundary at x = L):
$$ \frac{\partial^2 u}{\partial x^2} + \frac{\lambda}{(2\mu +\lambda)} \frac{\partial^2 u}{\partial y^2} = 0\tag{5}$$
$$ \frac{\partial^2 v}{\partial x^2}+\left(2+\frac{\lambda}{(2\mu +\lambda)}\right)\frac{\partial^2 v}{\partial y^2} = 0\tag{6}$$
If we follow the usual procedure of introducing a fictitious grid point at $x=L+\Delta x$, then we can use the boundary conditions to obtain finite difference approximations to the dependent variable values at the fictitious grid point:
$$u_{N+1,j}\approx u_{N-1,j}-2\Delta x\frac{\sigma_1}{(2\mu +\lambda)}-\frac{\lambda}{(2\mu +\lambda)} \frac{\Delta x}{\Delta y}(v_{N,j+1}-v_{N,j-1}) \tag{7}$$
$$v_{N+1,j}\approx v_{N-1,j}-\frac{\lambda}{(2\mu +\lambda)} \frac{\Delta x}{\Delta y}(u_{N,j+1}-u_{N,j-1})\tag{8} $$
If we substitute these into our finite difference approximations for $\partial^2 u/\partial x^2$ and $\partial^2 v/\partial x^2$, we obtain:
$$\frac{\partial^2u}{\partial x^2}\approx \frac{2(u_{N-1,j}-u_{N,j})}{(\Delta x)^2}-\frac{2}{\Delta x}\frac{\sigma_1}{(2\mu +\lambda)}-\frac{\lambda}{(2\mu +\lambda)}\frac{(v_{N,j+1}-v_{N,j-1})}{\Delta y} \tag{9}$$
$$\frac{\partial^2v}{\partial x^2}\approx \frac{2(v_{N-1,j}-v_{N,j})}{(\Delta x)^2}-\frac{\lambda}{(2\mu +\lambda)}\frac{(u_{N,j+1}-u_{N,j-1})}{\Delta y} \tag{10}$$

Comments?

 

[pyroknife]

Nice job.

HANDLING THE BOUNDARY CONDITION AT x= L
We set the grid up so that there are grid points along the boundaries x = L and y = L. We need to derive appropriate finite difference relationships that apply at these boundaries so that we can determine the displacements on these boundaries. At x = L, we showed that the boundary conditions are as follows:
$$\frac{\partial v}{\partial x}=-\frac{\partial u}{\partial y}\tag{1}$$ and
$$\frac{\partial u}{\partial x}=-\frac{\sigma_1}{(2\mu+\lambda)}-\frac{\lambda}{(2\mu+\lambda)}\frac{\partial v}{\partial y}\tag{2}$$
With no loss of generality, we can take the partial derivatives of these equations in the direction tangent to the boundary (i.e., the y direction) to obtain:
$$\frac{\partial^2 v}{\partial x\partial y}=-\frac{\partial^2 u}{\partial y^2}\tag{3}$$ and
$$\frac{\partial^2 u}{\partial x\partial y}=-\frac{\lambda}{(2\mu+\lambda)}\frac{\partial^2 v}{\partial y^2}\tag{4}$$
If we substitute these relationships into the stress equilibrium equations, we obtain (specifically for the boundary at x = L):
$$ \frac{\partial^2 u}{\partial x^2} + \frac{\lambda}{(2\mu +\lambda)} \frac{\partial^2 u}{\partial y^2} = 0\tag{5}$$
$$ \frac{\partial^2 v}{\partial x^2}+\left(2+\frac{\lambda}{(2\mu +\lambda)}\right)\frac{\partial^2 v}{\partial y^2} = 0\tag{6}$$
If we follow the usual procedure of introducing a fictitious grid point at $x=L+\Delta x$, then we can use the boundary conditions to obtain finite difference approximations to the dependent variable values at the fictitious grid point:
$$u_{N+1,j}\approx u_{N-1,j}-2\Delta x\frac{\sigma_1}{(2\mu +\lambda)}-\frac{\lambda}{(2\mu +\lambda)} \frac{\Delta x}{\Delta y}(v_{N,j+1}-v_{N,j-1}) \tag{7}$$
$$v_{N+1,j}\approx v_{N-1,j}-\frac{\lambda}{(2\mu +\lambda)} \frac{\Delta x}{\Delta y}(u_{N,j+1}-u_{N,j-1})\tag{8} $$
If we substitute these into our finite difference approximations for $\partial^2 u/\partial x^2$ and $\partial^2 v/\partial x^2$, we obtain:
$$\frac{\partial^2u}{\partial x^2}\approx \frac{2(u_{N-1,j}-u_{N,j})}{(\Delta x)^2}-\frac{2}{\Delta x}\frac{\sigma_1}{(2\mu +\lambda)}-\frac{\lambda}{(2\mu +\lambda)}\frac{(v_{N,j+1}-v_{N,j-1})}{\Delta y} \tag{9}$$
$$\frac{\partial^2v}{\partial x^2}\approx \frac{2(v_{N-1,j}-v_{N,j})}{(\Delta x)^2}-\frac{\lambda}{(2\mu +\lambda)}\frac{(u_{N,j+1}-u_{N,j-1})}{\Delta y} \tag{10}$$

Comments?

I agree! However, n this scenario, aren't the the stress BCs being applied at x,y=N? If so, then in Eqns (9) and (10), aren't the u&v values unknown at N?

 

[Chestermiller]

I agree! However, n this scenario, aren't the the stress BCs being applied at x,y=N? If so, then in Eqns (9) and (10), aren't the u&v values unknown at N?

Sure. We are solving for them along with the rest of ythe u's and v's at all the other grid points.

 

[pyroknife]

Sure. We are solving for them along with the rest of ythe u's and v's at all the other grid points.

I should have mentioned this early, but I am solving the transient equations, which includes the time derivative term. So instead of the 2 PDEs in post #22 being set equal to 0 on the RHS, they are instead set equal to $\rho \frac{\partial^2 u}{\partial t^2}$. In this case, I think I would need to know the u&v values at x&y=N?

 

[Chestermiller]

I should have mentioned this early, but I am solving the transient equations, which includes the time derivative term. So instead of the 2 PDEs in post #22 being set equal to 0 on the RHS, they are instead set equal to $\rho \frac{\partial^2 u}{\partial t^2}$. In this case, I think I would need to know the u&v values at x&y=N?

I don't think so. What is the nature of the transient loading?

 

[pyroknife]

I don't think so. What is the nature of the transient loading?

I am working on a research project for modeling thermal protection systems that involves 2 major components: (1) Thermochemical physics modeling and to a lesser extent (2) Structural mechanics modeling (linear elastic should be sufficient). The thermochemical modeling involves solving transient solid decomposition (arrhenius rate equation), continuity, and mixture-energy equations. I am working with a class of materials where the thermo properties such as thermal conductivity, porosity, permeability are expected to be a function of stress/strain.

I'm not too familiar with structural mechanics, but when would I be solving the steady state as opposed to the transient equation?

 

[Chestermiller]

I am working on a research project for modeling thermal protection systems that involves 2 major components: (1) Thermochemical physics modeling and to a lesser extent (2) Structural mechanics modeling (linear elastic should be sufficient). The thermochemical modeling involves solving transient solid decomposition (arrhenius rate equation), continuity, and mixture-energy equations. I am working with a class of materials where the thermo properties such as thermal conductivity, porosity, permeability are expected to be a function of stress/strain.

I'm not too familiar with structural mechanics, but when would I be solving the steady state as opposed to the transient equation?

You would be solving the steady state equation if things aren't changing with time. Otherwise, you need to include the time-dependence.

 

[pyroknife]

You would be solving the steady state equation if things aren't changing with time. Otherwise, you need to include the time-dependence.

Aren't the displacement fields always changing with time?

 

[Chestermiller]

Aren't the displacement fields always changing with time?

No way. If you exert a tensile force on a rod, is the displacement field always changing? If you hang a mass from a spring (at a point that it does not oscillate), is the displacement field of the points along the spring always changing?

 

[pyroknife]

No way. If you exert a tensile force on a rod, is the displacement field always changing? If you hang a mass from a spring (at a point that it does not oscillate), is the displacement field of the points along the spring always changing?

No. Hmmm. So if the boundary conditions are not a function of time, does that essentially mean the displacement field will not change?

Also, technically the density of the material I am modeling would change with time because as I mentioned the thermophysics part of the code models solid decomposition. The level of decomposition (level of density change) would depend on the material being modeled. So this could be a mechanism that could make the structural mechanics transient. Furthermore, I am also thinking about incorporating thermal strain into the structural computations, and thermal strain is a function of temperature, which is time-independent in the thermophysics part of the code. In addition, the internal pressure of the material may also be changing due to internal gas flow (the materials I am working with are usually highly porous), and this would change the stress I presume.

 

[Chestermiller]

OMG

No. Hmmm. So if the boundary conditions are not a function of time, does that essentially mean the displacement field will not change?

You are familiar with the term "static equilibrium," correct? We covered this in freshman physics. You also probably had a course in Statics and Dynamics, which again covered static equilibrium. Do you see how this relates to what you are asking?

Also, technically the density of the material I am modeling would change with time because as I mentioned the thermophysics part of the code models solid decomposition. The level of decomposition (level of density change) would depend on the material being modeled. So this could be a mechanism that could make the structural mechanics transient.

Where is the mass going?

Furthermore, I am also thinking about incorporating thermal strain into the structural computations, and thermal strain is a function of temperature, which is time-independent in the thermophysics part of the code.

So you are including transient heat conduction in your code? If so, then, yes, thermal strain makes the stresses and elastic strains transient. If the temperature is at steady state, however, they you are back to static equilibrium.

In addition, the internal pressure of the material may also be changing due to internal gas flow (the materials I am working with are usually highly porous), and this would change the stress I presume.

Yes. If you are going to be including this, then you need to familiarize yourself with the subject of poroelasticity. You need to learn about the Terzaghi effective stress, and how to describe the coupling between the fluid flow and the solid mechanics.

From what you describe above, I am beginning to sense a geophysical application, involving subsurface flow of liquids, and the effect of increased pore pressure on the solid mechanics (and the feedback of the solid mechanics on the storativity of the rock).

 

[pyroknife]

OMGYou are familiar with the term "static equilibrium," correct? We covered this in freshman physics. You also probably had a course in Statics and Dynamics, which again covered static equilibrium. Do you see how this relates to what you are asking?

Where is the mass going?

So you are including transient heat conduction in your code? If so, then, yes, thermal strain makes the stresses and elastic strains transient. If the temperature is at steady state, however, they you are back to static equilibrium.

Yes. If you are going to be including this, then you need to familiarize yourself with the subject of poroelasticity. You need to learn about the Terzaghi effective stress, and how to describe the coupling between the fluid flow and the solid mechanics.

From what you describe above, I am beginning to sense a geophysical application, involving subsurface flow of liquids, and the effect of increased pore pressure on the solid mechanics (and the feedback of the solid mechanics on the storativity of the rock).

Yes I know what static equilibrium means. Essentially some body at rest in which net forces sum to 0. I took a statics and dynamics class long ago. I'm still kind of confused on applied loading. With your tensile force applied to a rod example, if you apply this constant tensile force for 10minutes, wouldn't the displacement field within the rod change for the duration the force is applied?

The material is made up of a solid matrix where its pores are resin infused with some resin. At high temperatures, the resin decomposes and diffuses through the pores into the boundary layer. In addition, the solid matrix also ablates at the surface, resulting in mass law at the surface and a change in the computational domain.

Transient heat conduction is included. At some point the temperature does become steady state obviously.

I don't know if I want to go too in-depth with the structural modeling as my focus is in thermophysics&fluids and I don't know very much about structural mechanics at all (and frankly don't have nearly as much interest in it compared to thermophysics and fluids). Poroelasticity sounds like it would primarily involve microscale modeling, which could be interesting and microscale modeling is definitely something I think I will need to do in the future. I just briefly google'd poroelasticity.
"Modeling poroelasticity requires the coupling of two laws. The first of these is Darcy's law...The second law is the structural displacement of the porous matrix. Biot poroelasticity describes this coupled physics." I'm currently using Darcy's law to compute the velocity of the porous gas flow in the continuity equation. And the linear elastic equation solves for displacement.


This is actually an atmospheric entry application for space exploration. A lot of the thermal protection systems (e.g., PICA from NASA) found on modern atmospheric entry vehicles/probes are resin-infused.

 

[Chestermiller]

Yes I know what static equilibrium means. Essentially some body at rest in which net forces sum to 0. I took a statics and dynamics class long ago. I'm still kind of confused on applied loading. With your tensile force applied to a rod example, if you apply this constant tensile force for 10minutes, wouldn't the displacement field within the rod change for the duration the force is applied?

No. The solid would respond virtually instantaneously.

The material is made up of a solid matrix where its pores are resin infused with some resin. At high temperatures, the resin decomposes and diffuses through the pores into the boundary layer. In addition, the solid matrix also ablates at the surface, resulting in mass law at the surface and a change in the computational domain.

Transient heat conduction is included. At some point the temperature does become steady state obviously.

This all sounds very interesting. I'm a ChE, so this is the kind of thing I respond to.

I don't know if I want to go too in-depth with the structural modeling as my focus is in thermophysics&fluids and I don't know very much about structural mechanics at all (and frankly don't have nearly as much interest in it compared to thermophysics and fluids). Poroelasticity sounds like it would primarily involve microscale modeling, which could be interesting and microscale modeling is definitely something I think I will need to do in the future. I just briefly google'd poroelasticity.
"Modeling poroelasticity requires the coupling of two laws. The first of these is Darcy's law...The second law is the structural displacement of the porous matrix. Biot poroelasticity describes this coupled physics." I'm currently using Darcy's law to compute the velocity of the porous gas flow in the continuity equation. And the linear elastic equation solves for displacement.

OK. This helps. I don't think the pore pressures would be high enough to feed back into the solid mechanics.

This is actually an atmospheric entry application for space exploration. A lot of the thermal protection systems (e.g., PICA from NASA) found on modern atmospheric entry vehicles/probes are resin-infused.

Cool!