FEM modelling of a cantilever beam

[radou]

Due to complications I have with drawing and uploading pics at the moment, I'll simply describe the model instead of posting a pic; it's a cantilever beam (ignore the cross section and modulus of elasticity given in the pic, and the force P) of length L = 1, the coordinate axis x is along the beam, and the axis z is pointing downwards, with the origin at the fixed support on the left side of the beam. Imagine some continuous load q as a function of x distributed along the beam. This is the only load acting on the beam. The modulus of elasticity E and moment of inertia I are assumed to be constant along the beam. One needs to find the displacement w of the beam, i.e. solve the boundary problem:

$$EI \frac{d^4 w}{dx^4} = q$$ (1)

with boundary conditions

$$w(0) = 0$$ (vertical displacement at support equals zero)
$$w'(0) = 0$$ (slope at support equals zero - no rotation)
$$EI \frac{d^3 w}{dx^3} = 0$$ (shear force at x = 1 equals zero)
$$EI \frac{d^2 w}{dx^2} = 0$$ (torque at x = 1 equals zero).

This is where I need a push - I assume the first step is to find the variational formulation of the problem? Do I simply have to take some function satisfying the homogenous boundary condition and multiply equation (1) with it, integrate, and try to obtain a variational formulation that way? Further on, I read about the Galerkin approximation method. I'm interested in how exactly to construct the basis functions for the discretization space of the domain [0, 1] of the problem.

Anyone with some experience in FEM modeling of such simple mechanical systems - please help. The book I'm working with isn't specific enough, and it skips some steps. I'd like to solve this problem step by step in a clear and precise manner. Further on, I'm interesting in creating a Mathematica program, i.e. implementing this idea into an algorithm.

 

[PerennialII]

... using the Galerkin method is never boring I'll say ... great example how FEM actually comes to life and what's behind it all (+I prefer it over the other methods for some "clarity" ... probably a subjective view but anyways). Few things ... are you interested in the way FEM (for a beam element) is constructed in a "general engineering/physicist sense" or how one might construct the variational formulation in the "mathematician's" a tad more abstract way? The 1st would be the way leading to FEM and analysis of a beam the way it's commonly done for example in all major software packages, while the 2nd takes the "non-standard" way and works with the specific problem from the beginning using perhaps non-standard basis/shape/trial functions, fulfillment of the BCs of the PDE and so forth (I'm just thinking the derivations are typically somewhat different ... probably thinking too much here ).

Anyways, if following the Galerkin method to find the variational formulation the first steps would be writing of the residual eq/weak form with the weight function and all followed by the integration by parts routine. Pretty much how you've outlined it over there (take the PDE, multiply with an 'admissable' ("suitable") weight function and then proceed by integration by parts). Naturally what one faces pretty soon is the need to decide on what sort of an element one wishes to develop (the weight function). You're after the 'usual' Euler-Bernoulli beam element?

 

[radou]

Anyways, if following the Galerkin method to find the variational formulation the first steps would be writing of the residual eq/weak form with the weight function and all followed by the integration by parts routine.

OK, I'll do so and post my work here as soon as possible. I'd appreciate if you'd take a look at it then.

 

[radou]

So, let's take some function $v \in C^1([0,1])$ such that v(0) = 0, and multiply equation (1) with it. One easily obtains:
$$EI \int_{0}^1 \frac{d^4w(x)}{dx^4}v(x)dx = \int_{0}^1 q(x)v(x)dx$$
Now, after integrating by parts, I have
$$-EI \int_{0}^1 \frac{d^3w(x)}{dx^3}v'(x)dx = \int_{0}^1 q(x)v(x)dx$$,
but this seems to be wrong, since it's not what my book says. In the process of partial integration, I set
$$da = \frac{d^4w(x)}{dx^4} dx$$ and $$b = v$$
and used the partial integration formula $$ab |_{0}^1 - \int_{0}^1 a db$$.. (The notation isn't specially practical, since the v is already reserved for the weighting function, so I had to use a's and b's)

 

[PerennialII]

Looks good to me ... when you do the 2nd round of integration by parts does it look similar then (now you've moved one $\frac{d}{dx}$ to $v$, after "evening up the differentials" to get 2 on both v and w the weak form should be complete)? I'm wondering if that is the step your book skips ... ? Sure there is still the constant term from integration by parts which can keep there if don't want to introduce bcs "too early".

 

[radou]

Looks good to me ... when you do the 2nd round of integration by parts does it look similar then (now you've moved one $\frac{d}{dx}$ to $v$, after "evening up the differentials" to get 2 on both v and w the weak form should be complete)? I'm wondering if that is the step your book skips ... ? Sure there is still the constant term from integration by parts which can keep there if don't want to introduce bcs "too early".

Thanks for the hint, I just realized that I forgot to integrate by parts for the 2nd time!

So, we have (after integrating by parts of the left-hand side again):

$$EI w''(0)v'(0) + EI \int_{0}^1 w''(x)v''(x)dx = \int_{0}^1 q(x)v(x)dx$$

which represents the weak form of the problem, i.e. one needs to find, given the function q and the boundary conditions, the function w such that, for all functions v the equation above holds.

Now, the next step would be to exactly define the functional spaces we're dealing with here, right? Eg define where w and v are from, and formulate the problem in terms of finite dimensional functional spaces (i.e. subspaces of the original ones).

Perhaps it's more appropriate to switch to a different notation at this point, i.e. use

[tex]a(w, v) = EI \int_{0}^1 w''(x)v''(x)dx[/itex], and
$$(q, v) = \int_{0}^1 q(x)v(x)dx$$, so the weak form can be written as

$$a(w, v) = (q, v) - EI w''(0)v'(0)$$. (1')

Further on, for simplicity, could I assume EI = 1? Or perhaps that doesn't really matter?

 

[radou]

An addition to the previous post:

Let $$V = \left\{ v | \int_{0}^1 v'(x)^2 dx < \infty, v(0) = 0 \right\}$$.

So, the variational form of the problem would now be: for the given function q and boundary conditions, find w from V, such that for all v from V, the variational equation (1') from the post above holds. Am I correct on this one?

 

[PerennialII]

looks solid ... the constant term (which would then be substituted with the range of integration) should be something like:
$$\left(v \frac{d}{dx} EI \frac{d^{2}w}{dx^{2}}-\frac{dv}{dx} EI \frac{d^{2}w}{dx^{2}} \right)$$

... or well it could be if I didn't screw something up while doing a quick check.

Yeah, using the notation $a(v,w)$ and $(q,v)$ (or $l(v)$) is pretty common for the variational form, and a consistently accepted one. The defs for $v$ and $w$ could've been made already a tad earlier but yeah, and "forgetting" EI doesn't really matter much at this point. Having a twice differentiable $v$ will make it sound ($v \in H^{2}$) ... i.e. adding $\frac{d^{2} v}{d x^{2}}$ to V and the subspaces ought to be ok. Yep, the variational form even sounds ok.

And then off to interpolation.

 

[radou]

OK, so the next step is to create the Galerkin formulation of the problem.

So, we take a specific finite dimensional subspace V' of V, and write the variational equation (1') "in terms of it", e.g. we replace the elements from V with elements from V'.

Now, the problem is, how to choose this space V', and, further on, which basis to choose?

 

[PerennialII]

Yeah, this step has quite a bit of leeway in it. One can choose pretty much whatever and the variational formulation will yield a result of some sort (assuming the chosen v can come up with the bcs and belongs to V in general). So the decision is whether to choose a polynomial basis for a beam element (the number of dofs will then spec the type of polynomial) which will yield the 4 dof E-B beam element, or do something "exotic" and pick for example a harmonic basis or so. When developing a "general" element, the 1st is usually the route taken, although in reality this step is pretty free (illustrated well by math FEM texts, which compared to engineering ones are 'all over the place' at this point ).

 

[radou]

(illustrated well by math FEM texts, which compared to engineering ones are 'all over the place' at this point ).

Yeah, I know. That's why I use both. But the "math-ones", unfortunately, can become fairly complicated, and require additional efforts.

So, I could set V to be the space of all polynomials of degree less or equal to 3, and chose the standard basis {1, t, t^2, t^3}?

At this point I should look a bit more into the whole thing, but it seems I should now switch to the "element-point of view" and "solve" the variational equation for a generic element, and then merge the results into a global matrix with the help of transformation formulas between local and global coordinates.

This conception sounds beautiful and intuitively simple in theory, but obviously I'm pretty stuck right now.

 

[PerennialII]

Yeah, doing the polynomial general one is a good start. After that can create 'special' elements which solve specific problems with a single element with exceptional precision compared to general ones.

Using a general polynomial basic, 4 dofs -- the deflections and rotations on two nodes, and hence, $\lbrace 1, t, t^{2}, t^{3} \rbrace$ will do good.

Then it's about moving to work within a single element and using the polynomial basis to express (just thinking it'll 1st be done within an element, and the whole system is then when integrating over all elements)
$$w^{e} \left(x \right)=\sum_{i=1}^{4} u_{i}^{e} \phi_{i}^{e}$$
(used e to denote within element, i dof and u a generalized degree of freedom) for that element. The above has to also equal the polynomial expression for $w \left( x \right)$, and setting the polynomial expression equal to nodal values $u^{e}$ at nodal points will yield a linear system from where can solve the shape functions, $\phi_{i}^{e}$. Actually you don't need the expression $w^{e}$ for anything at this point, but after solving the linear system for the coefficients of the polynomial as a function of nodal dof values, that's the form that'll come out. Should get Hermite cubic splines out as a result. Making any sense?

 

[radou]

OK, so, essentially, this is the point: the variational approximation states that, for some finite dimensional subspace V' of V, we have to find w' from V', such that the variational equation (1') holds for every vector v from V' (of course, in the variational equation we substitute w with w'), i.e.

$$a(w', v) = (q, v) - EI w''(0)q(0)$$ (2)

(Note that I didn't replace the w'' in the constant term on the rhs with anything, I'm not sure if this is correct.)

Now since w' is from V' (let's say dim V' = n), and V' has a basis {$\phi_{1}, \dots, \phi_{n}$}, we can write w' = $\sum_{j=1}^n u_{j} \phi_{j}$. Now, the point is (at least this is how a Lemma about the existence of the solution of the variational approximation in general is proven in one of my texts), since (2) must hold for every v, we may "plug in" every basis element $\phi_{i}$ into (2), (along with w' written in terms of the basis vectors, and the fact that a(',') is a bilinear form), and hence we arrive at a system of n linear equations with n unknowns, i.e. the coefficients $u_{j}$.

More precisely, for v = $\phi_{i}$, and w' = $\sum_{j=1}^n u_{j} \phi_{j}$, (2) becomes:

$$\sum_{j = 1}^n u_{j} a(\phi_{j}, \phi_{i}) = (q, \phi_{i}) - EI w''(0)q(0)$$ (3)

for i = 1, ... n.

This is the system I was talking about, and, of course, now it is convenient to switch to matrix form and talk about the stiffness matrix, etc.

My first questions is: if I took V' to be (as mentioned earlier) the space of all polynomials of degree less or equal to 3 with the basis {1, t, t^2, t^3}, and if I went through the upper calculation, what exactly would I arrive at?

Second question: how many elements should there be, i.e. how many nodes? Should I actually write equation (3) for every element and include boundary conditions?

The problem is, at this state, practical and operative aspects are quite foggy to me.

 

[PerennialII]

...

This is the system I was talking about, and, of course, now it is convenient to switch to matrix form and talk about the stiffness matrix, etc.

My first questions is: if I took V' to be (as mentioned earlier) the space of all polynomials of degree less or equal to 3 with the basis {1, t, t^2, t^3}, and if I went through the upper calculation, what exactly would I arrive at?

You'd arrive at
$$w^{e} \left(x \right)=\sum_{i=1}^{4} u_{i}^{e} \phi_{i}^{e}$$

for w within an element (where $\phi_{i}$ are the result of the chosen basis), which will then be treated element-by-element in the variational formulation, $a \left(w, v \right) = \left(q, v \right)$, when the integration is carried out over the whole solution domain (=all elements).

Second question: how many elements should there be, i.e. how many nodes? Should I actually write equation (3) for every element and include boundary conditions?

At this point boundary conditions, when developing a general element, aren't required (sure could include them already, but then it wouldn't be general. One can keep all the dofs until using the element just for the sake of developing a general element for unspecified specific values of dof.).

The number of nodes and dofs is "input" to our selection of the basis (-> $\ldots t^{3}$ being the highest term), for the Euler-Bernoulli beam element with the variational formulation at hand, we need twice differentiable v:s to "get it running" (remembering V). Why 4 dofs? The resulting interpolation polynomial is in V (and will satisfy thus continuity conditions) and satisfies the essential boundary conditions (Dirichlet for displacement based FEM) of the problem (w, dw/dx at nodes ... so that the whole idea of interpolation using nodal values with the polynomial can be done consistently within an element ... 4 conditions, 4 parameter polynomial). (and this way can also nicely compute shear forces for the element btw, one d/dx more)

(extra* : sure in the end one can commit "variational crimes" and make up non-conforming elements for example, since there is really no reason to make "consistent" elements [well they usually work good due to their "solid" background... A reason ], but nothing is "really" forcing your hand in selecting these things... people have come up with great elements by breaking the "rules" and picking something inconsistent in terms of interpolation, fulfillment of bcs, continuity and so forth for particular elements. Getting bit ahead of myself :rofl: )