新浪博客

一个EDEM中MBD耦合案例解析

2017-07-14 21:45阅读:
MBD(多体动力学)与EDEM耦合大致分2类:第一,利用插件EAlink实现耦合;第二,通过软件耦合模块coupling server 编程耦合。今天简单介绍下第二种方法,以一个案例来大体分析下程序代码的各部分构成。
目前市面上关于MBD编程实现耦合的资料非常稀缺,仅有一些案例文件,大致关于自由体运动、可以启闭的档板、刮板运动、升降斗运动、铲斗与自卸车作用、两圆辊破碎颗粒等。下面以free body为例简单分析下源程序 cpp的结构组成。源程序如下:


// Copyright DEM Solutions
// Mark Cook
// April 2009
//
// Modified Bill Edwards
// April 2013
// EDEM Coupling include files
#include 'IEDEMCouplingV2_0_0.h'
// Additional libraries
#include
#include
#include
using namespace std;
using namespace NApiEDEM;
class CGeometry // A simple geometry class
{
public:
// ID
int id;

// Geometry properties
double mass;
C3x3Matrix momentOfInertia;
C3x3Matrix initialMomentOfInertia;
// External forces
C3dVector force;
C3dVector torque;
// Geometry accelerations & velocities
C3dVector acceleration;
C3dVector angularAcceleration;
C3dVector velocity;
C3dVector angularVelocity;
// Geometry position
C3dPoint originalCenterOfMass;
C3dPoint centerOfMass;
C3x3Matrix orientation;
C3dVector totalTranslation;
} ;
// Define the ending simulation time
const double ENDTIME = 1.0;
// Define the number of time-steps EDEM runs for between data exchanges
const int TIME_STEP_RATIO = 10;
// This controls whether output is written to the console or not
const bool CONSOLEOUT = false;
// This controls whether output is written to a file or not
const bool FILEOUT = false;
int main()
{
// Simulation time-step
double dt;
// Set the simulation settings. This example must be started from time 0 secs
double simTime = 0.0;
// Create a box geometry
CGeometry box;
// Set the box properties
box.mass = 0.02;
box.initialMomentOfInertia = C3x3Matrix(1E-6, 0, 0,
0, 1E-6, 0,
0, 0, 1E-6);
///////////////////////////// Coupling Initialisation /////////////////////////////
///////////////////////////////////////////////////////////////////////////////////

IEDEMCoupling coupling; // Create an instance of the EDEM Coupling to use
// Initialise the coupling
if (!coupling.initialiseCoupling())
{
cout << 'Can't intialise the EDEM Coupling Client' << endl;
exit(EXIT_FAILURE);
}
cout << 'EDEM Coupling Client initialised' << endl << 'Connecting to EDEM...' << endl;
// Connect to EDEM
if(!coupling.connectCoupling())
{
cout << 'Could not connect to EDEM' << endl;
exit(EXIT_FAILURE);
}
cout << 'Connection to EDEM successful' << endl;
////////////////////////// Obtain Simulation Variables /////////////////////////////
////////////////////////////////////////////////////////////////////////////////////
// Set the EDEM time to zero

// If you have your options set to load last time-step in the EDEM
// Simulator it will cause problems with reloading the last saved
// time-step.
// Either disable the option or save your EDEM deck at time-step zero
// before simulating
coupling.setEDEMTime(simTime);
coupling.getEDEMTimeStep(dt);
dt *= TIME_STEP_RATIO;
// Get the geometry ID
if(coupling.getGeometryId('box', box.id))
{
cout << 'Found the geometry' << endl;
}
// Get the initial coupling position and orientation from EDEM
coupling.getGeometryPosition(box.id, box.originalCenterOfMass, box.orientation);

// Prepare data write out file if required
if(FILEOUT == true)
{
//open file
fstream DataOut( 'Data_Output.csv', ios::out|ios::app );
// Add date and time
time_t rawtime;
time ( &rawtime );
DataOut << ctime (&rawtime) << endl;
// Add headers
DataOut << 'Time,Force(X),Force(Y),Force(Z),'
<< 'Torque(X),Torque(Y),Torque(Z),'
<< 'Vel(X),Vel(Y),Vel(Z),'
<< 'Ang Vel(X),Ang Vel(Y),Ang Vel(Z),'
<< 'Angle,'
<< 'Axis(X),Axis(Y),Axis(Z),'
<< endl;
DataOut.close();
}
///////////////////////////////Main Simulation Loop ////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////
while (simTime < ENDTIME)
{
simTime += dt;
// Get geometry forces
coupling.getGeometryForces(box.id, box.force, box.torque);
// Calculate the geometry translation
// Acceleration
box.acceleration = box.force/ box.mass;
// Velocity
box.velocity += box.acceleration * dt;
// linear translation = position + velocity * dt
box.totalTranslation += box.velocity * dt;
box.centerOfMass = box.originalCenterOfMass + box.totalTranslation;
// Calculate the geometry rotation
// Angular acceleration
box.angularAcceleration = box.torque * box.momentOfInertia.inv();
// Angular velocity
box.angularVelocity += box.angularAcceleration * dt;
// Calculate the axis about which roation occurs
C3dVector rotationAxis = box.angularVelocity/ box.angularVelocity.length();
// Angle rotated about the velocity axis
double rotationAngle = box.angularVelocity.length() * dt;
// Change in orientation
C3x3Matrix orientationChange = C3x3Matrix(rotationAxis, rotationAngle);
// Calculate the new orienation of the geometry
box.orientation *= orientationChange;
// Update the moment of inertia
box.momentOfInertia = box.orientation.transpose() * box.initialMomentOfInertia * box.orientation;
//////////////////////// Update EDEM with Calculated Motion /////////////////////////
/////////////////////////////////////////////////////////////////////////////////////

// Send motion data to EDEM
coupling.setGeometryMotion( box.id,
box.totalTranslation,
box.orientation,
box.velocity,
box.angularVelocity,
dt); // Action time should be a multiple
// of the EDEM time-step
// Tell EDEM to perform the simulation step and pause the simulation when done
coupling.simulate(dt, ENDTIME);
//////////////////////////////////// Reporting //////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////
if (CONSOLEOUT == true)
{
cout << 'Time = ' << simTime << endl;
cout << ' Input Force = ' << box.force.x() << ' ' << box.force.y()
<< ' ' << box.force.z() << endl;
cout << ' Input Moment = ' << box.torque.x() << ' ' << box.torque.y()
<< ' ' << box.torque.z() << endl;
cout << ' Velocity = ' << box.velocity.x() << ' ' << box.velocity.y()
<< ' ' << box.velocity.z() << endl;
cout << ' Ang Velocity = ' << box.angularVelocity.x() << ' '
<< box.angularVelocity.y() << ' ' << box.angularVelocity.z() << endl;
cout << ' Rot Angle = ' << rotationAngle << endl;
cout << ' Rot Axis = ' << rotationAxis.x() << ' '
<< rotationAxis.y() << ' ' << rotationAxis.z() << endl;
cout << ' Rotation Matrix = ' << endl;
cout << ' ' << box.orientation.XX() << ' ' << box.orientation.XY() << ' ' << box.orientation.XZ() << endl;
cout << ' ' << box.orientation.YX() << ' ' << box.orientation.YY() << ' ' << box.orientation.YZ() << endl;
cout << ' ' << box.orientation.ZX() << ' ' << box.orientation.ZY() << ' ' << box.orientation.ZZ() << endl;
cout << ' MoI Matrix = ' << endl;
cout << ' ' << box.orientation.XX() << ' ' << box.orientation.XY() << ' ' << box.orientation.XZ() << endl;
cout << ' ' << box.orientation.YX() << ' ' << box.orientation.YY() << ' ' << box.orientation.YZ() << endl;
cout << ' ' << box.orientation.ZX() << ' ' << box.orientation.ZY() << ' ' << box.orientation.ZZ() << endl;
}
if (FILEOUT == true)
{
//open file
fstream DataOut( 'Data_Output.csv', ios::out|ios::app );
cout << simTime << ',';
cout << box.force.x() << ',' << box.force.y() << ',' << box.force.z() << ',';
cout << box.torque.x() << ',' << box.torque.y() << ',' << box.torque.z() << ',';
cout << box.velocity.x() << ',' << box.velocity.y() << ',' << box.velocity.z() << ',';
cout << box.angularVelocity.x() << ',' << box.angularVelocity.y() << ',' << box.angularVelocity.z() << ',';
cout << rotationAngle << ',';
cout << rotationAxis.x() << ',' << rotationAxis.y() << ',' << rotationAxis.z() << endl;
//close file
DataOut.close();
}
}
return 0; // Program exited normally
}

1.关于结构

1.1整个程序只有一个main()主函数。
程序自上而下,首先在开头定义了几个要调用的函数库,如#include等。
1.2接着是class CGeometry{},这句是声明几何体属性变量名。如:几何体id(name),质量,力,力矩,速度,加速度,转动惯量等。
1.3main()函数。
1.3.1coupling initialisation
该部分是用来探测EDEM coupling server模块是否开启,如果开启,但执行下一步。如果没有,并在提示窗口返回could not connect to EDEM
1.3.2 obtain simulation variables.
该部分内容用来定义仿真初始状态时系统参数设置及几何体参数设置情况。比如,可在这部分按程序语句定义重力、几何体质量、速度(角速度、力矩)等参数。之后coupling server读取几何体id,并将程序语句中的相关设置参数传输到EDEM中,同时,将几何体的相关初始信息(如:centerofmass, orientation等传输到程序中,以便随后进行相关计算。)
1.3.3 main simulation loop.
这是程序进行计算的主要循环部分。可以对一个几何体或者多个几何体进行位移速度等参数计算。程序通过sim Time+=dt实现每个时间步长的迭加,通过语句coupling.getGeometryForces()来获取颗粒对几何体的力与力矩作用。随后再通过相关语句(牛顿第二定律)实现加速度、速度、位移等的计算。
1.3.4 update EDEM with Calculated Motion.
这部分是将1.3.3所计算出来的数据更新到EDEM中,实现几何体的位置、速度等的变化。

2.关于编程
市面上已有MBD案例程序的内容大体都是以上提到的这几个部分构成,不同之处主要表现在几何体数量、运动方式、预读txt文件等的不同。在做自己项目相关的仿真时,完全可以参考已有案例中的语句进行编程。



注:个人学习经验,难免有不到之处,仅供交流学习!

我的更多文章

下载客户端阅读体验更佳

APP专享