aerodc81.h File Reference

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	outa_t
struct	vam_t
struct	c81_data

Typedefs
typedef struct outa_t	outa_t
typedef struct vam_t	vam_t
typedef struct c81_data	c81_data

Enumerations
enum	{   OUTA_ALPHA = 1, OUTA_GAMMA = 2, OUTA_MACH = 3, OUTA_CL = 4,   OUTA_CD = 5, OUTA_CM = 6, OUTA_CLALPHA = 7, OUTA_ALF1 = 8,   OUTA_ALF2 = 9, OUTA_DAN = 10, OUTA_DAM = 11, OUTA_DCN = 12,   OUTA_DCM = 13, OUTA_LAST = 20 }
enum	{   VAM_DENSITY = 0, VAM_SOUND_CELERITY, VAM_CHORD, VAM_FORCE_POSITION,   VAM_BC_POSITION, VAM_TWIST, VAM_LAST }

Functions
int	c81_aerod2 (doublereal *W, const vam_t *VAM, doublereal *TNG, outa_t *OUTA, c81_data *data)
int	c81_aerod2_u (doublereal *W, const vam_t *VAM, doublereal *TNG, outa_t *OUTA, c81_data *data, long unsteadyflag)

Variables
const outa_t	outa_Zero

Typedef Documentation

typedef struct c81_data c81_data

typedef struct outa_t outa_t

typedef struct vam_t vam_t

Enumeration Type Documentation

anonymous enum

Enumerator
OUTA_ALPHA
OUTA_GAMMA
OUTA_MACH
OUTA_CL
OUTA_CD
OUTA_CM
OUTA_CLALPHA
OUTA_ALF1
OUTA_ALF2
OUTA_DAN
OUTA_DAM
OUTA_DCN
OUTA_DCM
OUTA_LAST

Definition at line 39 of file aerodc81.h.

     {
        // 0 unused?
        OUTA_ALPHA      = 1,
        OUTA_GAMMA      = 2,
        OUTA_MACH       = 3,
        OUTA_CL         = 4,
        OUTA_CD         = 5,
        OUTA_CM         = 6,
        // 7 unused?
        OUTA_CLALPHA    = 7,
        OUTA_ALF1       = 8,
        OUTA_ALF2       = 9,
        OUTA_DAN        = 10,
        OUTA_DAM        = 11,
        OUTA_DCN        = 12,
        OUTA_DCM        = 13,
        // 14-19 unused?

        OUTA_LAST       = 20
};

anonymous enum

Enumerator
VAM_DENSITY
VAM_SOUND_CELERITY
VAM_CHORD
VAM_FORCE_POSITION
VAM_BC_POSITION
VAM_TWIST
VAM_LAST

Definition at line 96 of file aerodc81.h.

     {
        VAM_DENSITY = 0,
        VAM_SOUND_CELERITY,
        VAM_CHORD,
        VAM_FORCE_POSITION,
        VAM_BC_POSITION,
        VAM_TWIST,

        VAM_LAST
};

Function Documentation

int c81_aerod2	(	doublereal *	W,
		const vam_t *	VAM,
		doublereal *	TNG,
		outa_t *	OUTA,
		c81_data *	data
	)

Definition at line 100 of file aerodc81.c.

References c81_data::ad, c81_data::al, outa_t::alpha, c81_data::am, grad::atan2(), vam_t::bc_position, outa_t::cd, vam_t::chord, outa_t::cl, outa_t::clalpha, outa_t::cm, grad::cos(), vam_t::density, grad::fabs(), vam_t::force_position, outa_t::gamma, get_coef(), get_dcla(), M_PI, outa_t::mach, c81_data::md, c81_data::ml, c81_data::mm, c81_data::NAD, c81_data::NAL, c81_data::NAM, c81_data::NMD, c81_data::NML, c81_data::NMM, vam_t::sound_celerity, grad::sqrt(), and c81_data::stall.

{
        /* 
         * velocita' del punto in cui sono calcolate le condizioni al contorno
         */
        doublereal v[3];
        doublereal vp, vp2, vtot;
        doublereal rho = VAM->density;
        doublereal cs = VAM->sound_celerity;
        doublereal chord = VAM->chord;
        
        doublereal cl = 0., cl0 = 0., cd = 0., cd0 = 0., cm = 0.;
        doublereal alpha, gamma, cosgam, mach, q;
        doublereal dcla;
        
        doublereal ca = VAM->force_position;
        doublereal c34 = VAM->bc_position;
        
        const doublereal RAD2DEG = 180.*M_1_PI;
        const doublereal M_PI_3 = M_PI/3.;

        enum { V_X = 0, V_Y = 1, V_Z = 2, W_X = 3, W_Y = 4, W_Z = 5 };
        
        /* 
         * porta la velocita' al punto di calcolo delle boundary conditions 
         */
        v[V_X] = W[V_X];
        v[V_Y] = W[V_Y] + c34*W[W_Z];
        v[V_Z] = W[V_Z] - c34*W[W_Y];
        
        vp2 = v[V_X]*v[V_X] + v[V_Y]*v[V_Y];
        vp = sqrt(vp2);
        
        vtot = sqrt(vp2 + v[V_Z]*v[V_Z]);
        
        /*
         * non considera velocita' al di sotto di 1.e-6
         * FIXME: rendere parametrico?
         */
        if (vp/cs < 1.e-6) {
                TNG[V_X] = 0.;
                TNG[V_Y] = 0.;
                TNG[V_Z] = 0.;
                TNG[W_X] = 0.;
                TNG[W_Y] = 0.;
                TNG[W_Z] = 0.;
                
                OUTA->alpha = 0.;
                OUTA->gamma = 0.;
                OUTA->mach = 0.;
                OUTA->cl = 0.;
                OUTA->cd = 0.;
                OUTA->cm = 0.;
                OUTA->clalpha = 0.;
                
                return 0;
        }
        
        /*
         * FIXME: gestire quali angoli indicano cosa
         *
         * Idea di base: capire in base ad un  criterio di linearita'
         * a quale angolo di incidenza il profilo stalla, quindi
         * oltre quell'angolo prendere la correzione per flusso trasverso
         * in un certo modo, altrimenti in un altro ?!?
         */
        alpha = atan2(-v[V_Y], v[V_X]);
        OUTA->alpha = alpha*RAD2DEG;  
        gamma = atan2(-v[V_Z], fabs(v[V_X]));   /* come in COE0 (aerod2.f) */
        /* gamma = atan2(-v[V_Z], vp); */               /* secondo me (?!?) */
        OUTA->gamma = gamma*RAD2DEG;
        
        if (fabs(gamma) > M_PI_3) {
                /* tanto ne viene preso il coseno ... */
                gamma = M_PI_3;
        }
        
        cosgam = cos(gamma);
        mach = (vtot*sqrt(cosgam))/cs;
        OUTA->mach = mach;
        
        /*
         * Note: all angles in c81 files MUST be in degrees
         */
        get_coef(data->NML, data->ml, data->NAL, data->al,
                        OUTA->alpha, mach, &cl, &cl0);
        get_coef(data->NMD, data->md, data->NAD, data->ad,
                        OUTA->alpha, mach, &cd, &cd0);
        get_coef(data->NMM, data->mm, data->NAM, data->am,
                        OUTA->alpha, mach, &cm, NULL);

        dcla = get_dcla(data->NML, data->ml, data->stall, mach);
        
/*
 * da COE0 (aerod2.f):
 * 
        ASLRF = ASLOP0
        IF(DABS(ALFA).LT.1.D-6) GOTO 10
        ASLRF = CLIFT/(ALFA*COSGAM)
        IF(ASLRF.GT.ASLOP0) ASLRF = ASLOP0
     10 CLIFT = ASLRF*ALFA
 *
 */

        /*
         * in soldoni: se si e' oltre il tratto lineare, prende la
         * secante con l'angolo corretto per la freccia (fino a 60 
         * gradi) e poi ricalcola il cl con l'angolo vero; in questo
         * modo il cl viene piu' grande di circa 1/cos(gamma) ma solo
         * fuori del tratto lineare.
         */
        dcla *= RAD2DEG;
        if (fabs(alpha) > 1.e-6) {
                doublereal dclatmp = (cl - cl0)/(alpha*cosgam);
                if (dclatmp < dcla) {
                        dcla = dclatmp;
                }
        }
        cl = cl0 + dcla*alpha;
        
        OUTA->cl = cl;
        OUTA->cd = cd;
        OUTA->cm = cm;
        OUTA->clalpha = dcla;

        q = .5*rho*chord*vp2;

        TNG[V_X] = -q*(cl*v[V_Y] + cd*v[V_X])/vp;
        TNG[V_Y] = q*(cl*v[V_X] - cd*v[V_Y])/vp;
        TNG[V_Z] = -q*cd0*v[V_Z]/vp;
        TNG[W_X] = 0.;
        TNG[W_Y] = -ca*TNG[V_Z];
        TNG[W_Z] = q*chord*cm + ca*TNG[V_Y];
        
        /* 
         * Radial drag (TNG[V_Z]) consistent with Harris, JAHS 1970
         * and with CAMRAD strip theory section forces 
         */

        return 0;
}

Here is the call graph for this function:

int c81_aerod2_u	(	doublereal *	W,
		const vam_t *	VAM,
		doublereal *	TNG,
		outa_t *	OUTA,
		c81_data *	data,
		long	unsteadyflag
	)

Definition at line 243 of file aerodc81.c.

References c81_data::ad, c81_data::al, outa_t::alf1, outa_t::alf2, outa_t::alpha, c81_data::am, grad::atan2(), vam_t::bc_position, outa_t::cd, vam_t::chord, outa_t::cl, outa_t::clalpha, outa_t::cm, grad::cos(), outa_t::dam, outa_t::dan, outa_t::dcm, outa_t::dcn, vam_t::density, grad::exp(), grad::fabs(), vam_t::force_position, outa_t::gamma, get_coef(), get_dcla(), M_PI, outa_t::mach, c81_data::md, c81_data::ml, c81_data::mm, c81_data::mstall, c81_data::NAD, c81_data::NAL, c81_data::NAM, c81_data::NMD, c81_data::NML, c81_data::NMM, grad::pow(), vam_t::sound_celerity, grad::sqrt(), and c81_data::stall.

Referenced by TheodorsenAeroData::AssJac(), C81AeroData::GetForces(), C81MultipleAeroData::GetForces(), and C81InterpolatedAeroData::GetForces().

{
        /* 
         * velocita' del punto in cui sono calcolate le condizioni al contorno
         */
        doublereal v[3];
        doublereal vp, vp2, vtot;
        doublereal rho = VAM->density;
        doublereal cs = VAM->sound_celerity;
        doublereal chord = VAM->chord;
        
        doublereal cl = 0., cl0 = 0., cd = 0., cd0 = 0., cm = 0.;
        doublereal alpha, gamma, cosgam, mach, q;
        doublereal dcla;
        
        doublereal ca = VAM->force_position;
        doublereal c34 = VAM->bc_position;

        const doublereal RAD2DEG = 180.*M_1_PI;
        const doublereal M_PI_3 = M_PI/3.;

        enum { V_X = 0, V_Y = 1, V_Z = 2, W_X = 3, W_Y = 4, W_Z = 5 };
        
        /* 
         * porta la velocita' al punto di calcolo delle boundary conditions 
         */
        v[V_X] = W[V_X];
        v[V_Y] = W[V_Y] + c34*W[W_Z];
        v[V_Z] = W[V_Z] - c34*W[W_Y];
        
        vp2 = v[V_X]*v[V_X] + v[V_Y]*v[V_Y];
        vp = sqrt(vp2);
        
        vtot = sqrt(vp2 + v[V_Z]*v[V_Z]);
        
        /*
         * non considera velocita' al di sotto di 1.e-6
         * FIXME: rendere parametrico?
         */

        if (vp/cs < 1.e-6) {
                TNG[V_X] = 0.;
                TNG[V_Y] = 0.;
                TNG[V_Z] = 0.;
                TNG[W_X] = 0.;
                TNG[W_Y] = 0.;
                TNG[W_Z] = 0.;
                
                OUTA->alpha = 0.;
                OUTA->gamma = 0.;
                OUTA->mach = 0.;
                OUTA->cl = 0.;
                OUTA->cd = 0.;
                OUTA->cm = 0.;
                OUTA->clalpha = 0.;
                
                return 0;
        }
        
        /*
         * FIXME: gestire quali angoli indicano cosa
         *
         * Idea di base: capire in base ad un  criterio di linearita'
         * a quale angolo di incidenza il profilo stalla, quindi
         * oltre quell'angolo prendere la correzione per flusso trasverso
         * in un certo modo, altrimenti in un altro ?!?
         */
        alpha = atan2(-v[V_Y], v[V_X]);
        OUTA->alpha = alpha*RAD2DEG;  
        gamma = atan2(-v[V_Z], fabs(v[V_X]));   /* come in COE0 (aerod2.f) */
        /* gamma = atan2(-v[V_Z], vp); */               /* secondo me (?!?) */
        OUTA->gamma = gamma*RAD2DEG;
        
        if (fabs(gamma) > M_PI_3) {
                /* tanto ne viene preso il coseno ... */
                gamma = M_PI_3;
        }
        
        cosgam = cos(gamma);
        mach = (vtot*sqrt(cosgam))/cs;
        OUTA->mach = mach;

        /*
         * mach cannot be more than .99 (see aerod.f)
         */
        if (mach > .99) {
                mach = .99;
        }

        /*
         * Compute cl, cd, cm based on selected theory
         */
        switch (unsteadyflag) {
        case 0: 

                /*
                 * Note: all angles in c81 files MUST be in degrees
                 */
                get_coef(data->NML, data->ml, data->NAL, data->al, 
                                OUTA->alpha, mach, &cl, &cl0);
                get_coef(data->NMD, data->md, data->NAD, data->ad, 
                                OUTA->alpha, mach, &cd, &cd0);
                get_coef(data->NMM, data->mm, data->NAM, data->am, 
                                OUTA->alpha, mach, &cm, NULL);

                dcla = get_dcla(data->NML, data->ml, data->stall, mach);
        
/*
 * da COE0 (aerod2.f):
 * 
        ASLRF = ASLOP0
        IF(DABS(ALFA).LT.1.D-6) GOTO 10
        ASLRF = CLIFT/(ALFA*COSGAM)
        IF(ASLRF.GT.ASLOP0) ASLRF = ASLOP0
     10 CLIFT = ASLRF*ALFA
 *
 */

                /*
                 * in soldoni: se si e' oltre il tratto lineare, prende la
                 * secante con l'angolo corretto per la freccia (fino a 60 
                 * gradi) e poi ricalcola il cl con l'angolo vero; in questo
                 * modo il cl viene piu' grande di circa 1/cos(gamma) ma solo
                 * fuori del tratto lineare.
                 */
                dcla *= RAD2DEG;
                if (fabs(alpha) > 1.e-6) {
                        doublereal dclatmp = (cl - cl0)/(alpha*cosgam);
                        if (dclatmp < dcla) {
                                dcla = dclatmp;
                                cl = cl0 + dcla*alpha;
                        }
                }
                break;

        case 1:
                return -1;

        case 2: {

                /*
                 * Constants from unsteady theory
                 * synthetized by Richard L. Bielawa,
                 * 31th A.H.S. Forum Washington D.C. 
                 * May 1975
                 */
                doublereal A, B, A2, B2, ETA, ASN, ASM, 
                        SGN, SGM, SGMAX, 
                        DAN, DCN, DAM, DCM, 
                        S2, alphaN, alphaM, C1,
                        dcma, dclatan, ALF1, ALF2,
                        cn;
                
                const doublereal PN[] = { 
                        -3.464003e-1, 
                        -1.549076e+0, 
                        4.306330e+1, 
                        -5.397529e+1,
                        5.781402e+0,
                        -3.233003e+1,
                        -2.162257e+1,
                        1.866347e+1,
                        4.198390e+1,
                        3.295461e+2,
                };
        
                const doublereal QN[] = {
                        1.533717e+0,
                        6.977203e+0,
                        1.749010e+3,
                        1.694829e+3,
                        -1.771899e+3,
                        -3.291665e+4,
                        2.969051e+0,
                        -3.632448e+1,
                        -2.268578e+3,
                        6.601995e+3,
                        -9.654208e+3, 
                        8.533930e+4,
                        -1.492624e+0,
                        1.163661e+1
                };

                const doublereal PM[] = {
                        1.970065e+1,
                        -6.751639e+1,
                        7.265269e+2,
                        4.865945e+4,
                        2.086279e+4,
                        6.024672e+3,
                        1.446334e+2,
                        8.586896e+2,
                        -7.550329e+2,
                        -1.021613e+1,
                        2.247664e+1,
                };
        
                const doublereal QM[] = {
                        -2.322808e+0,
                        -1.322257e+0,
                        -2.633891e+0,
                        -2.180321e-1,
                        4.580014e+0,
                        3.125497e-1,
                        -2.828806e+1,
                        -4.396734e+0,
                        2.565870e+2,
                        -1.204976e+1,
                        -1.157802e+2,
                        8.612138e+0,
                };

                enum {
                        U_1 = 0,
                        U_2 = 1,
                        U_3 = 2,
                        U_4 = 3,
                        U_5 = 4,
                        U_6 = 5,
                        U_7 = 6,
                        U_8 = 7,
                        U_9 = 8,
                        U10 = 9,
                        U11 = 10,
                        U12 = 11,
                        U13 = 12,
                        U14 = 13
                };

                /*
                 * This is the static stall angle for Mach = 0
                 * (here a symmetric airfoil is assumed; the real
                 * _SIGNED_ static stall should be considered ...)
                 */
                const doublereal ASN0 = .22689, ASM0 = .22689;

                ALF1 = OUTA->alf1;
                ALF2 = OUTA->alf2;

                A = .5*chord*ALF1/vp;
                B = .25*chord*chord*ALF2/vp2;

                ETA = sqrt(pow(A/.048, 2) + pow(B/.016, 2));

                if (alpha < 0.) {
                        A = -A;
                        B = -B;
                }

                if (ETA > 1.) {
                        A /= ETA;
                        B /= ETA;
                }

                A2 = A*A;
                B2 = B*B;

                ASN = ASN0*(1. - mach);
                ASM = ASM0*(1. - mach);

                SGN = fabs(alpha/ASN);
                SGM = fabs(alpha/ASM);
                SGMAX = 1.839 - 70.33*fabs(B);
                if (SGMAX > 1.86) {
                        SGMAX = 1.86;
                }
                if (SGN > SGMAX) {
                        SGN = SGMAX;
                }
                if (SGM > SGMAX) {
                        SGM = SGMAX;
                }

                DAN = (A*(PN[U_1] + PN[U_5]*SGN) + B*(PN[U_2] + PN[U_6]*SGN)
                        + exp(-1072.52*A2)*(A*(PN[U_3] + PN[U_7]*SGN)
                                +A2*(PN[U_9] + PN[U10]*SGN))
                        + exp(-40316.42*B2)*B*(PN[U_4] + PN[U_8]*SGN))*ASN;
                DCN = A*(QN[U_1] + QN[U_3]*A2 + SGN*(QN[U_7] + QN[U_9]*A2 + QN[U13]*SGN)
                                + B2*(QN[U_5] + QN[U11]*SGN))
                        + B*(QN[U_2] + QN[U_4]*A2
                                        + SGN*(QN[U_8] + QN[U10]*A2 + QN[U14]*SGN)
                                        + B2*(QN[U_6] + QN[U12]*SGN));
                DAM = (A*(PM[U_1] + PM[U_3]*A2 + PM[U_5]*B2 + PM[U10]*SGM + PM[U_7]*A)
                                + B*(PM[U_2] + PM[U_4]*B2 + PM[U_6]*A2
                                        + PM[U11]*SGM + PM[U_8]*B + PM[U_9]*A))*ASM;

                S2 = SGM*SGM;

                DCM = A*(QM[U_2] + QM[U_8]*A + SGM*(QM[U_4] + QM[U10]*A)
                                + S2*(QM[U_6] + QM[U12]*A))
                        + B*(QM[U_1] + QM[U_7]*B + SGM*(QM[U_3] + QM[U_9]*B)
                                        + S2*(QM[U_5] + QM[U11]*B));

                OUTA->dan = DAN*RAD2DEG;
                OUTA->dam = DAM*RAD2DEG;
                OUTA->dcn = DCN;
                OUTA->dcm = DCM;

                /* 
                 * I think I need to apply this contribution 
                 * with the sign of alpha, because otherwise 
                 * it gets discontinuous as alpha changes sign;
                 * I definitely need the original reference :(
                 */
                if (alpha < 0.) {
                        DAN = -DAN;
                        DCN = -DCN;
                        DAM = -DAM;
                        DCM = -DCM;
                }

                alphaN = (alpha - DAN)*RAD2DEG;
                get_coef(data->NML, data->ml, data->NAL, data->al, 
                                alphaN, mach, &cl, &cl0);
                get_coef(data->NMD, data->md, data->NAD, data->ad, 
                                alphaN, mach, &cd, &cd0);

                alphaM = (alpha - DAM)*RAD2DEG;
                get_coef(data->NMM, data->mm, data->NAM, data->am, 
                                alphaM, mach, &cm, NULL);

                dcla = get_dcla(data->NML, data->ml, data->stall, mach);
                dcma = get_dcla(data->NMM, data->mm, data->mstall, mach);

                /* note: cl/alpha in 1/deg */
                dclatan = dcla;
                if (fabs(alphaN) > 1.e-6) {
                        dclatan = (cl - cl0)/(alphaN*cosgam);
                }
                cl = cl0 + dclatan*alphaN;

                /* back to 1/rad */
                dcla *= RAD2DEG;
                dcma *= RAD2DEG;

                C1 = .9457/sqrt(1. - mach*mach);

                /* 
                 * the unsteady correction is "cn", 
                 * so split it in "cl" and "cd"
                 * (note: if "vp" is too small the routine exits
                 * earlier without computing forces)
                 */
                cn = dcla*DAN + DCN*C1;
                cl += cn*v[V_X]/vp;     /* cos(x) */
                cd -= cn*v[V_Y]/vp;     /* sin(x) */
                cm += dcma*DAM + DCM*C1;

                break;
        }
        }

        /*
         * Save cl, cd, cm for output purposes
         */
        OUTA->cl = cl;
        OUTA->cd = cd;
        OUTA->cm = cm;
        OUTA->clalpha = dcla;

        /*
         * Local dynamic pressure
         */
        q = .5*rho*chord*vp2;

        /*
         * airfoil forces and moments in the airfoil frame
         */
        TNG[V_X] = -q*(cl*v[V_Y] + cd*v[V_X])/vp;
        TNG[V_Y] = q*(cl*v[V_X] - cd*v[V_Y])/vp;
        TNG[V_Z] = -q*cd0*v[V_Z]/vp;
        TNG[W_X] = 0.;
        TNG[W_Y] = -ca*TNG[V_Z];
        TNG[W_Z] = q*chord*cm + ca*TNG[V_Y];
        
        /* 
         * Radial drag (TNG[V_Z]) consistent with Harris, JAHS 1970
         * and with CAMRAD strip theory section forces 
         */

        return 0;
}

Here is the call graph for this function:

Variable Documentation

const outa_t outa_Zero

Definition at line 87 of file aerodc81.c.