一个EDEM中MBD耦合案例解析_沙发老师

http://blog.sina.com.cn/u/2144762345

首页博文目录关于我

个人资料

微博

加好友发纸条

写留言加关注

博客等级：
博客积分：

博客访问：
关注人气：
获赠金笔：0支
赠出金笔：0支
荣誉徽章：

正文字体大小：大中小

一个EDEM中MBD耦合案例解析

(2017-07-14 21:45:25)

标签：

杂谈

分类： CAE

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

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文件等的不同。在做自己项目相关的仿真时，完全可以参考已有案例中的语句进行编程。

注：个人学习经验，难免有不到之处，仅供交流学习！

阅读┊ 收藏 ┊ 喜欢 ▼ ┊打印┊举报/Report

前一篇：workflow一键开启蚂蚁森林

后一篇：EDEM最新版耦合接口编译过程中问题说明

新浪BLOG意见反馈留言板　欢迎批评指正