main.cpp 2.73 KB
Newer Older
1
// This code simulates a potential flow (subsonic flow) around a horizontal NACA0012 airfoil.
Alexandre Halbach's avatar
Alexandre Halbach committed
2
// The problem is nonlinear because the air density depends on the air speed.
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27


#include "sparselizardbase.h"


using namespace mathop;

void sparselizard(void)
{	
    // The domain regions as defined in 'airfoil2D.geo':
    int air = 1, airfoil = 2, downstream = 3, upstream = 4;
    
    // Load the 'airfoil2D' mesh: 
    mesh mymesh("airfoil2D.msh");
    
    // Define the whole outer boundary:
    int gammaouter = regionunion({upstream, downstream});
    
    // Define the velocity potential field 'phi' with standard nodal shape functions ("h1").
    // grad(phi) is the fluid velocity. Field x is the x coordinate field.
    field phi("h1"), x("x");
    
    // Specific weight of air under some circumstances:
    double gamma = 1.4;
    
Alexandre Halbach's avatar
Alexandre Halbach committed
28
    // Define the air density 'rho' and the Mach number:
29
    expression rho = pow(1 + (gamma-1)/2.0 * 0.7 * 0.7 * (1-grad(phi)*grad(phi)), 1.0/(gamma-1));
Alexandre Halbach's avatar
Alexandre Halbach committed
30
    expression machnumber = sqrt(grad(phi)*grad(phi)) / sqrt(1.0/(0.7*0.7) + 0.5*(gamma-1) * (1-grad(phi)*grad(phi)));
31 32 33 34 35 36 37 38 39 40
    
    // We suppose outside the air domain a uniform speed of 1 with direction left to right.
    // Since grad(phi) is the speed we have that grad(x) indeed gives us what we want.
    phi.setconstraint(gammaouter, x);
    
    // Define the potential flow formulation:
    formulation pf;
    
    // On the airfoil boundary the default Neumann condition dphi/dnormal = 0 
    // automatically ensures that there is no fluid entering the airfoil.
Alexandre Halbach's avatar
Alexandre Halbach committed
41 42 43 44 45 46 47 48 49
    // We thus do not need anything else than this:
    pf += integral(air, rho * grad(dof(phi)) * grad(tf(phi)) );
    
    // Start the nonlinear iteration with an all zero guess:
    vec sol(pf);
    
    double relres = 1;
    while (relres > 1e-7)
    {
Alexandre Halbach's avatar
Alexandre Halbach committed
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
        // Generate the formulation:
        pf.generate();
        // Get A and b to solve Ax = b:
        mat A = pf.A();
        vec b = pf.b();
        
        // Compute a relative residual:
        relres = (b - A*sol).norm() / b.norm();
        
        // Solve Ax = b:
        sol = solve(A, b);
        
        // Transfer the data from the solution vector to field 'phi' on the whole 'air' region:
        phi.setdata(air, sol);
        
        std::cout << "Current iteration has relative residual: " << relres << std::endl;
    }
    
Alexandre Halbach's avatar
Alexandre Halbach committed
68
    // Write the fluid speed (i.e. grad(phi)) and the Mach number:
Alexandre Halbach's avatar
Alexandre Halbach committed
69
    grad(phi).write(air, "flowspeed.pos");
Alexandre Halbach's avatar
Alexandre Halbach committed
70
    machnumber.write(air, "machnumber.pos");
Alexandre Halbach's avatar
Alexandre Halbach committed
71 72
    
    // Code validation line. Can be removed.
Alexandre Halbach's avatar
Alexandre Halbach committed
73
    std::cout << (machnumber.integrate(air, 3) < 62.4149 && machnumber.integrate(air, 3) > 62.4145);
74 75 76 77 78 79 80 81 82 83 84 85 86 87
}

int main(void)
{	
    SlepcInitialize(0,{},0,0);

    sparselizard();

    SlepcFinalize();

    return 0;
}