nonlift.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.


## Compute non-lifting (linear or nonlinear) potential flow around a NACA 0012
# Adrien Crovato
#
# Test the nonlinear shock-capturing capability
#
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.

import dart.utils as floU
import dart.default as floD
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import math

def main():
    # timer
    tms = fwk.Timers()
    tms['total'].start()

    # define flow variables
    alpha = 0*math.pi/180 # must be zero for this testfile
    M_inf = 0.8
    c_ref = 1
    dim = 2

    # mesh the geometry
    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
    tms['msh'].start()
    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.0075, 'msLe' : 0.0075}
    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
    tms['msh'].stop()

    # set the problem and add fluid, initial/boundary conditions
    tms['pre'].start()
    pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil')
    tms['pre'].stop()

    # solve problem
    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
    tms['solver'].start()
    solver = floD.newton(pbl)
    solver.run()
    solver.save(gmshWriter)
    tms['solver'].stop()

    # extract Cp
    Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
    # display results
    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
    print('       M    alpha       Cl       Cd       Cm')
    print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alpha*180/math.pi, solver.Cl, solver.Cd, solver.Cm))

    # display timers
    tms['total'].stop()
    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
    print('CPU statistics')
    print(tms)

    # visualize solution and plot results
    floD.initViewer(pbl)
    floD.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})

    # check results
    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
    tests = CTests()
    if M_inf == 0:
        tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.405, 5e-2))
    elif M_inf == 0.7:
        tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.63, 5e-2))
    elif M_inf == 0.8:
        tests.add(CTest('iteration count', solver.nIt, 13, 3, forceabs=True))
        tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.89, 5e-2))
        tests.add(CTest('Cd', solver.Cd, 0.0058, 5e-4, forceabs=True))
    else:
        raise Exception('Test not defined for this flow')
    tests.add(CTest('Cl', solver.Cl, 0., 1e-2))
    tests.add(CTest('Cm', solver.Cm, 0., 1e-2))
    tests.run()

    # eof
    print('')

if __name__ == "__main__":
    main()