0
$\begingroup$

I have developed an Octave script to solve the nonlinear Euler-Bernoulli beam equations with linearized von Karman-strains, i.e. higher-order terms are dropped. The simulation results agree with analytical results in small and moderately large displacements in static analysis. However, when doing a dynamic simulation with large displacements I get incorrect results.

I used this site as a reference for assembling the stiffness matrix and the tangent stiffness matrix. I should note that I think the site has a typo for the stiffness coefficients $K_{21}$. I think the 0.5 should be dropped from the $K_{21}$ expression and also for the tangent stiffness matrix there should be $T_{21}=K_{21}$ instead of $T_{21}=2K_{21}$ (that is, if you don't symmetrize the matrices).

Now as the static linear analysis seems to work fine, I tried a dynamic nonlinear analysis. For this I used a simple pendulum with one beam element (length L=1m) starting from rest and a horizontal position. For the first few time steps everything seems fine until the beam element starts elongating (see figures below). At t=0.5 the length of the beam is almost 1.5m, which should not be possible. The internal force vector is initially pointing vertically upward but starts oscillating along with the inertial force vector and is never really aligned with the length of the beam as you would excpect for a tension force. Also if I do a static analysis and rotate the beam 90 degrees with Dirichlet BCs, I get large moments acting on the beam even though we are simulating rigid body rotation.

Pendulum at t=0.32s, Pendulum at t=0.5s

So my question is, is it possible to simulate large rotations in the total Lagrangian description (assuming small strains) and using the von Karman strains? Or is there some other reason why the simulation is not working? Or is my time integration scheme not suitable for this problem?

For the dynamic analysis I used the Newmark method for time integration (trapezoid rule), with boundary conditions set to u(0)=0 and w(0)=0. The stiffness matrices are integrated with gauss quadrature with reduced integration for the nonlinear terms, i.e. $K_{12}$, $K_{21}$, $T_{12}$, $T_{21}$ and $K_{22}^{NL}$ and full integration for the linear terms.

Improve this question
$\endgroup$
3
  • $\begingroup$ I have implemented a version of the geometrically exact beam in MATLAB/Octave. This follows the formulation given in the FEM book by Zienkiewicz (vol2). You can check it at this Github repository $\endgroup$
    – Chenna K
    Jan 9, 2022 at 16:43
  • $\begingroup$ I have used this code for the examples presented in paper1. $\endgroup$
    – Chenna K
    Jan 9, 2022 at 16:48
  • $\begingroup$ Thanks, thats a great paper. Will have to look closer into it! $\endgroup$
    – Tepa
    Jan 9, 2022 at 20:17

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.