Hi all, I am simulating a very simple cable element under traction (I also tested pure flexion) : one node is fixed, and a force is applied to the end node.
For the code attached, the elongation of the beam is 9.61 mm, while, analytically, it should be 12.7 mm. I did the same simulation with BeamEuler element, and the result is indeed 12.73 mm. With cable ANCF, even with lower force (divided by 10), the error is significant (1.11 mm vs 1.27 mm). The problem is the same in flexion and the error does not seem to be linear. When I increase the number of nodes in the builder, the error decreases but it remains significant (10.88 mm vs 12.7 mm with 1000 nodes instead of 1). I am using these elements wrong ? I thought they were suitable for large displacement and could replace Beam Euler if no twisting or shear were present ? Thanks a lot for the help, Solenne -- You received this message because you are subscribed to the Google Groups "ProjectChrono" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/projectchrono/bc489328-64ec-48de-b541-9d295dab0a83n%40googlegroups.com.
// ============================================================================= // PROJECT CHRONO - http://projectchrono.org // // Copyright (c) 2014 projectchrono.org // All rights reserved. // // Use of this source code is governed by a BSD-style license that can be found // in the LICENSE file at the top level of the distribution and at // http://projectchrono.org/license-chrono.txt. // // ============================================================================= // Authors: Dario Mangoni, Radu Serban // ============================================================================= // // FEA for 3D beams of 'cable' type (ANCF gradient-deficient beams) // uses the Chrono MKL module // // ============================================================================= #include "chrono/physics/ChSystemSMC.h" #include "chrono/timestepper/ChTimestepper.h" #include "chrono_pardisomkl/ChSolverPardisoMKL.h" #include "chrono/solver/ChIterativeSolverLS.h" #include "chrono_irrlicht/ChIrrApp.h" #include "FEAcables.h" int main(int argc, char* argv[]) { GetLog() << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << "\n\n"; // Select solver type (MKL or MINRES). ChSolver::Type solver_type = ChSolver::Type::MINRES; // Create a Chrono::Engine physical system ChSystemSMC my_system; my_system.Set_G_acc(0); // Create a mesh, that is a container for groups of elements and // their referenced nodes. auto my_mesh = chrono_types::make_shared<ChMesh>(); // Create the model (defined in FEAcables.h) // double beam_L = 0.1; // double beam_diameter = 0.015; double beam_L = 100; double beam_diameter = 1; double beam_E = 1e3; double beam_nu = 0.3; double beam_RayleighDamping = 0.001; double beam_density = 6.5e-9; ChVector<> force(100, 0, 0); // nodes coord ChVector<> A(0.0, 0.0, 0.0); ChVector<> B(beam_L, 0.0, 0.0); // Create a section, i.e. thickness and material properties // for beams. This will be shared among some beams. auto msection_cable = chrono_types::make_shared<ChBeamSectionCable>(); msection_cable->SetDiameter(beam_diameter); msection_cable->SetYoungModulus(beam_E); msection_cable->SetDensity(beam_density); msection_cable->SetBeamRaleyghDamping(beam_RayleighDamping); // This ChBuilderCableANCF helper object is very useful because it will // subdivide 'beams' into sequences of finite elements of beam type, ex. // one 'beam' could be made of 5 FEM elements of ChElementBeamANCF_3333 class. // If new nodes are needed, it will create them for you. ChBuilderCableANCF builder; // Now, simply use BuildBeam to create a beam from a point to another: builder.BuildBeam(my_mesh, // the mesh where to put the created nodes and elements msection_cable, // the ChBeamSectionCable to use for the ChElementBeamANCF_3333 elements 1, // the number of ChElementBeamANCF_3333 to create A, // the 'A' point in space (beginning of beam) B); // the 'B' point in space (end of beam) // After having used BuildBeam(), you can retrieve the nodes used for the beam, // For example say you want to fix both pos and dir of A end and apply a force to the B end: builder.GetLastBeamNodes().front()->SetFixed(true); builder.GetLastBeamNodes().back()->SetForce(force); // Remember to add the mesh to the system! my_system.Add(my_mesh); // Set solver and solver settings switch (solver_type) { case ChSolver::Type::PARDISO_MKL: { std::cout << "Using PardisoMKL solver" << std::endl; auto solver = chrono_types::make_shared<ChSolverPardisoMKL>(); my_system.SetSolver(solver); solver->UseSparsityPatternLearner(false); solver->LockSparsityPattern(false); solver->SetVerbose(false); break; } case ChSolver::Type::MINRES: { std::cout << "Using MINRES solver" << std::endl; auto solver = chrono_types::make_shared<ChSolverMINRES>(); my_system.SetSolver(solver); solver->SetMaxIterations(200); solver->SetTolerance(1e-12); solver->EnableDiagonalPreconditioner(true); solver->EnableWarmStart(true); // IMPORTANT for convergence when using EULER_IMPLICIT_LINEARIZED solver->SetVerbose(false); break; } default: { std::cout << "Solver type not supported." << std::endl; break; }} my_system.Update(); // Set integrator my_system.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED); // SIMULATION my_system.DoStaticNonlinear(20,true); // my_system.DoStaticNonlinear(); std::cout << " Cable Elongation = " << builder.GetLastBeamNodes().back()->GetPos().x() - B[0] << std::endl; return 0; }
