polar.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.


## Polar
# Adrien Crovato
#
# Compute lift and polar curves

import math
import dart.utils as dU
import tbox.utils as tU
from fwk.coloring import ccolors
from .config import Config

class Polar(Config):
    """Utility to compute the polar of a lifting surface

    Attributes
    ----------
    dim : int
        problem dimensions
    format : str
        output format type
    slice : array of float
        arrays of y-coordinates to perform slices (for 3D only)
    tag : int
        index of physical tag to slice (see gmsh)
    mach : float
        freestream Mach number
    alphas : array of float
        range of angles of attack to compute
    Cls : array of float
        lift coefficients for the different AoAs
    Cds : array of float
        drag coefficients for the different AoAs
    Cms : array of float
        moment coefficients for the different AoAs
    """
    def __init__(self, p):
        """Setup and configure

        Parameters
        ----------
        p : dict
            dictionary of parameters to configure DART
        """
        # save some parameters for later use
        self.dim = p['Dim']
        self.format = p['Format']
        try:
            self.slice = p['Slice']
            self.tag = p['TagId']
        except:
            self.slice = None
            self.tag = None
        self.mach = p['M_inf']
        self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step'])
        # basic init
        Config.__init__(self, p)

    def run(self):
        """Compute the polar
        """
        # initialize the loop
        ac = 0 if self.alphas[0] == self.alphas[-1] else 1
        self.Cls = []
        self.Cds = []
        self.Cms = []
        print(ccolors.ANSI_BLUE + 'Sweeping AoA from', self.alphas[0], 'to', self.alphas[-1], ccolors.ANSI_RESET)
        for alpha in self.alphas:
            # define current angle of attack
            alpha = alpha*math.pi/180
            # update problem
            self.pbl.update(alpha)
            # run the solver and save the results
            print(ccolors.ANSI_BLUE + 'Running the solver for', alpha*180/math.pi, 'degrees', ccolors.ANSI_RESET)
            self.tms["solver"].start()
            self.solver.run()
            self.tms["solver"].stop()
            self.solver.save(self.mshWriter, ac)
            print('\n')
            # extract Cp
            if self.dim == 2:
                Cp = dU.extract(self.bnd.groups[0].tag.elems, self.solver.Cp)
                tU.write(Cp, 'Cp_airfoil_'+str(ac)+'.dat', '%1.5e', ', ', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, z, Cp', '')
            elif self.dim == 3 and self.format == 'vtk' and self.slice:
                dU.writeSlices(self.msh.name, self.slice, self.tag, ac)
            # extract force coefficients
            self.Cls.append(self.solver.Cl)
            self.Cds.append(self.solver.Cd)
            self.Cms.append(self.solver.Cm)
            # end loop
            ac+=1

    def disp(self):
        """Display the results and draw the polar
        """
        # display results
        print(ccolors.ANSI_BLUE + 'Airfoil polar' + ccolors.ANSI_RESET)
        print("       M    alpha       Cl       Cd       Cm")
        i = 0
        while i < len(self.alphas):
            print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}".format(self.mach, self.alphas[i], self.Cls[i], self.Cds[i], self.Cms[i]))
            i = i+1
        print('\n')
        # display timers
        self.tms["total"].stop()
        print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET)
        print(self.tms)
        # plot results
        if self.alphas[0] != self.alphas[-1]:
            tU.plot(self.alphas, self.Cls, "alpha", "Cl", "")
            tU.plot(self.alphas, self.Cds, "alpha", "Cd", "")
            tU.plot(self.alphas, self.Cms, "alpha", "Cm", "")