ERF_AdvanceGeneralAD.cpp File Reference

#include <ERF_GeneralAD.H>
#include <ERF_IndexDefines.H>
#include <ERF_Interpolation_1D.H>
#include <cmath>
#include <ERF_NumericalConstants.H>

Include dependency graph for ERF_AdvanceGeneralAD.cpp:

Functions
AMREX_FORCE_INLINE AMREX_GPU_DEVICE int	find_rad_loc_index (const Real rad, const Real *bld_rad_loc, const int n_bld_sections)
AMREX_FORCE_INLINE AMREX_GPU_DEVICE std::array< Real, 2 >	compute_source_terms_Fn_Ft (const Real rad, const Real avg_vel, const Real *bld_rad_loc, const Real *bld_twist, const Real *bld_chord, int n_bld_sections, const Real *bld_airfoil_aoa, const Real *bld_airfoil_Cl, const Real *bld_airfoil_Cd, const int n_pts_airfoil, const Real *velocity, const Real *rotor_RPM, const Real *blade_pitch, const int n_spec_extra)

Function Documentation

compute_source_terms_Fn_Ft()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE std::array<Real,2> compute_source_terms_Fn_Ft	(	const Real	rad,
		const Real	avg_vel,
		const Real *	bld_rad_loc,
		const Real *	bld_twist,
		const Real *	bld_chord,
		int	n_bld_sections,
		const Real *	bld_airfoil_aoa,
		const Real *	bld_airfoil_Cl,
		const Real *	bld_airfoil_Cd,
		const int	n_pts_airfoil,
		const Real *	velocity,
		const Real *	rotor_RPM,
		const Real *	blade_pitch,
		const int	n_spec_extra
	)

{
 
    Real rpm = interpolate_1d(velocity, rotor_RPM, avg_vel, n_spec_extra);
    Real pitch = interpolate_1d(velocity, blade_pitch, avg_vel, n_spec_extra);
 
    Real Omega = rpm/Real(60.0)*two*PI;
    Real rho = Real(1.226);
 
    Real B = three;
    Real rhub = two;
    Real rtip = Real(63.5);
 
    Real twist = interpolate_1d(bld_rad_loc, bld_twist, rad, n_bld_sections);
    Real c = interpolate_1d(bld_rad_loc, bld_chord, rad, n_bld_sections);
 
    // Iteration procedure
 
    Real s = myhalf*c*B/(PI*rad);
 
    Real at, an, V1, Vt, Vr, psi, L, D, Cn, Ct;
    Real ftip, fhub, F, Cl, Cd, at_new, an_new;
 
    at = Real(0.1);
    an = Real(0.1);
 
    bool is_converged = false;
 
    for(int i=0;i<100;i++) {
        V1 = avg_vel*(1-an);
        Vt = Omega*(one+at)*rad;
        Vr = std::pow(V1*V1+Vt*Vt,myhalf);
 
        psi = std::atan2(V1,Vt);
 
        Real aoa = psi*Real(180.0)/PI - twist + pitch;
 
        Cl = interpolate_1d(bld_airfoil_aoa, bld_airfoil_Cl, aoa, n_pts_airfoil);
        Cd = interpolate_1d(bld_airfoil_aoa, bld_airfoil_Cd, aoa, n_pts_airfoil);
 
        //Cl = Real(1.37);
        //Cd = Real(0.014);
 
        //printf("rad, aoa, Cl, Cd    = %0.15g %0.15g %0.15g %0.15g\n", rad, aoa, Cl, Cd);
 
        Cn = Cl*std::cos(psi) + Cd*std::sin(psi);
        Ct = Cl*std::sin(psi) - Cd*std::cos(psi);
 
        ftip = B*(rtip-rad)/(two*rad*std::sin(psi)+Real(1e-10));
        fhub = B*(rad-rhub)/(two*rad*std::sin(psi)+Real(1e-10));
 
        AMREX_ALWAYS_ASSERT(std::fabs(std::exp(-fhub))<=one);
        AMREX_ALWAYS_ASSERT(std::fabs(std::exp(-ftip))<=one);
 
        F = two/PI*(std::acos(std::exp(-ftip)) + std::acos(std::exp(-fhub)) );
 
        at_new = one/ ( Real(4.0)*F*std::sin(psi)*std::cos(psi)/(s*Ct+Real(1e-10)) - one );
        an_new = one/ ( one + Real(4.0)*F*amrex::Math::powi<2>(std::sin(psi))/(s*Cn + Real(1e-10)) );
        at_new = std::max(Real(0.), at_new);
 
        if(std::fabs(at_new-at) < Real(1e-5) and std::fabs(an_new-an) < Real(1e-5)) {
            //printf("Converged at, an = %d %0.15g %0.15g %0.15g\n",i, at, an, psi);
            at = at_new;
            an = an_new;
            is_converged = true;
            break;
        }
        at = at_new;
        an = an_new;
        //printf("Iteration, at, an = %0.15g %0.15g %0.15g\n",at, an, psi);
    }
 
    if(!is_converged) {
        Abort("The iteration procedure for the generalized actuator disk did not converge. Exiting...");
    }
 
    // Iterations converged. Now compute Ft, Fn
 
    L = myhalf*rho*Vr*Vr*c*Cl;
    D = myhalf*rho*Vr*Vr*c*Cd;
 
    Real Fn = L*std::cos(psi) + D*std::sin(psi);
    Real Ft = L*std::sin(psi) - D*std::cos(psi);
 
     //printf("Fn and Ft %0.15g %0.15g %0.15g %0.15g\n", L, D, std::cos(psi), std::sin(psi));
 
    std::array<Real, 2> Fn_and_Ft;
    Fn_and_Ft[0] = Fn;
    Fn_and_Ft[1] = Ft;
 
    return Fn_and_Ft;
 
    //exit(0);
}

Referenced by GeneralAD::source_terms_cellcentered().

Here is the call graph for this function:

Here is the caller graph for this function:

find_rad_loc_index()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE int find_rad_loc_index	(	const Real	rad,
		const Real *	bld_rad_loc,
		const int	n_bld_sections
	)

{
      // Find the index of the radial location
    int index=-1;
    Real rhub = two;
    Real rad_from_hub = rad - rhub;
    if(rad_from_hub < zero) {
        index = 0;
    }
    else {
        for(int i=0;i<n_bld_sections;i++){
            if(bld_rad_loc[i] > rad) {
                index = i;
                break;
            }
        }
    }
    if(index == -1 and rad > bld_rad_loc[n_bld_sections-1]) {
        index = n_bld_sections-1;
    }
    if(index == -1) {
        //printf("The radial section is at %0.15g m\n",rad);
        Abort("Could not find index of the radial section.");
    }
 
    return index;
}

Referenced by GeneralAD::source_terms_cellcentered().

Here is the caller graph for this function: