projectchrono/chrono

FEA mesh and rigid body contacts, there is a difference between the theoretical solution and the results.

695729709 opened this issue · 1 comments

Dear all,
I am solving the problem of a sliding block sliding down a slope under the action of gravity. But when I tried to solve the small-scale slider model, I encountered a problem where there was a significant difference between the theoretical solution and the results.
The FEA mesh is C3D4, which from abaqus. And the rigid plate is built by chrono which rotate 30 degrees along z-axis.
First,the length, width, and height of this slider model are 1.4m,0.52m and1.4m respectively. The calculation results of acceleration of node are as follows: ax is -2.0427, ay is -1.1785, and az is 1.2e-5. The theoretical solution are as follows: ax is -2.0385,ay is -1.1768, az is zero.
Compare two results, acceleration of node is quite similar to the theoretical solution when stable contacts occured.
And then ,I tried sclae the FEA mesh to one tenth.But the calculation results of acceleration acceleration of node when stable contacts occured are as follows: ax is -1.3986, ay is -0.7985, and az is 0.0045. The difference is significant between the theoretical solution and the results.

And this is my code.
`// =============================================================================
//
// FEA contacts using non-smooth contacts.
//
// =============================================================================

#include "chrono/physics/ChSystemSMC.h"
#include "chrono/physics/ChSystemNSC.h"
#include "chrono/physics/ChBodyEasy.h"
#include "chrono/physics/ChLoadContainer.h"
#include "chrono/geometry/ChTriangleMeshConnected.h"
#include "chrono/solver/ChIterativeSolverLS.h"
#include "chrono_pardisomkl/ChSolverPardisoMKL.h"

#include "chrono/fea/ChMesh.h"
#include "chrono/fea/ChMeshFileLoader.h"
#include "chrono/fea/ChContactSurfaceMesh.h"
#include "chrono/fea/ChContactSurfaceNodeCloud.h"
#include "chrono/assets/ChVisualShapeFEA.h"
#include "chrono/fea/ChElementCableANCF.h"
#include "chrono/fea/ChBuilderBeam.h"
#include "chrono/solver/ChSolverADMM.h"

#include "chrono/fea/ChLinkPointFrame.h"

#include "chrono_irrlicht/ChVisualSystemIrrlicht.h"

using namespace chrono;
using namespace chrono::geometry;
using namespace chrono::fea;
using namespace chrono::irrlicht;

// read cube_hole mesh
// set link for node and body
void test_1()
{
GetLog() << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << "\n\n";

// Create a Chrono::Engine physical system
ChSystemNSC sys;

sys.SetNumThreads(ChOMP::GetNumProcs(), 0, 1);

//
// CREATE THE PHYSICAL SYSTEM
//

collision::ChCollisionModel::SetDefaultSuggestedEnvelope(0.0025);
collision::ChCollisionModel::SetDefaultSuggestedMargin(0.0025);

// Set default effective radius of curvature for all SCM contacts.
collision::ChCollisionInfo::SetDefaultEffectiveCurvatureRadius(1);

// Use this value for an outward additional layer around meshes, that can improve
// robustness of mesh-mesh collision detection (at the cost of having unnatural inflate effect)
double sphere_swept_thickness = 0.002;

// Create the surface material, containing information
// about friction etc.
// It is a NSC non-smooth contact material that we will assign to
// all surfaces that might generate contacts.

auto mysurfmaterial = chrono_types::make_shared<ChMaterialSurfaceNSC>();
mysurfmaterial->SetFriction(0.3f);
mysurfmaterial->SetRestitution(0);

// Create a floor:

auto ground = chrono_types::make_shared<ChBody>();
ground->SetBodyFixed(true);
sys.Add(ground);


// floor as a simple collision primitive:
auto mfloor = chrono_types::make_shared<ChBodyEasyBox>(2, 0.1, 2, 2700, true, true, mysurfmaterial);
mfloor->SetPos(ChVector<>(0, -0.1, 0));
mfloor->SetRot(Q_from_AngAxis(CH_C_PI / 6.0, VECT_Z));
mfloor->SetBodyFixed(true);
mfloor->GetVisualShape(0)->SetTexture(GetChronoDataFile("textures/concrete.jpg"));
sys.Add(mfloor);


// create cylinder body for link to node
double radius_cylinder = 0.25/10.0, height_cylinder = 0.4/10.0;
auto cylinder = chrono_types::make_shared<ChBodyEasyCylinder>(geometry::ChAxis::Z, radius_cylinder, height_cylinder, 1000);
cylinder->SetPos(ChVector<>(0, 0.8/10.0, 0));
cylinder->GetVisualShape(0)->SetColor(ChColor(0.5f, 0.5f, 0));
sys.Add(cylinder);

// Create a mesh. We will use it for tetrahedrons.

auto my_mesh = chrono_types::make_shared<ChMesh>();

// Create a continuum material, that must be assigned to each solid element in the mesh,
// and set its parameters

auto mmaterial = chrono_types::make_shared<ChContinuumElastic>();
mmaterial->Set_E(200e9);  // rubber 0.01e9, steel 200e9
mmaterial->Set_v(0.3);
mmaterial->Set_RayleighDampingK(0.01);
mmaterial->Set_density(1000);

//
// Example: a tetrahedral tire
//

std::map<std::string, std::vector<std::shared_ptr<ChNodeFEAbase>>> node_sets;
ChCoordsys<> cpos(ChVector<>(0, 0.8/10.0, 0), Q_from_AngAxis(0, VECT_Y));
ChMeshFileLoader::FromAbaqusFile(my_mesh, GetChronoDataFile("models/tractor_wheel/cube_hole.inp").c_str(),
                                 mmaterial, node_sets, cpos.pos,ChMatrix33<>(1.0/10.0));

// Create the contact surface(s).
// In this case it is a ChContactSurfaceMesh, that allows mesh-mesh collsions.

auto mcontactsurf = chrono_types::make_shared<ChContactSurfaceMesh>(mysurfmaterial);
my_mesh->AddContactSurface(mcontactsurf);
mcontactsurf->AddFacesFromBoundary(sphere_swept_thickness);  // do this after my_mesh->AddContactSurface

// Remember to add the mesh to the system!
sys.Add(my_mesh);

//
// set link
//

auto nodeset_sel = "NODE_CONNECT";
for (auto i = 0; i < node_sets.at(nodeset_sel).size(); ++i) {
    auto mlink = chrono_types::make_shared<ChLinkPointFrame>();
    mlink->Initialize(std::dynamic_pointer_cast<ChNodeFEAxyz>(node_sets[nodeset_sel][i]), cylinder);
    sys.Add(mlink);
}


//
// Optional...  visualization
//

// Visualization of the FEM mesh.
// This will automatically update a triangle mesh (a ChTriangleMeshShape
// asset that is internally managed) by setting  proper
// coordinates and vertex colors as in the FEM elements.
// Such triangle mesh can be rendered by Irrlicht or POVray or whatever
// postprocessor that can handle a colored ChTriangleMeshShape).
auto mvisualizemesh = chrono_types::make_shared<ChVisualShapeFEA>(my_mesh);
mvisualizemesh->SetFEMdataType(ChVisualShapeFEA::DataType::NODE_SPEED_NORM);
mvisualizemesh->SetColorscaleMinMax(0.0, 5.50);
mvisualizemesh->SetSmoothFaces(true);
my_mesh->AddVisualShapeFEA(mvisualizemesh);

auto mvisualizemeshcoll = chrono_types::make_shared<ChVisualShapeFEA>(my_mesh);
mvisualizemeshcoll->SetFEMdataType(ChVisualShapeFEA::DataType::CONTACTSURFACES);
mvisualizemeshcoll->SetWireframe(true);
mvisualizemeshcoll->SetDefaultMeshColor(ChColor(0.2, 0.2, 0.2));
my_mesh->AddVisualShapeFEA(mvisualizemeshcoll);



// Create the Irrlicht visualization system
auto vis = chrono_types::make_shared<ChVisualSystemIrrlicht>();
vis->AttachSystem(&sys);
vis->SetWindowSize(800, 600);
vis->SetWindowTitle("FEA contacts");
vis->Initialize();
vis->AddLogo();
vis->AddSkyBox();
vis->AddTypicalLights();
vis->AddCamera(ChVector<>(0.0, 0.0, 2.4));
vis->EnableContactDrawing(ContactsDrawMode::CONTACT_DISTANCES);

// SIMULATION LOOP

// Use the ADMM solver: it has the capability of handling both FEA and NSC!
auto solver = chrono_types::make_shared<ChSolverADMM>(chrono_types::make_shared<ChSolverPardisoMKL>());
solver->EnableWarmStart(true);
solver->SetMaxIterations(60);
solver->SetToleranceDual(1e-8);
solver->SetTolerancePrimal(1e-8);
solver->SetRho(1);
solver->SetStepAdjustPolicy(ChSolverADMM::AdmmStepType::BALANCED_UNSCALED);
sys.SetSolver(solver);

sys.SetSolverForceTolerance(1e-10);

auto node_output = std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(100));

while (vis->Run()) {
    vis->BeginScene();
    vis->Render();
    vis->RenderFrame(ChFrame<>(ground->GetCoord()));

    vis->EndScene();
    sys.DoStepDynamics(0.005);
    std::cout << "Time: " << sys.GetChTime() << std::endl;
    std::cout << "Ncontacts: " << sys.GetNcontacts() << std::endl;
    std::cout << "Pos: " << node_output->GetPos() << std::endl;
    std::cout << "Vel: " << node_output->GetPos_dt() << std::endl;
    std::cout << "Acc: " << node_output->GetPos_dtdt() << std::endl;
    std::cout << "CylinderPos: " << cylinder->GetPos() << std::endl;
    std::cout << "CylindeVel: " << cylinder->GetPos_dt() << std::endl;
    std::cout << "CylindeAcc: " << cylinder->GetPos_dtdt() << std::endl;
}

}

int main(int argc, char* argv[])
{

test_1();


return 0;

}
`

the results of the large slider is as follows:
large slider

the results of the small slider is as follows:
small slider