Incorporating moment-rotation springs to 1D FEA

[LRJ]

Does anyone know how to include moment-rotation springs into a 1D FEA model, or know of a reference which explains this? I've got my own program to model a laterally loaded pile with p-y springs, but now I want to incorporate M-θ springs and was wondering how to implement them. Specifically, the effect of the M-θ springs on the load vectors: can I include the M-θ spring effect as an additional set of loads to the element load vector? Do all terms in the element load vector (F1, M1, F2, M2) get affected, or just the moment terms?

Thanks in advance for any help.

 

[Shu Jiang PEng SE]

I think you can add the spring stiffness of each degree of freedom to the corresponding element of the element stiffness matrix.

 

[LRJ]

Are you saying I should simply add another stiffness term in front of the element mass matrix, [m_M], and multiply against [m_M] accordingly?

So:

(k_p-y + k_M-θ).[m_M]

Rather than the current:

k_p-y.[m_M]

???

Would the units be correct? Moreover, the units of the M-θ gradient are different to the p-y stiffness.

 

[Shu Jiang PEng SE]

[K]{X}={F}
The element stiffness matrix without spring is:
k11 k12 k13 k14 k15 k16
k21 k22 k23 k24 k25 k26
k31 k32 k33 k34 k35 k36
k41 k42 k43 k44 k45 k45
k51 k52 k53 k54 k55 k56
k61 k62 k63 k64 k65 k66
Add spring stiffness the matrix will be like
k11* k12 k13 k14 k15 k16
k21 k22* k23 k24 k25 k26
k31 k32 k33* k34 k35 k36
k41 k42 k43 k44* k45 k45
k51 k52 k53 k54 k55* k56
k61 k62 k63 k64 k65 k66*

k11*=k11+spring k22*=k22+spring

 

[LRJ]

To clarify, I am considering a 1D FEA model of a laterally loaded pile. This has 4 degrees of freedom rather than 6.

Nonetheless, from memory, the governing equation is:

EI.d⁴w/dx⁴ + P.d²w/dx² + k_p-y.dw/dx = q

Where:
EI = Flexural rigidity of pile
k_p-y = Stiffness of p-y curve
P = Axial load

When solving with FEA and using matrix notation:

([k_M] + [m_M] + [g_M]).{y} = {Q}

Where:
[k_M] = Element stiffness matrix (for 'internal' pile forces caused by EI)
[m_M] = Element mass matrix (for 'external' pile forces caused by soil reaction - currently just p-y springs incorporated; p values are distributed loads across the element length)
[g_M] = Element geometric matrix (for inclusion of axial force across elements)
{y} = Displacement vector
{Q} = Load vector

Where do moment-rotation springs fit into the above governing equation? Moreover, what is the governing equation when including moment-rotation springs?

The units of the element stiffness matrix appear incompatible with the stiffness of a moment-rotation curve (units of moment-rotation stiffness would be kNm per radian for a lumped spring or kNm/m per radian for a distributed moment reaction)

 

[gravityandinertia]

I think the easiest way to think about this, is that for a 1D element they will act at the ends. So going to the fundamental equation of [K][D] = [F]

Your moment rotation springs become an external force dependent on displacement, follow me for a single line of the matrix system

k_n*d_n = f_n
k = stiffness
d = displacement
f = force
n is the line in the matrix being represented.

f_n however is a function of d_n
f_n = m_n*d_n
m_n = rotational stiffness

Now substitute
k_n*d_n = m_n*d_n

There is now a problem. You can't solve this system since you need all d_n terms on the left side. So we rearrange
k_n*d_n - m_n*d_n = 0

combine terms to make your stiffness term more apparent
(k_n - m_n)*d_n = 0

So, to include your rotational stiffness you need to subtract off the stiffness of your M-θ springs.

 

[LRJ]

I agree that the effect would be deducted. However, as I mentioned previously, the units of the M-θ stiffness are incompatible with the element stiffness matrix - as I defined it in my previous post - so I do not believe it is simply a case of deducting the gradient of the M-θ spring from the matrix.

Do you have a reference explaining the basis for what you discuss? I think you are both coming at this from a structural engineering perspective where there is only one element matrix (the element stiffness matrix). To model distributed loads along a pile you need other matrices, as I mentioned in my previous post.

 

[gravityandinertia]

I'm confused as to your implementation. Having 4 DOF instead of 6, have you removed the rotational degrees-of-freedom? If so, then you are right, you can't model this with inconsistent units. You need a 6 DOF element to implement this properly.

Also, I'm not coming at this from a single element perspective,the concept is the same for all fea, whether performing a single element or many elements. The global stiffness matrix for a many element system is built entirely from local stiffness matrices of single elements.

 

[gravityandinertia]

I'm re-reading your posts more carefully and I don't really understand your implementation. Where does a mass matrix come into play if you aren't including dynamic forces? If it's to include gravitational forces why would those not just be converted into the external force vector?

 

[LRJ]

Apologies if it was not clear.

I didn't say you were coming at this from a single element perspective. I presumed that you were more familiar with structural analysis where there is typically just the element stiffness matrix rather than the other two matrices I mentioned in my 3rd post in this thread.

The element mass matrix is not just used for dynamic analyses. A beam on an elastic foundation also has an element mass matrix too. For a 4DOF laterally loaded beam-column (pile):

[m_M] =

 156   22L    54   -13L
 22L   4L[sup]2[/sup]   13L   -3L[sup]2[/sup]
 54   13L   156   -22L
-13L  -3L   -22L    4L[sup]2[/sup]

multiplied by k_p-y.L/420

The above is pretty standard for what I'm looking at and is referenced in several places, including 'Programming the Finite Element Method' by Smith, Griffiths and Margetts.

I suspect I could include the M-θ reaction in the element load vector, though I need to know how to resolve a distributed moment. I do wonder whether I need to add another matrix similar to the above, however.

 

[LRJ]

After doing a bit more digging I have found a reference which explains how to incorporate the moment-rotation springs. Essentially you use another element geometric matrix (which I have already used to include the effects of axial load - see my third post of this thread). The matrix is of the form in the attachment to this post. In the attachment, k_θ is the stiffness of the moment-rotation response, B is the breadth of the pile and l is the length of the element.

Reference:
DAVIDSON, H. L., CASS, P. G., KHILJI, K. H. and MCQUADE, P. V., 1982. Laterally Loaded Drilled Pier Research, Volume 1: Design Methodology. Monroeville: GAI Consultants Inc.

Contacting the Electric Power Research Institute (whose address was 3412 Hillview Avenue, Palo Alto, California 94304 in 1982) will result in an auto-reply with a link to their portfolio. The above reference can be accessed from there.

 

[gravityandinertia]

Sorry for misunderstanding, I see what you are saying. In my past experience, I've never seen a mass matrix implemented like that and it seems to have two problems for me:

1. If you are combining it with your stiffness matrix, you have inconsistent units
2. If you aren't combining with your stiffness matrix, you now must increase the size of all matrices involved to make the addition of the compatible with each other

Maybe we're perhaps speaking two different languages, but trying to say the same thing. I've always seen what you are calling the mass matrix called the foundational stiffness matrix, not mass. For the element equation you have shown above Mass*displacement =/= force, so you have inconsistent units there if this isn't a stiffness matrix of some sort.

I'm also not sure about your governing equation. It seems to me, if I'm thinking correctly on you p-y springs and they are foundational springs, then their stiffness should be multiplied by w instead of d²w/dx² since it is related to the amount the pile displaces. Additionally, I would have to spend more time with it, but the axial term looks correct for a bar that is loaded only axially, but looks incorrect for a beam in bending.

 

[LRJ]

The units of the mass matrix are consistent with the element stiffness matrix since the units of EI (kN/m²) are the same as k_p-y ( p is in kN/m; y is in m; therefore: p/y = kN/m/m = kN/m²).

I think this is almost certainly a case of different engineering fields having the same names for slightly different things. Indeed that was what I was alluding to before. Indeed mass doesn't actually come into the [m_M] matrix at all for this sort of analysis, so I think the name was borrowed.

You are also quite right that the equation I gave in my third post in this thread should've actually been:

EI.d⁴w/dx⁴ + P.d²w/dx² + k_p-y.dw/dx = q

I have corrected the post above to reflect this. Thanks for pointing it out - I had implemented it as above, though I didn't write it down here correctly.

The axial term is for axial loading of an element. Is there a different term for axial + bending? Is this where I need to upgrade it to a Timoshenko beam (i.e. to account for shear deformations)? That is on my to-do list. Or is there some other matrix that is used instead?

 

[gravityandinertia]

This is a good discussion so far. See my added thoughts below.

EI has units of kN-m² not kN/m².

E = kN/m²
I = m⁴
EI = kN-m²

I also don't understand your division by y on the mass matrix though it may be my lack of understanding of p-y springs. It may also be not understanding the terms you have contained in your y vector. Going back to your governing equation w(x) has units of m dw/dx is unitless since it's an angle. Units look something like this for derivatives of w(x):

w(x) = m
dw/dx = unitless
d²w/dx² = 1/m
d³w/dx³ = 1/m²
d⁴w/dx⁴ = 1/m³

From what I see in your governing equation you have the first term having units of kN/m, the second term, if P = AE/L, has units of kN/m² and the third term has units of kN/m² based on what you said about the k-y springs.

As for the timoshenko beam, it's been a while since I looked at it, and it's not as easy to remember as classical beam theory, so I would have to revisit that, which I can't do at the moment.

 

[LRJ]

Sorry, you're quite right about the units of EI. My mistake!

The governing equation and matrices are published in a number of places, so please don't rely on what I've put in here - there's every chance I've made other typos.