Perform thermal analysis (demo_FEA_thermal.cpp)

Tutorial that teaches how to use the FEA module to compute the scalar temperature field in a solid (Poisson problem) given boundary conditions.

// =============================================================================
// 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: Alessandro Tasora
// =============================================================================
//
// FEA visualization using Irrlicht
//
// =============================================================================
#include "chrono/physics/ChSystemSMC.h"
#include "chrono/solver/ChIterativeSolverLS.h"
#include "chrono/fea/ChElementBar.h"
#include "chrono/fea/ChElementTetraCorot_4.h"
#include "chrono/fea/ChElementTetraCorot_10.h"
#include "chrono/fea/ChElementHexaCorot_8.h"
#include "chrono/fea/ChElementHexaCorot_20.h"
#include "chrono/fea/ChContinuumThermal.h"
#include "chrono/fea/ChContinuumElectrostatics.h"
#include "chrono/fea/ChNodeFEAxyzP.h"
#include "chrono/fea/ChMesh.h"
#include "chrono/fea/ChMeshFileLoader.h"
#include "chrono/fea/ChLinkNodeFrame.h"
#include "chrono/assets/ChVisualShapeFEA.h"
#include "chrono_irrlicht/ChVisualSystemIrrlicht.h"
// Remember to use the namespace 'chrono' because all classes
// of Chrono::Engine belong to this namespace and its children...
using namespace chrono;
using namespace chrono::fea;
using namespace chrono::irrlicht;
int main(int argc, char* argv[]) {
std::cout << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << std::endl;
// Create a Chrono::Engine physical system
// Create a mesh, that is a container for groups
// of elements and their referenced nodes.
auto my_mesh = chrono_types::make_shared<ChMesh>();
// Create a material, that must be assigned to each element,
// and set its parameters
auto mmaterial = chrono_types::make_shared<ChContinuumThermal>();
mmaterial->SetSpecificHeatCapacity(2);
mmaterial->SetThermalConductivity(200);
//
// Add some TETAHEDRONS:
//
// Load a .node file and a .ele file from disk, defining a complicate tetrahedron mesh.
// This is much easier than creating all nodes and elements via C++ programming.
// You can generate these files using the TetGen tool.
try {
ChMeshFileLoader::FromTetGenFile(my_mesh, GetChronoDataFile("fea/beam.node").c_str(),
GetChronoDataFile("fea/beam.ele").c_str(), mmaterial);
} catch (std::exception myerr) {
std::cerr << myerr.what() << std::endl;
return 0;
}
for (unsigned int inode = 0; inode < my_mesh->GetNumNodes(); ++inode) {
if (auto mnode = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(inode))) {
mnode->SetPos(mnode->GetPos() * ChVector3d(3, 1, 3));
}
}
//
// Set some BOUNDARY CONDITIONS on nodes:
//
// Impose load on the 180th node
auto mnode3 = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(180));
mnode3->SetLoad(20); // thermal load: heat flux [W] into node
// Impose field on two top nodes (remember the SetFixed(true); )
auto mnode1 = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(my_mesh->GetNumNodes() - 1));
mnode1->SetFixed(true);
mnode1->SetFieldVal(0.5); // field: temperature [K]
auto mnode2 = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(my_mesh->GetNumNodes() - 2));
mnode2->SetFixed(true);
mnode2->SetFieldVal(0.5); // field: temperature [K]
// Impose field on the base points:
for (unsigned int inode = 0; inode < my_mesh->GetNumNodes(); ++inode) {
if (auto mnode = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(inode))) {
if (mnode->GetPos().y() < 0.01) {
mnode->SetFixed(true);
mnode->SetFieldVal(10); // field: temperature [K]
}
}
}
// Remember to add the mesh to the system!
sys.Add(my_mesh);
// Visualization of the FEM mesh.
// This will automatically update a triangle mesh (a ChVisualShapeTriangleMesh
// 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 ChVisualShapeTriangleMesh).
// Do not forget AddVisualShapeFEA() at the end!
// Paint the colored mesh with temperature scale (NODE_FIELD_VALUE is the scalar field of the Poisson problem)
auto mvisualizemesh = chrono_types::make_shared<ChVisualShapeFEA>(my_mesh);
mvisualizemesh->SetFEMdataType(ChVisualShapeFEA::DataType::NODE_FIELD_VALUE);
mvisualizemesh->SetColorscaleMinMax(-1, 12);
mvisualizemesh->SetShrinkElements(false, 0.85);
mvisualizemesh->SetSmoothFaces(true);
my_mesh->AddVisualShapeFEA(mvisualizemesh);
// This will paint the wireframe
auto mvisualizemeshB = chrono_types::make_shared<ChVisualShapeFEA>(my_mesh);
mvisualizemeshB->SetFEMdataType(ChVisualShapeFEA::DataType::SURFACE);
mvisualizemeshB->SetWireframe(true);
my_mesh->AddVisualShapeFEA(mvisualizemeshB);
// This will paint the heat flux as line vectors
auto mvisualizemeshC = chrono_types::make_shared<ChVisualShapeFEA>(my_mesh);
mvisualizemeshC->SetFEMdataType(ChVisualShapeFEA::DataType::NONE);
mvisualizemeshC->SetFEMglyphType(ChVisualShapeFEA::GlyphType::ELEM_VECT_DP);
mvisualizemeshC->SetSymbolsScale(0.003);
mvisualizemeshC->SetDefaultSymbolsColor(ChColor(0.1f, 0.2f, 0.2f));
mvisualizemeshC->SetZbufferHide(false);
my_mesh->AddVisualShapeFEA(mvisualizemeshC);
// Create the Irrlicht visualization system
auto vis = chrono_types::make_shared<ChVisualSystemIrrlicht>();
vis->SetWindowSize(800, 600);
vis->SetWindowTitle("FEM thermal");
vis->Initialize();
vis->AddLogo();
vis->AddSkyBox();
vis->AddLight(ChVector3d(+20, 20, +20), 90, ChColor(0.5, 0.5, 0.5));
vis->AddLight(ChVector3d(-20, 20, -20), 90, ChColor(0.7f, 0.8f, 0.8f));
vis->AddCamera(ChVector3d(0, 0.7, -1), ChVector3d(0, 0.4, 0));
vis->AttachSystem(&sys);
// SIMULATION LOOP
// Use MINRES solver to handle stiffness matrices.
auto solver = chrono_types::make_shared<ChSolverMINRES>();
sys.SetSolver(solver);
solver->SetMaxIterations(150);
solver->SetTolerance(1e-6);
solver->EnableDiagonalPreconditioner(true);
solver->EnableWarmStart(true); // IMPORTANT for convergence when using EULER_IMPLICIT_LINEARIZED
solver->SetVerbose(false);
sys.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED); // fast, less precise
// Note: if you are interested only in a single LINEAR STATIC solution
// (not a transient thermal solution, but rather the steady-state solution),
// at this point you can uncomment the following line:
//
// sys.DoStaticLinear();
//
// Also, in the following while() loop, remove application.DoStep();
// so you can spin the 3D view and look at the solution.
while (vis->Run()) {
vis->BeginScene();
vis->Render();
vis->EndScene();
sys.DoStepDynamics(0.01);
if (sys.GetChTime() > 5)
break;
}
// Print some node temperatures..
for (unsigned int inode = 0; inode < my_mesh->GetNumNodes(); ++inode) {
if (auto mnode = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(inode))) {
if (mnode->GetPos().x() < 0.01) {
std::cout << "Node at y=" << mnode->GetPos().y() << " has T=" << mnode->GetFieldVal() << "\n";
}
}
}
return 0;
}
std::string GetChronoDataFile(const std::string &filename)
Get the full path to the specified filename, given relative to the Chrono data directory (thread safe...
Definition: ChGlobal.cpp:37
void Add(std::shared_ptr< ChPhysicsItem > item)
Attach an arbitrary ChPhysicsItem (e.g.
Definition: ChSystem.cpp:196
virtual void Render() override
Draw all 3D shapes and GUI elements at the current frame.
Definition: ChVisualSystemIrrlicht.cpp:570
virtual void Initialize() override
Initialize the visualization system.
Definition: ChVisualSystemIrrlicht.cpp:181
void AddSkyBox(const std::string &texture_dir=GetChronoDataFile("skybox/"))
Add a sky box in a 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:357
static void FromTetGenFile(std::shared_ptr< ChMesh > mesh, const char *filename_node, const char *filename_ele, std::shared_ptr< ChContinuumMaterial > my_material, ChVector3d pos_transform=VNULL, ChMatrix33<> rot_transform=ChMatrix33<>(1))
Load tetrahedrons from .node and .ele files as saved by TetGen.
Definition: ChMeshFileLoader.cpp:39
Namespace with classes for the Irrlicht module.
Definition: ChApiIrr.h:47
virtual void EndScene() override
End the scene draw at the end of each animation frame.
Definition: ChVisualSystemIrrlicht.cpp:560
Definition of a visual color.
Definition: ChColor.h:30
virtual void BeginScene() override
Perform any necessary operations at the beginning of each rendering frame.
Definition: ChVisualSystemIrrlicht.cpp:543
virtual int AddCamera(const ChVector3d &pos, ChVector3d targ=VNULL) override
Add a camera in an Irrlicht 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:290
int DoStepDynamics(double step_size)
Advance the dynamics simulation by a single time step of given length.
Definition: ChSystem.cpp:1632
Class for a physical system in which contact is modeled using a smooth (penalty-based) method.
Definition: ChSystemSMC.h:30
ChVector3< double > ChVector3d
Alias for double-precision vectors.
Definition: ChVector3.h:283
void SetWindowTitle(const std::string &win_title)
Set the windoiw title (default "").
Definition: ChVisualSystemIrrlicht.cpp:138
virtual void AttachSystem(ChSystem *sys) override
Attach another Chrono system to the run-time visualization system.
Definition: ChVisualSystemIrrlicht.cpp:165
virtual void SetSolver(std::shared_ptr< ChSolver > newsolver)
Attach a solver (derived from ChSolver) for use by this system.
Definition: ChSystem.cpp:319
virtual bool Run() override
Run the Irrlicht device.
Definition: ChVisualSystemIrrlicht.cpp:243
irr::scene::ILightSceneNode * AddLight(const ChVector3d &pos, double radius, ChColor color=ChColor(0.7f, 0.7f, 0.7f))
Add a point light to the scene.
Definition: ChVisualSystemIrrlicht.cpp:402
Main namespace for the Chrono package.
Definition: ChCamera.cpp:17
double GetChTime() const
Get the simulation time of this system.
Definition: ChSystem.h:154
void SetTimestepperType(ChTimestepper::Type type)
Set the method for time integration (time stepper type).
Definition: ChSystem.cpp:412
void AddLogo(const std::string &logo_filename=GetChronoDataFile("logo_chronoengine_alpha.png"))
Add a logo in a 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:336
Namespace for FEA classes.
Definition: ChVisualShapeFEA.h:28
void SetWindowSize(unsigned int width, unsigned int height)
Set the window size (default 640x480).
Definition: ChVisualSystemIrrlicht.cpp:134