Adjusting mass flow-inlet boundary to achieve a target thrust

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!