How can we achieve target thrust at a mass-flow-inlet boundary?

For external aero applications involving the modeling of a full aircraft with engine and pylon, client typically knows the massflow rate incoming to the engine fan inlet and also the thrust produced by the engine.

The mfo UDF/Scheme (or the new mfo in Fluent 6.1) can be used at the pressure-outlet boundary representing the engine fan inlet to achieve the target massflow rate. Typically, at the engine exit boundary (incoming to the external fluid domain), the massflow inlet boundary type is used where user specify the massflow rate (known) and the total temperature. While the total temperature at the engine exit is usually known, client often prefers to specify the thrust instead.
How to do this ?
User can manually adjust the total temperature at the massflow-inlet boundary until the desired thrust is achieved. However, this manual procedure is tedious and time consuming. A UDF similar to the mfo procedure has been written to adjust the total temperature automatically to acquire the target thrust.

As with the mfo UDF, to avoid non-convergence, this UDF should only be applied when the actual thrust is about 20% from the target value.

/**************************************************************************/
/* UDF to achieve target Thrust at a mass-flow-inlet b.c. */
/* How to use this UDF ? */
/* - For 2D and 3D problem, set the global parameter */
/* Factor to 1.0 */
/* - For 2D-Axisymmetric, set the global parameter */
/* Factor to 2.0*Pi */
/* - Modify the values of Thrust_target, alpha, Ttot_max and */
/* Ttot_min inside the DEFINE_PROFILE udf */
/* Thrust_target : targer thrust rate in Newton */
/* alpha : under-relaxation (default=0.5) */
/* Ttot_min : lower bound of maximum allowed */
/* total temp. change in K */
/* Ttot_max : upper bound of maximum allowed */
/* total temp. change in K */
/* - Compile this UDF according to the architecture and */
/* Fluent version that you are using */
/* - Attach the compiled library and hook the UDF Thrust_fix */
/* to a mass-flow-inlet boundary */
/************************************************************************/

#include "udf.h"

#define Pi 3.14159265359
#define Factor 1
#define gamma 1.4
#define R 287.0
#define solver 0 /* 0 for segregated, 1 for coupled */

float dTCompute(float Thrust_target,float alpha,
float Ttot_min,float Ttot_max,int zoneID)
{
cell_t c0;
face_t f;
float A[ND_ND],area,MFR,Temp,Vmag,VmagSq,VsoundSq;
float Ttot,dTtot,Ttot_new;
float Thrust_present;
float sumarea,sumMFR,sumTtot,sumT,sumVmag;
Thread *tc0,*tf;
Domain *domain;

domain = Get_Domain(1);
tf = Lookup_Thread(domain,zoneID);

sumarea = 0.0;
sumMFR = 0.0;
sumT = 0.0;
sumTtot = 0.0;
sumVmag = 0.0;

begin_f_loop(f,tf)
{
c0 = F_C0(f,tf);
tc0 = F_C0_THREAD(f,tf);
F_AREA(A,f,tf);
area = Factor*NV_MAG(A);

/* check for SOLVER */
if ( solver==0 )
MFR = F_FLUX(f,tf);
else
MFR = C_R(c0,tc0)*( C_U(c0,tc0)*A[0] + C_V(c0,tc0)*A[1] + C_W(c0,tc0)*A[2] );

MFR = Factor*fabs(MFR);
Vmag = MFR/(area*C_R(c0,tc0));
VmagSq = Vmag*Vmag;
VsoundSq = gamma*R*C_T(c0,tc0);

sumarea += area;
sumMFR += MFR;
sumT += C_T(c0,tc0)*MFR;
sumTtot += (C_T(c0,tc0)*(1+(0.5*(gamma-1)*VmagSq/VsoundSq)))*MFR;
sumVmag += Vmag*MFR;
}
end_f_loop(f,tf)

area = sumarea;
MFR = sumMFR;
Temp = sumT/MFR;
Ttot = sumTtot/MFR;
Vmag = sumVmag/MFR;
Thrust_present = sumVmag;

dTtot = alpha*Temp*( Thrust_target - Thrust_present)/Thrust_present;

Ttot_new = Ttot + dTtot;

if ( Ttot_new >= Ttot_max ) Ttot_new = Ttot_max;
if ( Ttot_new <= Ttot_min ) Ttot_new = Ttot_min;

Message("Mass Flow =%f n",MFR);
Message("Thrust actual =%f n",Thrust_present);
Message("Thrust target =%f n",Thrust_target);
Message("Total-Temp Prevoius =%f n",Ttot);
Message("Total-Temp change =%f n",dTtot);
Message("Total-Temp Now =%f n",Ttot_new);

return Ttot_new;
}


DEFINE_PROFILE(Thrust_fix,tf,position)
{
face_t f;
int zoneID ;
float Thrust_target,alpha,Ttot_min,Ttot_max,Ttot_new;

Thrust_target = 900.0; /* Target thrust in Newton */
alpha = 0.5; /* Relaxation */
Ttot_min = 20.0; /* Minimum total temperature in K */
Ttot_max = 1000.0; /* Maximum total temperature in K */

zoneID = THREAD_ID(tf);

Ttot_new = dTCompute(Thrust_target,alpha,Ttot_min,Ttot_max,zoneID);

begin_f_loop(f,tf)
{
F_PROFILE(f,tf,position) = Ttot_new;
}
end_f_loop(f,tf)
}





Show Form
No comments yet. Be the first to add a comment!