El problema “Pressure over airfoil surface with OpenFOAM” resuelto. Diccionario sampleDict.

En un principio pudo parecer más sencillo de lo que es.

Paraview nos permite realizar muchas tareas de post-procesado una vez tenemos todo calculado, sin embargo esto no es muy cómodo cuando cuando lo que se pretende es extraer datos para trabajar con ellos con el objetivo de evolucionar el modelo como es mi caso. De acuerdo, se pueden usar macros e incluso se pueden escribir scripts de post-procesado mediante Python para paraview… pero creo que esto tiene más ventajas y es más asequible para todo aquel ajeno a la programación.

OpenFOAM nos ofrece utilidades de post-procesado, el problema es que están documentadas al estilo OpenFOAM, es decir, dan pocas pistas sobre cómo usarlas.

Principalmente usaremos el comando “sample” durante este post porque no quisiera emplear software de terceros para mi programa si es posible (link a la documentación de sample). En cualquier caso vosotros podéis echarle un ojo a swak4Foam que permite hacer bastantes más cosas.

Vamos con sample. Para emplear el comando sample necesitamos un diccionario, esto es, un documento llamado sampleDict dentro de la carpeta system del caso. Una vez que uno ve la pinta que tiene el diccionario cuando encuentra un fichero de ejemplo entre todo el código de OpenFOAM las cosas marchan mucho mejor. El fichero pinta tal que así:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      sampleDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Set output format : choice of
//      xmgr
//      jplot
//      gnuplot
//      raw
//      vtk
//      csv
setFormat raw;

// Surface output format. Choice of
//      null        : suppress output
//      ensight     : Ensight Gold format, one field per case file
//      foamFile    : separate points, faces and values file
//      dx          : DX scalar or vector format
//      vtk         : VTK ascii format
//      raw         : x y z value format for use with e.g. gnuplot 'splot'.
//
// Note:
// other formats such as obj, stl, etc can also be written (by proxy)
// but without any values!
surfaceFormat vtk;

// optionally define extra controls for the output formats
formatOptions
{
    ensight
    {
        format  ascii;
    }
}

// interpolationScheme. choice of
//      cell          : use cell-centre value only; constant over cells (default)
//      cellPoint     : use cell-centre and vertex values
//      cellPointFace : use cell-centre, vertex and face values.
//      pointMVC      : use point values only (Mean Value Coordinates)
//      cellPatchConstrained : use cell-centre except on boundary faces where
//        it uses the boundary value. For use with e.g. patchCloudSet.
// 1] vertex values determined from neighbouring cell-centre values
// 2] face values determined using the current face interpolation scheme
//    for the field (linear, gamma, etc.)
interpolationScheme cellPoint;

// Fields to sample.
fields
(
    p
    U
);

// Set sampling definition: choice of
//      uniform             evenly distributed points on line
//      face                one point per face intersection
//      midPoint            one point per cell, inbetween two face intersections
//      midPointAndFace     combination of face and midPoint
//
//      polyLine            specified points, not nessecary on line, uses
//                          tracking
//      cloud               specified points, uses findCell
//      triSurfaceMeshPointSet  points of triSurface
//
// axis: how to write point coordinate. Choice of
// - x/y/z: x/y/z coordinate only
// - xyz: three columns
//  (probably does not make sense for anything but raw)
// - distance: distance from start of sampling line (if uses line) or
//             distance from first specified sampling point
//
// type specific:
//      uniform, face, midPoint, midPointAndFace : start and end coordinate
//      uniform: extra number of sampling points
//      polyLine, cloud: list of coordinates
//      patchCloud: list of coordinates and set of patches to look for nearest
sets
(
    lineX1
    {
        type        uniform;
        axis        distance;

        //- cavity. Slightly perturbed so not to align with face or edge.
        start       (0.0201 0.05101 0.00501);
        end         (0.0601 0.05101 0.00501);
        nPoints     10;
    }

    lineX2
    {
        type        face;
        axis        x;

        //- cavity
        start       (0.0001 0.0525 0.00501);
        end         (0.0999 0.0525 0.00501);
    }

    somePoints
    {
        type    cloud;
        axis    xyz;
        points  ((0.049 0.049 0.00501)(0.051 0.049 0.00501));
    }

    somePatchPoints
    {
        // Sample nearest points on selected patches. Looks only up to
        // maxDistance away. Any sampling point not found will get value
        // pTraits::max (usually VGREAT)
        // Use with interpolations:
        // - cell (cell value)
        // - cellPatchConstrained (boundary value)
        // - cellPoint (interpolated boundary value)
        type        patchCloud;
        axis        xyz;
        points      ((0.049 0.099 0.005)(0.051 0.054 0.005));
        maxDistance 0.1;    // maximum distance to search
        patches     (".*Wall.*");
    }
);

He puesto solamente la parte que se ocupa de extraer datos bidimensionales, dejo el resto para otra entrada.

Lo primero que nos encontramos es la opción setFormat, las distintas configuraciones están bastante claras, yo siempre trabajo en raw. Prefiero escribir mis scripts de Gnuplot adaptándome a ese formato.


setFormat raw;

La opción formatOptions no nos hace falta por ahora. Podemos borrar el campo si queremos, no hay problema.

En el campo fields debemos introducir los campos de los cuales pretendemos extraer datos, dejemos por ahora la presión;


fields ( p );

Dentro del campo sets {} se especifican los distintos campos de datos que buscamos extraer.

Ahora tenemos varias opciones en cuanto a la forma de especificar cómo queremos que se distribuyan los puntos de control. LineX1 obtiene los puntos de control dividiendo la recta que especificamos con el punto de inicio y de fin en tantos segmentos como nPoints. Podéis jugar con las distintas opciones y seguramente os deis cuenta de que es difícil extraer datos de una superficie curva.

El problema a la hora de extraer los datos de los puntos que deja la intersección de un plano sobre una superficie es que no hay una manera eficiente de hacerlo (o no la he encontrado). No podemos especificar los puntos como la intersección de un plano y un patch, por lo que lo único que se me ha ocurrido es suministrarle al diccionario los puntos de forma manual, suerte que tengo esas coordenadas calculadas a la hora de hacer la malla…

Tenemos dos opciones para hacer esto, o bien somePoints o bien somePathPoints. Quizás la segunda sea más segura, pero puede no ser muy exacta, debemos de tener cuidado y comprobar los datos. Hay que prestar especial atención a la hora de elegir los métodos de interpolación teniendo en cuenta el mallado.

La opción somePoints nos da cierto margen de error a la hora de calcular nuestros puntos.

En cualquiera de las dos opciones debemos especificar los puntos. Como os digo, yo los tengo calculados porque con ellos genero las mallas. La forma de hacerlo sería entonces algo así..


/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
version     2.0;
format      ascii;
class       dictionary;
object      sampleDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

setFormat raw;
interpolationScheme cellPoint;

fields ( p );

sets
(
    upper
    {
        type        patchCloud;
        axis        xyz;
        points (

      (0.000252 0.000000 -0.002807)
      (0.001007 0.000000 -0.005576)
      (...      ...      -...     ) //No coloco todos los puntos aquí..
      (0.999748 0.000000 -0.000037)
      (1.000000 0.000000  0.000000)
      );
        maxDistance 0.001;    // Distancia maxima de búsqueda
        patches     ("airfoil");
    }
);

Una vez tenemos el diccionario tendremos que resolver el caso con el solver pertinente y, una vez resuelto, ejecutamos desde la terminal el comando “sample”. Si todo va bien nos creará una carpeta llamada sets que contiene distintas carpetas (para los distintos intervalos de iteraciones/tiempo que hemos calculado).
Sólo nos queda ahora hacer un pequeño script para Gnuplot con el que mostrar los datos. Ahí va un ejemplo:

set title 'Pressure over airfoil surface.'
set timestap
set ylabel 'Pressure (Pa)'
set xlabel 'chord (m)'
set grid
plot "sets/500/upper_p.xy" using 1:4 with linespoints pi 10 title 'Upper surface.',\
     "sets/500/lower_p.xy" using 1:4 with linespoints pi 10 title 'Lower surface.'

Veamos primero en paraView el campo de presiones:

Campo de presiones para el perfil en paraView

Y ahora podemos ver la salida:

Salida de Gnuplot para las presiones sobre la superficie del perfil bidimensional.

Parece que no ha ido del todo mal. He probado con un perfil simétrico y las presiones en la superficie superior coinciden con la inferior cuando la corriente es paralela a la cuerda.

He hecho la malla con CMeshFoil y todo esto es parte del post-procesado de datos que necesito para escribir el módulo de optimización de perfiles. Os dejo aquí la malla que cmesfoil ha generado en unos pocos segundos😉

Malla generada por cmeshfoil.

Cuando termine los exámenes intentaré escribir una función que permita extraer datos de la intersección de un plano de corte con un patch.

Espero que os sea de ayuda. Ahora ya tenemos más datos con los que trabajar. Have fun!

Esta entrada fue publicada en Aerodinamica, cfd, CMeshFoil, Gnuplot, OpenFOAM y etiquetada , , , , , . Guarda el enlace permanente.

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s