Commit 197346e7 authored by Stefan Pfeifer's avatar Stefan Pfeifer

Add beam element and tests

parent f2b96606
......@@ -9,6 +9,7 @@ HEADERS += src/gui/MainWindow.hpp \
src/fem/System.hpp \
src/fem/elements/Element.hpp \
src/fem/elements/BarElement.hpp \
src/fem/Node.hpp \
src/fem/elements/MassElement.hpp \
src/fem/elements/BeamElement.hpp \
src/fem/Node.hpp \
src/fem/View.hpp
......@@ -65,6 +65,15 @@ public:
elements.push_back(&element);
}
template<class ElementContainer>
void add_elements(ElementContainer& container)
{
for(auto it = container.begin(); it != container.end(); ++it)
{
add_element(*it);
}
}
size_t dofs() const
{
return u.size();
......@@ -132,6 +141,24 @@ public:
p(dof.index) = val;
}
// Euclidean distance between node_a and node_b
double get_distance(Node node_a, Node node_b)
{
double dx = get_u()(node_b.x) - get_u()(node_a.x);
double dy = get_u()(node_b.y) - get_u()(node_a.y);
return std::hypot(dx, dy);
}
// Angle between the x-axis and the line connecting node_a and node_b
double get_angle(Node node_a, Node node_b)
{
double dx = get_u()(node_b.x) - get_u()(node_a.x);
double dy = get_u()(node_b.y) - get_u()(node_a.y);
return std::atan2(dy, dx);
}
void solve_statics_lc()
{
v.setZero();
......@@ -148,7 +175,7 @@ public:
get_tangent_stiffness(K);
double norm_0 = (p - q).norm();
if(norm_0 == 0)
if(norm_0 == 0) // Todo: Introduce absolute tolerance for this
{
return; // System already in equilibrium
}
......@@ -164,14 +191,14 @@ public:
u += stiffness_dec.solve(p - q);
update_element_states();
get_internal_forces(q);
get_tangent_stiffness(K);
double norm_i = (p - q).norm();
if(norm_i/norm_0 < epsilon)
{
return;
}
get_internal_forces(q);
get_tangent_stiffness(K);
}
throw std::runtime_error("Maximum number of iterations exceeded");
......@@ -202,6 +229,8 @@ public:
Eigen::VectorXd e = Eigen::VectorXd::Zero(dofs());
e(dof.index) = 1.0;
Eigen::LDLT<Eigen::MatrixXd> stiffness_dec;
// Todo: Think of a better control flow
auto solve_equilibrium = [&](double displacement)
{
......@@ -218,8 +247,7 @@ public:
get_tangent_stiffness(K);
delta_q = p - q;
// Todo: Mode decomposition out of the main loop
Eigen::LDLT<Eigen::MatrixXd> stiffness_dec(K);
stiffness_dec.compute(K);
if(stiffness_dec.info() != Eigen::Success)
{
throw std::runtime_error("Decomposition of the stiffness matrix failed");
......@@ -246,12 +274,12 @@ public:
return iteration;
};
double displacement = u(dof.index);
double initial_displacement = displacement;
double init_displacement = u(dof.index);
double delta_displacement = target_displacement - init_displacement;
for(unsigned i = 0; i < n_steps; ++i)
{
displacement += initial_displacement + (target_displacement - initial_displacement)*double(i)/double(n_steps - 1);
double displacement = init_displacement + delta_displacement*double(i)/double(n_steps - 1);
solve_equilibrium(displacement);
}
}
......
#pragma once
#include "Element.hpp"
#include <cmath>
class BeamElement: public Element
{
private:
std::array<Dof, 6> dofs;
double phi_ref_0;
double phi_ref_1;
double rhoA;
double L;
double dx;
double dy;
Eigen::Matrix<double, 3, 6> J;
Eigen::Matrix<double, 3, 1> e;
Eigen::Matrix<double, 3, 3> C;
public:
BeamElement(Node nd0, Node nd1, double Cee, double Ckk, double Cek, double rhoA, double L)
: dofs{{nd0.x, nd0.y, nd0.phi, nd1.x, nd1.y, nd1.phi}},
phi_ref_0(0.0),
phi_ref_1(0.0),
rhoA(rhoA),
L(L)
{
C << Cee, -Cek, Cek,
-Cek, 4.0*Ckk, 2.0*Ckk,
Cek, 2.0*Ckk, 4.0*Ckk;
}
void set_reference_angles(double phi_ref_0, double phi_ref_1)
{
this->phi_ref_0 = phi_ref_0;
this->phi_ref_1 = phi_ref_1;
}
virtual void set_state(const VectorView<Dof> u, const VectorView<Dof> v)
{
dx = u(dofs[3]) - u(dofs[0]);
dy = u(dofs[4]) - u(dofs[1]);
double alpha = std::atan2(dy, dx);
// Elastic coordinates
double sin_e1 = std::sin(u(dofs[2]) + phi_ref_0 - alpha);
double cos_e1 = std::cos(u(dofs[2]) + phi_ref_0 - alpha);
double sin_e2 = std::sin(u(dofs[5]) + phi_ref_1 - alpha);
double cos_e2 = std::cos(u(dofs[5]) + phi_ref_1 - alpha);
e(0) = std::hypot(dx, dy) - L;
e(1) = std::atan(sin_e1/cos_e1); // Todo: Replace atan(sin(x)/cos(x)) with a cheaper function
e(2) = std::atan(sin_e2/cos_e2);
// Jacobian of the elastic coordinates
double a0 = std::pow(dx*dx + dy*dy, -0.5);
double a1 = 1.0/(dx*dx + dy*dy);
double j0 = a0*dx;
double j1 = a0*dy;
double j2 = a1*dx;
double j3 = a1*dy;
J << -j0, -j1, 0.0, j0, j1, 0.0,
-j3, j2, 1.0, j3, -j2, 0.0,
-j3, j2, 0.0, j3, -j2, 1.0;
}
virtual void get_masses(VectorView<Dof> M) const override
{
double alpha = 0.02; // Todo: Magic number
double m = 0.5*rhoA*L;
double I = alpha*rhoA*std::pow(L, 3);
M(dofs[0]) += m;
M(dofs[1]) += m;
M(dofs[2]) += I;
M(dofs[3]) += m;
M(dofs[4]) += m;
M(dofs[5]) += I;
}
virtual void get_internal_forces(VectorView<Dof> q) const override
{
q(dofs) += 1.0/L*J.transpose()*C*e;
}
virtual void get_tangent_stiffness(MatrixView<Dof> K) const override
{
// Todo: Introduce some more variables for common expressions
double a0 = std::pow(dx*dx + dy*dy, -0.5);
double a1 = 1.0/(dx*dx + dy*dy);
double a2 = std::pow(dx*dx + dy*dy, -1.5);
double a3 = std::pow(dx*dx + dy*dy, -2.0);
double b0 = a2*dx*dx - a0;
double b1 = a2*dy*dy - a0;
double b2 = a2*dx*dy;
double b3 = 2.0*a3*dx*dx - a1;
double b4 = 2.0*a3*dy*dy - a1;
double b5 = 2.0*a3*dx*dy;
Eigen::Matrix<double, 3, 6> dJ0;
dJ0 << -b0, -b2, 0.0, b0, b2, 0.0,
-b5, b3, 0.0, b5, -b3, 0.0,
-b5, b3, 0.0, b5, -b3, 0.0;
Eigen::Matrix<double, 3, 6> dJ1;
dJ1 << -b2, -b1, 0.0, b2, b1, 0.0,
-b4, b5, 0.0, b4, -b5, 0.0,
-b4, b5, 0.0, b4, -b5, 0.0;
Eigen::Matrix<double, 6, 6> Kn = Eigen::Matrix<double, 6, 6>::Zero();
Kn.col(0) = 1.0/L*dJ0.transpose()*C*e;
Kn.col(1) = 1.0/L*dJ1.transpose()*C*e;
Kn.col(3) = -1.0/L*dJ0.transpose()*C*e;
Kn.col(4) = -1.0/L*dJ1.transpose()*C*e;
K(dofs) += Kn + 1.0/L*J.transpose()*C*J;
}
virtual double get_potential_energy() const override
{
return 0.5/L*e.transpose()*C*e;
}
};
......@@ -18,7 +18,7 @@ public:
}
virtual void set_state(const VectorView<Dof> u, const VectorView<Dof> v)
virtual void set_state(const VectorView<Dof>, const VectorView<Dof>)
{
}
......@@ -27,15 +27,15 @@ public:
{
M(node.x) += m;
M(node.y) += m;
M(node.phi) += m;
M(node.phi) += I;
}
virtual void get_internal_forces(VectorView<Dof> q) const override
virtual void get_internal_forces(VectorView<Dof>) const override
{
}
virtual void get_tangent_stiffness(MatrixView<Dof> K) const override
virtual void get_tangent_stiffness(MatrixView<Dof>) const override
{
}
......
......@@ -8,5 +8,6 @@ CONFIG -= qt
SOURCES += tests/main.cpp \
tests/harmonic_oscillator.cpp \
tests/tangent_stiffness.cpp \
tests/bar_trusses.cpp
tests/bar_trusses.cpp \
tests/large_deformation_beams.cpp
#include "../src/fem/System.hpp"
#include "../src/fem/elements/BeamElement.hpp"
#include <catch.hpp>
#include <vector>
#include <iostream>
// "Large deformation of a cantilever beam subjected to a vertical tip load" as described in [1]
// [1] On the correct representation of bending and axial deformation in the absolute nodal coordinate formulation with an elastic line approach
// Johannes Gerstmayr, Hans Irschik, Journal of Sound and Vibration 318 (2008) 461-487
TEST_CASE("Large deformation cantilever")
{
unsigned N = 15; // Number of elements
double L = 2.0;
double b = 0.1;
double h = 0.1;
double I = b*h*h*h/12.0;
double A = b*h;
double E = 2.07e11;
double F0 = 3.0*E*I/(L*L);
System system;
std::vector<Node> nodes;
std::vector<BeamElement> elements;
// Create nodes
for(unsigned i = 0; i < N+1; ++i)
{
bool active = (i != 0);
nodes.push_back(system.create_node({{double(i)/double(N)*L, 0.0, 0.0}}, {{active, active, active}}));
}
// Create elements
for(unsigned i = 0; i < N; ++i)
{
elements.push_back(BeamElement(nodes[i], nodes[i+1], E*A, E*I, 0.0, 0.0, L/double(N)));
}
system.add_elements(elements);
system.set_external_force(nodes[N].y, F0);
system.solve_statics_lc();
// Tip displacements
double ux_num = L - system.get_u()(nodes[N].x);
double uy_num = system.get_u()(nodes[N].y);
// Reference displacements according to [1]
double ux_ref = 0.5089228704;
double uy_ref = 1.2083981311;
double error_x = std::abs((ux_num - ux_ref)/ux_ref);
double error_y = std::abs((uy_num - uy_ref)/uy_ref);
// Todo: Why is the precision not better?
REQUIRE(error_x < 1.1e-3);
REQUIRE(error_y < 1.1e-3);
}
// Static test "Bending of a pre-curved beam into a full circle" from [1]
// [1] On the correct representation of bending and axial deformation in the absolute nodal coordinate formulation with an elastic line approach
// Johannes Gerstmayr, Hans Irschik. Journal of Sound and Vibration 318 (2008) 461-487
TEST_CASE("Large deformation circular beam")
{
unsigned N = 20; // Number of elements
double R = 1.0;
double EA = 5.0e8;
double EI = 10.0e5;
System system;
std::vector<Node> nodes;
std::vector<BeamElement> elements;
// Create nodes
for(unsigned i = 0; i < N+1; ++i)
{
bool active = (i != 0);
double phi = double(i)/double(N)*M_PI;
nodes.push_back(system.create_node({{R*std::sin(phi), R*(std::cos(phi) - 1.0), 0.0}}, {{active, active, active}}));
}
// Create elements
for(unsigned i = 0; i < N; ++i)
{
double dist = system.get_distance(nodes[i], nodes[i+1]);
double angle = system.get_angle(nodes[i], nodes[i+1]);
BeamElement element(nodes[i], nodes[i+1], EA, EI, 0.0, 0.0, dist);
element.set_reference_angles(angle - system.get_u()(nodes[i].phi),
angle - system.get_u()(nodes[i+1].phi));
elements.push_back(element);
}
system.add_elements(elements);
system.solve_statics_dc(nodes[N].phi, -M_PI, 15, [&]()
{
//std::cout << "phi = " << system.get_u()(nodes[N].phi) << "\n";
});
double M_num = -system.get_external_force(nodes[N].phi);
double M_ref = EI*R;
double error_M = std::abs((M_num - M_ref)/M_ref);
double error_x = std::abs(system.get_u()(nodes[N].x));
double error_y = std::abs(system.get_u()(nodes[N].y));
// Todo: Why is the error in the torque so much higher than the displacement error?
REQUIRE(error_M < 1.1e-3);
REQUIRE(error_x < 1.0e-14);
REQUIRE(error_y < 1.0e-14);
}
#include "../src/fem/System.hpp"
#include "../src/fem/elements/BarElement.hpp"
#include "../src/fem/elements/BeamElement.hpp"
#include <catch.hpp>
#include <iostream>
......@@ -46,3 +47,35 @@ TEST_CASE("Tangent stiffness matrix: BarElement")
check_at(-1.5, -1.5, 0.0, 0.0);
}
TEST_CASE("Tangent stiffness matrix: BeamElement")
{
auto check_at = [](double x0, double y0, double phi0, double x1, double y1, double phi1)
{
double L = 1.0;
double EA = 100.0;
double EI = 10.0;
System system;
Node node0 = system.create_node({{x0, y0, phi0}}, {{true, true, true}});
Node node1 = system.create_node({{x1, y1, phi1}}, {{true, true, true}});
BeamElement element01(node0, node1, EA, EI, 0.0, 0.0, L);
system.add_element(element01);
system.update_element_states();
system.update_element_states();
system.update_element_states();
check_stiffness_matrix(system);
};
check_at(0.0, 0.0, 0.0, 1.0, 0.0, 0.0);
check_at( 0.0, -6.0, 1.0, -7.0, -2.0, 3.0);
check_at( 4.0, -7.0, -5.0, -9.0, 6.0, -6.0);
check_at( 9.0, -4.0, 0.0, 4.0, 2.0, -7.0);
check_at( 2.0, 0.0, 5.0, 10.0, -1.0, 8.0);
check_at(-6.0, 5.0, 4.0, -6.0, -6.0, -5.0);
check_at(-4.0, -6.0, -4.0, 5.0, 0.0, 5.0);
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment