Skip to content
Permalink
Fetching contributors…
Cannot retrieve contributors at this time
87 lines (60 sloc) 2.73 KB
// This code simulates a potential flow (subsonic flow) around a horizontal NACA0012 airfoil.
// The problem is nonlinear because the air density depends on the air speed.
#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;
// Define the air density 'rho' and the Mach number:
expression rho = pow(1 + (gamma-1)/2.0 * 0.7 * 0.7 * (1-grad(phi)*grad(phi)), 1.0/(gamma-1));
expression machnumber = sqrt(grad(phi)*grad(phi)) / sqrt(1.0/(0.7*0.7) + 0.5*(gamma-1) * (1-grad(phi)*grad(phi)));
// 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.
// 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)
{
// 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;
}
// Write the fluid speed (i.e. grad(phi)) and the Mach number:
grad(phi).write(air, "flowspeed.pos");
machnumber.write(air, "machnumber.pos");
// Code validation line. Can be removed.
std::cout << (machnumber.integrate(air, 3) < 62.4149 && machnumber.integrate(air, 3) > 62.4145);
}
int main(void)
{
SlepcInitialize(0,{},0,0);
sparselizard();
SlepcFinalize();
return 0;
}
You can’t perform that action at this time.