Nonlinear response of a cantilever beam

[aeroaero]

Hi all,

I am working on the geometrically nonlinear dynamic response of a cantilever beam using an in-house FEM code. I have a 16 m cantilever beam whose fundamental frequency is around 5 Hz.

In a particular simulation, I discretized the beam into 50 elements and tried to use a time-step of 0.0001 sec for a tip excitation of 50sin20t. I get some kind of a singular matrix error in my code. I tried the same excitation with the same time-step but with around 20 finite elements. It then worked. It did'nt work for 30 elements again. I may be exciting higher modes which can't be handled by my solver. I am not sure, though.

However, I am sure there is a rule as to how small my time-step could be for a particular discretization of my structure, just like how we use CFL number in the numerical simulation of fluid problems.

I am wondering if someone can point me to a reference which explains the relation between the FEM discretization and time-step to be taken for geometrically nonlinear dynamic simulations of the kind I described above. Also, it would help if someone can comment on my observations above.

Thanks in advance

Regards,
Aero^2

 

[GregLocock]

No references sorry, but practically 1e-4 time steps is ample in any non linear code or sim I have written, for a 5 Hz mode. Typically I would use a 1 ms step for a model with a frequency response of 30 Hz, and in theory that is excessive.






Cheers

Greg Locock

Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[Drej]

First of all, I would check your units. If you applying a 50sin20t fixed displacement you are forcing the cantilever at a single frequency (looking at your equation, it doesn't appear to be at 5 Hz though), so you should know whether or not you're exciting this mode in your model. Also, if your forcing amplitude is 50 units, and these are metres (you refer to the beam being 16m long), then you will obviously expect problems, even without forcing the beam at it's fundamental frequency (fixed displacement of 50m), unless it is over a very short period of time.


------------
See faq569-1083 for details on how to make best use of Eng-Tips.com

 

[aeroaero]

Hi Greg and Drej,

Thanks for responding so promptly.

Drej, let me correct myself. The forcing is 5sin30t and not 50sin20t. The problem does'nt seem to be the forcing frequency at all. It is with my finite-element discretization. There seems to be a connection between the number of beam finite elements I choose and the time-step of the simulation.

So, just to summarize the problem, the beam is of length 16m split into 50 finite elements, fundamental mode is at 5 Hz, the tip forcing function is 5sin30t, and the time-step I choose for this is 0.0001 which did'nt work. I get a singular matrix error for this choice of time-step. However, when I chose 10 finite elements, this time-step worked.

When I chose 0.001 as the time-step, even the 50 beam finite elements case worked. Basically, any time-step above 0.001 one worked for me for any number of finite elements. For some reason, I am unable to go below a certain time-step when I choose more number of finite elements. I am almost certain that there is no physics that I am violating here. Also, I am sure that my nonlinear beam solver is working fine.

So, please let me know what you think of this problem and the relation between time-step and the number of fem elements.

Thanks in advance

Regards,
Aero^2

 

[GregLocock]

So are you applying a force or a displacement? How quickly (ie at what elapsed time) do you get the singularity error message?

What damping have you used?




Cheers

Greg Locock

Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[GregLocock]

Oh, and what material are you using, and have you restrained it against out of plane motion?

Cheers

Greg Locock

Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[aeroaero]

Hi Greg,

I am applying a force at the tip as I mentioned (5sin30t). I get this singularity in my Newton-Raphson jacobian at the second time-step itself which is around 0.0002 s (when the time-step is 0.0001 s).


Thanks,
Satish

 

[aeroaero]

Hi,

The EI (flat-wise) is 4x10^6 N-m^2 and the mass per unit length of the beam is 0.769 kg/m. Again, as mentioned, this is a cantilever beam with a sinusoidal force at the tip. And, I have a diagonal stiffness matrix which uncouples one motion from the other.

It is definitely restrained out of plane. It has a cantilevered boundary condition on one end.

As mentioned before, the same set of things work for time-step 0.0001 and 10 finite elements and not for 0.0001 and 50 finite elements. A time-step of 0.001 and above works for any number of finite elements.


Thanks,
AeroAero

 

[aeroaero]

Oops.. forgot to mention that I do not have any structural damping in the model.

 

[GregLocock]

So in the real world you are applying a very high frequency impulse to an undamped system. Shear waves or longitudinal waves might become important.



Cheers

Greg Locock

Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[aeroaero]

Greg,

Thanks. So, by choosing such a small time-step, am I applying a high frequency impulse ?? Also, what would shear waves have to do with the finite element discretization ? Please clarify if possible.

Thanks,
Aero^2

 

[GregLocock]

You may start to excite modal behaviour in directions you aren't expecting, such as axially or shear. I don't know any more than that, but I was trying to think of possible exlanations for your model's behaviour.

To be honest I don't see how you can say "Also, I am sure that my nonlinear beam solver is working fine. ", but if it is working properly then there must be a physical solution.



Cheers

Greg Locock

Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[mastroit2002]

I agree with GregLocock: the finite element beam solver may be the kernel of the question.

I tried to figure out the kind of problem You are solving and the way chosen for solution: I think that You are solving a transient analysis of a cantilever beam whose fixed point undergoes from time 0 a specified dispacement function of time.

As I remember from bernulli like models (correct me if different) you are solving a problem of this form:

(EJ w,xx),xx= mu w,tt [x,t space and time dimension, EJ beam stiffness, w beam local displacement]

it follows (from dimension calculation) that first modal frequency should be near:

omega^2=EJ/(mu l^4) -> near 9 rad/s (1.42 Hz)

I wonder on what a kind of non linearity You are dealing with (EJ,m non constant or time dependent) that would give a very different value for first natural frequency.

" CFL number in the numerical simulation of fluid problems" exists even for these structural problems but is central only for pure explicit integration schemes.

It shouldn't affect implicit like solution methods (Newton-Raphson jacobian for time derivative).

The time step shouldn't be the problem but I think that the problem could be the way the initial condition is set: you may have a non finite at end spatial derivative under the presence of step of super imposed displacement.

Are you using any numerical tool like Matlab to perform this task? I still prefer explicit methods on implicit (too mathematically complex) methods [Matrix multiplication vs linear system solution].

With best regards,

mastro

 

[aeroaero]

Dear Mastro,

What I solve is a set of nonlinear beam equations derived based on no kinematic assumptions to the beam reference line as is done in our typical E-B/Timoshenko models.
You were a bit off in the calculation of the fundamental mode. There is a factor of 3.52 after you take the square root for omega in your formula.

I now realize that time-step should not be a problem for implicit methods. Did'nt think when I made the CFL number statement in my previous post. Thanks.

Now, I don't understand the following statement in your post:

"The time step shouldn't be the problem but I think that the problem could be the way the initial condition is set: you may have a non finite at end spatial derivative under the presence of step of super imposed displacement."

What initial condition are we talking about here. I have zero initial condition to my solution vector. And what is this non finite at end spatial derivative .. ?? Would you please elaborate this a bit ?

I use Fortran90 and not MATLAB or anything like that.


Thanks,
Aero^2

 

[mastroit2002]

>You were a bit off in the calculation of the fundamental >mode. There is a factor of 3.52 after you take the square >root for omega in your formula.

Probably you are right: that was only for my by hand calculation (that is interesting that simply using all beam stiffness and mass reported to the free side there should be an 73% of error).

I don't think you can trust me as a printed book when I was giving a meaning to your question.

From Bernulli simplified:
omega_n^2= 3 * EJ/M*l^3= 3 * EJ/Mu*l^4

The magic number 3.52 recalls me the solution of trascendental equation "cos(alfa)*cosh(alfa)+1=0" alfa^2=3.52 I've got from mechanics book.

omega_n^2= 3.52 * EJ/Mu*l^4

I don't know what a kind of method You are using: finite element methods necessarly make an assuption on a the spatial field function.

I claimed that "but I think that the problem could be the way the initial condition is set: you may have a non finite at end spatial derivative under the presence of step of super imposed displacement."

You answered that "I have zero initial condition to my solution vector".

The problem can follow due to time solution for non finite derivatives. When performing the simulation there could be a by zero division when you are starting the simulation due to the fact of having a beam end with a time imposed displacement.

Hw, number of point of discretization and frequency are connected: every doubling of points should increment four times the frequency not for stability but for resolution.

I don't think this could be really of help since there several topics I cannot figure out from Your problem: hw feel free to contact me for more insight.

Regards,

Mastro

 

[GBor]

Aero^2,

I just ran this on a commercial package with no problems. Fundamental frequency of around 4.95 Hz. Time stepping at 0.0001. I had to make up E and I independently to yield the proper stiffness and I discretized the forcing function, but it ran...I realize that doesn't help much, but perhaps your non-linear processor has its limitations.

GBor

 

[toog]

AeroAero,

Probably a dumb question, but are you using double precision in your Fortran code? Single precision FEM codes lose all significant digits for about 30+ degrees of freedom which seems to be about where you get to.

 

[aeroaero]

Thanks Mastro, GBor, and toog for your responses.

GBor, thanks for trying it with a commericial package. Are you trying with ANSYS/NASTRAN ? Just curious. Do you have the test-case parameters so that I can see it working here ? Please don't be obligated to do this.

I set-up a case in MSC.NASTRAN myself and I am yet to get the result. The case is still running. Hopefully, it gives junk at the end of the day so that I don't have to doubt my beam solver for anything.

Toog - I am using double precision (real(8) is what I use everywhere in the code). So, that's not the problem.


Thanks guys again. Will keep you ppl posted.

Regards,
Aero^2

 

[GBor]

I used Algor. I'm not in the office today (a wife and four daughters...Valentine's Day becomes a holiday!), but the parameters were something like:

Round beam 16 m long with a 1 meter radius (yes, I know this is ridiculous, but I was trying to do this quickly). I seem to recall this giving an I of approximately 0.74 (sorry, too lazy to recalculate). I set the material properties to give the EI that you required. I used a poisson's of 0.3 and calculated G for a linear material. I used a mass density that provided the weight load that you specified. broke the beam into 50 elements. Since you said you had constrained it to the in-plane directions, my boundary conditions were fixed at one end and no out-of-plane translations or rotations. Basically this would be for a beam in the xy plane: fixed on one end, no z translation, no x or y rotations.

I set up the forcing function by placing a unit load on the other end and exciting it with the sine wave that you specified.

I think that's about it.

GBor

 

[pja]

1) I am assuming you have check that you get the correct linear results..ie your eigenvalues converge to the correct answer as you increase the number of elements, etc
2) If you are sure about above then the fact that you get an answer with a smaller timestep tells me that you are having problems with an initial transient, as someone else mentioned.

However, this being said, your result should get better in this case as you add more elements not worse. I would check your nonlinear solver against some published results. There are plenty to chose from.