**Appendix A**

#include "udf.h" #include "stdio.h" #include "time.h" #include "stdlib.h" #include "math.h" static real speed[ND\_ND]; static real weiyi\_before[ND\_ND],sudu\_before[ND\_ND],jiasudu\_before[ND\_ND]; static real speed\_nz; static real weiyi\_before\_nz,sudu\_before\_nz,jiasudu\_before\_nz; static real number\_1 = 0; static real number\_2 = 0; DEFINE\_CG\_MOTION(first,dt,vel,omega,time,dtime) { real m[ND\_ND],k[ND\_ND],zunibi[ND\_ND],zunixishu[ND\_ND],freq\_rad[ND\_ND]; real beta,gama; real k\_equ[ND\_ND],inertia\_equ[ND\_ND],damping\_equ[ND\_ND],force \_equ[ND\_ND];

```
real code_1,code_2,code_3,code_4,code_5,code_6,code_7;
real weiyi[ND_ND],sudu[ND_ND],jiasudu[ND_ND];
real m_nz,k_nz,zunibi_nz,zunixishu_nz,freq_rad_nz;
real beta_nz,gama_nz;
real k_equ_nz,inertia_equ_nz,damping_equ_nz,force_equ_nz;
real code_1nz,code_2nz,code_3nz,code_4nz,code_5nz,code_6nz,code_7nz;
real weiyi_nz,sudu_nz,jiasudu_nz;
real hezai,force,moment,li[ND_ND],liju[ND_ND];
real cog[ND_ND];
int x, y, z;
#if !RP_NODE
FILE *fp_1;
#endif
#if !RP_HOST
Domain *d;
Thread *t;
face_t f;
cell_t c;
real NV_VEC (A);
d = Get_Domain(1);
t = DT_THREAD(dt);
for(x=0; x<ND_ND; x++)
      {
           cog[x] = DT_CG(dt)[x];
      }
Compute_Force_And_Moment(d, t, cog, li, liju, FALSE);
hezai = li[0];
force = li[1];
moment = liju[2];
#endif
node_to_host_real_4(hezai,force,moment,cog[1]);
if (number_1 == 0)
      {
           weiyi_before[1] = 0;
           sudu_before[1] = 0;
           jiasudu_before[1] = 0;
           weiyi_before[0] = 0;
           sudu_before[0] = 0;
           jiasudu_before[0] = 0;
           weiyi_before_nz = 0;
           sudu_before_nz = 0;
           jiasudu_before_nz = 0;
           number_1 = number_1 + 1;
      }
m[1] = 1; /**the mass of vertical movement, Kg**/
k[1] =1; /**the stiffness of vertical movement, N/m**/
zunixishu[1] =1;
gama = 0.5;
beta = 0.25;
code_1 = gama/(beta*dtime);
code_2 = 1/(beta*dtime*dtime);
code_3 = (1-gama)*dtime;
code_4 = gama/beta;
code_5 = gama*dtime*(0.5-beta)/beta;
```
code\_6 = 1/(beta\*dtime); code\_7 = (0.5-beta)/beta; k\_equ[1] = k[1] + code\_1\*zunixishu[1] + code\_2\*m[1]; inertia\_equ[1] = m[1] \* (code\_2\*weiyi\_before[1] + code\_6\*sudu\_before[1] + code\_7\*jiasudu\_before[1]); damping\_equ[1] = zunixishu[1] \* (-1\*sudu\_before[1] - code\_3\*jiasudu\_before[1] + code\_1\*weiyi\_before[1] + code\_4\*sudu\_before[1] + code\_5\*jiasudu\_before[1]); force\_equ[1] = force + inertia\_equ[1] + damping\_equ[1]; weiyi[1] = force\_equ[1]/k\_equ[1]; sudu[1] = sudu\_before[1] + code\_3\*jiasudu\_before[1] + code\_1\*(weiyi[1]-weiyi\_ before[1]) - code\_4\*sudu\_before[1] - code\_5\*jiasudu\_before[1]; jiasudu[1] = (force - zunixishu[1]\*sudu[1] - k[1]\*weiyi[1])/m[1]; weiyi\_before[1] = weiyi[1]; sudu\_before[1] = sudu[1]; jiasudu\_before[1] = jiasudu[1]; speed[1] = sudu[1]; #if !RP\_NODE fp\_1 = fopen("Fluent\_couple.dat","a+"); fprintf(fp\_1, "%.16lf %.16lf %.16lf %.16lf %.16lf %.16lf %.16lf \n", time, hezai, force, moment, weiyi[1], sudu[1], jiasudu[1]); /\*\*1\_shijian, 2\_zuli, 3\_shengli, 4\_liju, 5\_niuzhuanjiao(hudu), 6\_jiaosudu, 7\_jiaojiasudu, 8\_shuxiangweiyi(m), 9\_shuxiangsudu, 10\_shuxiangjiasudu, 11\_zhongxinshuxiangweizhi\*\*/ fclose(fp\_1); #endif host\_to\_node\_real\_3(sudu[0],sudu[1],sudu\_nz); #if !RP\_HOST NV\_S(vel, =, 0.0); NV\_S(omega, =, 0.0); if (!Data\_Valid\_P()) return; vel[1] = sudu[1]; #endif } DEFINE\_CG\_MOTION(second,dt,vel,omega,time,dtime) { host\_to\_node\_real\_3(speed[0],speed[1], speed\_nz); #if !RP\_HOST vel[1] = speed[1]; #endif }
