Commit f2b96606 authored by Stefan Pfeifer's avatar Stefan Pfeifer

Add displacement control

parent 35fc35be
......@@ -118,6 +118,14 @@ public:
}
// Todo: Should this be done via a View instead? Problem: Add vs set
double get_external_force(Dof dof)
{
if(dof.active)
return p(dof.index);
else
return 0.0;
}
void set_external_force(Dof dof, double val)
{
if(dof.active)
......@@ -169,6 +177,85 @@ public:
throw std::runtime_error("Maximum number of iterations exceeded");
}
void solve_statics_dc(Dof dof,
double target_displacement,
unsigned n_steps,
const std::function<void()>& callback)
{
v.setZero();
a.setZero();
unsigned max_iter = 50; // Todo: Magic number
double epsilon = 1e-8; // Todo: Magic number
if(!dof.active)
{
throw std::runtime_error("Displacement control not possible for a fixed DOF");
}
Eigen::VectorXd q(dofs());
Eigen::MatrixXd K(dofs(), dofs());
Eigen::VectorXd delta_q(dofs());
Eigen::VectorXd alpha(dofs());
Eigen::VectorXd beta(dofs());
Eigen::VectorXd e = Eigen::VectorXd::Zero(dofs());
e(dof.index) = 1.0;
// Todo: Think of a better control flow
auto solve_equilibrium = [&](double displacement)
{
unsigned iteration = 0;
while(iteration < max_iter)
{
if(iteration == max_iter - 1)
{
throw std::runtime_error("Maximum number of iterations exceeded");
}
// Calculate stiffness matrix and out of balance loads
get_internal_forces(q);
get_tangent_stiffness(K);
delta_q = p - q;
// Todo: Mode decomposition out of the main loop
Eigen::LDLT<Eigen::MatrixXd> stiffness_dec(K);
if(stiffness_dec.info() != Eigen::Success)
{
throw std::runtime_error("Decomposition of the stiffness matrix failed");
}
alpha = stiffness_dec.solve(delta_q);
beta = stiffness_dec.solve(e);
double delta_f = (displacement - u(dof.index) - alpha(dof.index))/beta(dof.index);
u += alpha + beta*delta_f;
p(dof.index) += delta_f;
update_element_states();
if((alpha + delta_f*beta).norm()/double(dofs()) < epsilon) // Todo: Better convergence criterion
{
callback();
return iteration;
}
++iteration;
}
callback();
return iteration;
};
double displacement = u(dof.index);
double initial_displacement = displacement;
for(unsigned i = 0; i < n_steps; ++i)
{
displacement += initial_displacement + (target_displacement - initial_displacement)*double(i)/double(n_steps - 1);
solve_equilibrium(displacement);
}
}
void solve_dynamics(double timestep_factor, const std::function<bool()>& callback)
{
VectorXd M(dofs());
......
......@@ -71,3 +71,38 @@ TEST_CASE("Small deformation bar truss")
REQUIRE(std::abs(s_numeric - s_analytic) < 1e-6);
}
// Todo: Why does the displacement control not allow passing the point 0.5*H?
// Read section on displacement control in 'Nonlinear Finite Element Analysis of Solids and Structures (René De Borst,Mike A. Crisfield,Joris J. C.)
TEST_CASE("Large deformation bar truss")
{
double H = 1.0;
double EA = 10000.0;
double F = 1500.0;
System system;
Node node01 = system.create_node({{ 0.0, 0.0, 0.0}}, {{false, false, false}});
Node node02 = system.create_node({{ H, H, 0.0}}, {{false, true, false}});
Node node03 = system.create_node({{2.0*H, 0.0, 0.0}}, {{false, false, false}});
BarElement element_01_02(node01, node02, M_SQRT2*H, EA, 0.0, 0.0);
BarElement element_02_03(node02, node03, M_SQRT2*H, EA, 0.0, 0.0);
system.add_element(element_01_02);
system.add_element(element_02_03);
// Solve for static equilibrium and get vertical deflection of node 02
system.solve_statics_dc(node02.y, 0.6*H, 100, [&]()
{
double s = system.get_u()(node02.y);
double f_num = system.get_external_force(node02.y);
double L0 = M_SQRT2*H;
double L = std::hypot(H, s);
double f_ref = 2.0*EA*s*(L - L0)/(L*L0);
double error = std::abs(f_num - f_ref); // Absolute error because f_ref = 0 initially
REQUIRE(error < 1e-9);
});
}
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