[转载]fluent学习----关于six dof的神作求解释

2017-01-14 by:CAE仿真在线来源:互联网

原文地址:fluent学习----关于six dof的神作求解释作者:松松

本人在中文流体网上见了一篇6dof的资料,但对里面的last_call代表什么意思不能理解,还有在编译时出错是怎么回事?好像不认识boolean这个运算符 ,我算的是单刚体,是否可以把last_call去掉
....src6dof.c(195) : error C2146: syntax error : missing ';' before identifier 'last_call'
....src6dof.c(195) : error C2065: 'last_call' : undeclared identifier
望大虾不吝赐教
程序如下:
#include "udf.h"
#define ZONE_ID1 8

#define MASS 907.185
#define IXX 27.116
#define IYY 488.094
#define IZZ 488.094
#define IXZ 0.0
#define IXY 0.0
#define IYZ 0.0
#define NUM_CALLS 2
#define R2D 180./M_PI

#if !RP_NODE
static void
write_velocities (real *v_old, real *om_old, real *euler)
{
FILE *file_velo;
register int i;
if ((file_velo = fopen ("velocities_3d", "w")) == NULL)
Message ("nCannot open velocities_3d file");
for (i = 0; i < ND_ND; i++)
fprintf (file_velo, "%e ", v_old);
for (i = 0; i < 3; i++)
fprintf (file_velo, "%e ", om_old);
for (i = 0; i < 3; i++)
fprintf (file_velo, "%e ", euler);
fclose(file_velo);
}
#endif
static void
read_velocities (real *v_old, real *om_old, real *euler)
{
float read_v_old[3], read_om_old[3], read_euler[3];
register int i;
#if !RP_NODE
FILE *file_velo;
if ((file_velo = fopen ("velocities_3d", "r")) == NULL)
{
Message ("nCannot open velocities_3d file...creating one!");
if ((file_velo = fopen ("velocities_3d", "w")) == NULL)
Message ("nCannot create velocities_3d file...dying!");
else
{
for (i = 0; i < 9; i++)
fprintf (file_velo, "%e ", 0.0);
fclose (file_velo);
file_velo = fopen ("velocities_3d", "r");
}
}
for (i = 0; i < ND_ND; i++)
fscanf(file_velo, "%e", &read_v_old);
# if RP_2D
read_v_old[2] = 0.0;
# endif
for (i = 0; i < 3; i++)
fscanf(file_velo, "%e", &read_om_old);
for (i = 0; i < 3; i++)
fscanf(file_velo, "%e", &read_euler);
fclose(file_velo);
#endif
#if PARALLEL
{
real exchange[9];
# if RP_HOST
for (i = 0; i < 3; i++)
{
exchange= (real)read_v_old;
exchange[i + 3] = (real)read_om_old;
exchange[i + 6] = (real)read_euler;
}
# endif
host_to_node_real (exchange, 9);
# if RP_NODE
for (i = 0; i < 3; i++)
{
read_v_old= (float)exchange;
read_om_old= (float)exchange[i + 3];
read_euler= (float)exchange[i + 6];
}
# endif
}
#endif
for (i = 0; i < 3; i++)
{
v_old= (real)read_v_old;
om_old= (real)read_om_old;
euler= (real)read_euler;
}
}

DEFINE_CG_MOTION(six_dof, dt, cg_vel, cg_omega, time, dtime)

{
real f_glob[3];
real af_glob[3];
real af_body[3];
real m_glob[3];
real m_body[3];
real v_old_glob[3];
real om_old_body[3];
real om_body[3];
real al_body[3];
real x_cg[3];
real euler_angle[3];
real euler_angle_old[3];
real euler_rate[3];
real grav_acc[3];
real cf, ct, cs;
real sf, st, ss;
real RC[3][3];
real RD[3][3];

real wx, wy, wz;
real r1, r2, r3;
real s1, s2, s3;
real D;
real a11, a21, a31, a12, a22, a32, a13, a23, a33;
register int i;
Domain *domain = Get_Domain (1);
Thread *tf1 = Lookup_Thread (domain, ZONE_ID1);
static int calls = 0;
static boolean last_call = FALSE;

grav_acc[0] = 0.0;
grav_acc[1] = 0.0;
grav_acc[2] = 9.81;

read_velocities (v_old_glob, om_old_body, euler_angle);

for (i = 0; i < ND_ND; i++)
x_cg= DT_CG(dt);

if ((++calls) == NUM_CALLS)
last_call = TRUE;

Compute_Force_And_Moment (domain, tf1, x_cg, f_glob, m_glob, TRUE);
for (i = 0; i < ND_ND; i++)
af_glob= f_glob;

for (i = 0; i < ND_ND; i++)
f_glob+= grav_acc*MASS;
add_injector_forces (x_cg, euler_angle, f_glob, m_glob);

for (i = 0; i < ND_ND; i++)
cg_vel= f_glob* dtime / MASS + v_old_glob;

for (i = 0; i < ND_ND; i++)
x_cg= x_cg+ cg_vel* dtime;

cf = cos(euler_angle[0]);
ct = cos(euler_angle[1]);
cs = cos(euler_angle[2]);
sf = sin(euler_angle[0]);
st = sin(euler_angle[1]);
ss = sin(euler_angle[2]);
RC[0][0] = ct*cs;
RC[0][1] = sf*st*cs-cf*ss;
RC[0][2] = cf*st*cs+sf*ss;
RC[1][0] = ct*ss;
RC[1][1] = sf*st*ss+cf*cs;
RC[1][2] = cf*st*ss-sf*cs;
RC[2][0] = -st;
RC[2][1] = sf*ct;
RC[2][2] = cf*ct;
RD[0][0] = 1.0;
RD[0][1] = sf*st/ct;
RD[0][2] = cf*st/ct;
RD[1][0] = 0.0;
RD[1][1] = cf;
RD[1][2] = -sf;
RD[2][0] = 0.0;
RD[2][1] = sf/ct;
RD[2][2] = cf/ct;

NV_ROTN_V (m_body, = , m_glob, RC);
NV_ROTN_V (af_body, = , af_glob, RC);

wx = om_old_body[0];
wy = om_old_body[1];
wz = om_old_body[2];
r1 = IXX * wx - IXY * wy - IXZ * wz;
r2 = -IXY * wx + IYY * wy - IYZ * wz;
r3 = -IXZ * wx - IYZ * wy + IZZ * wz;
s1 = m_body[0] - (wy*r3-wz*r2);
s2 = m_body[1] - (wz*r1-wx*r3);
s3 = m_body[2] - (wx*r2-wy*r1);
D = IXX * IYY * IZZ - IXX * IYZ * IYZ -
IYY * IXZ * IXZ - IZZ * IXY * IXY -
2.* IXY * IXZ * IYZ;
a11 = IYY * IZZ - IYZ * IYZ;
a12 = IXY * IZZ + IYZ * IXZ;
a13 = IXY * IYZ + IXZ * IYY;
a21 = IXY * IZZ + IYZ * IXZ;
a22 = IXX * IZZ - IXZ * IXZ;
a23 = IXX * IYZ + IXZ * IXY;
a31 = IXY * IYZ + IXZ * IYY;
a32 = IXX * IYZ + IXZ * IXY;
a33 = IXX * IYY - IXY * IXY;

al_body[0] = (a11 * s1 + a12 * s2 + a13 * s3) / D;
al_body[1] = (a21 * s1 + a22 * s2 + a23 * s3) / D;
al_body[2] = (a31 * s1 + a32 * s2 + a33 * s3) / D;

for (i = 0; i < 3; i++)
om_body= om_old_body+ al_body* dtime;

NV_ROTP_V (cg_omega, = ,om_body, RC);

for (i = 0; i < 3; i++)
euler_angle_old= euler_angle;
NV_ROTP_V (euler_rate, =, om_body, RD);
N3V_S (euler_rate, *=, dtime);
N3V_V (euler_angle, +=, euler_rate);

#if !RP_NODE

if (last_call)
{
FILE *file_diag;

if ((file_diag = fopen ("diagnostics_3d", "a")) == NULL)
Message ("nCannot open output file");
fprintf (file_diag, "%ft" , time);
fprintf (file_diag, "%ft%ft%ft", x_cg[0], x_cg[1], x_cg[2]);
fprintf (file_diag, "%ft%ft%ft", cg_vel[0], cg_vel[1], cg_vel[2]);
fprintf (file_diag, "%ft%ft%ft", om_body[0], om_body[1], om_body[2]);
fprintf (file_diag, "%ft%ft%ft",
euler_angle[0], euler_angle[1], euler_angle[2]);
fprintf (file_diag, "%ft%ft%ft", af_body[0], af_body[1], af_body[2]);
fprintf (file_diag, "%ft%ft%fn", m_body[0], m_body[1], m_body[2]);
fclose(file_diag);

write_velocities (cg_vel, om_body, euler_angle);

calls = 0;
last_call = FALSE;
}
#endif
}

来自http://www.baisi.net/thread-2278951-1-1.html

http://www.efluid.com.cn/bbs/showtopic-4153.aspx

学习six-dof其他资料

http://http://forums.caenet.cn/showtopic-531777.aspx

http://forums.caenet.cn/showtopic-532886.aspx

开放分享：优质有限元技术文章,助你自学成才