# Boundary Layer¶

## Description¶

Analyze wall-resolved boundary layers.

## Construction¶

import antares
myt = antares.Treatment('bl')


## Parameters¶

• base: antares.Base

The input base.

• coordinates: list(str)

The variable names that define the set of coordinates.

• bl_type: str, default= aerodynamic

Type of boundary layer:

• ‘aerodynamic’: aerodynamic boundary layer

• ‘thermal’: thermal boundary layer

• edge_crit: int, default= 7

Criteria used for the boundary layer edge detection. Add the integer value associated to each criterion to get the chosen combination.

1. $$dV/ds < 1\%$$

1. Stock & Haase

1. isentropic velocity

1. isentropic Mach number

if edge_crit includes 8, then the edge results from the closest point from the wall above all criteria used. Otherwise, the edge comes from the average of all edges computed independent from each criteria.

e.g.: ‘edge_crit=5’ means a mix between the criteria 1 and 4.

• wall_type: str, default= resolved

Type of wall resolution (‘resolved’ or ‘modeled’).

• mesh_type: str in [2D, 2.5D, 3D], default= 3D

Dimensionality of the mesh.

• max_bl_pts: int, default= 100

Maximum number of points in the boundary layer.

The number of points inside the boundary layer = number of calculated boundary layer points + added_bl_pts

• smoothing_cycles: int, default= 0

Number of smoothing cycles for boundary layer results.

• families: dict(str: int) or list(str)

Surface parts on which the boundary layer is to be computed.

If the value is a dictionary: The key value is a family name. This family name must be refered by the boundary condition objects. The boundary condition associated to this family name will be processed. The value associated to a key is a marker of the corresponding family. Two families may have the same marker.

As an example,

t[‘families’] = {‘WingExt’: 3, ‘Flap’: 4, ‘WingInt’: 3}

Both ‘WingExt’ and ‘WingInt’ families have the same marker. Then, all corresponding mesh points will be processed alltogether. They will be concatenated in the output base.

If the value is a list, then a different marker will be assigned to each item of the list.

As an example,

t[‘families’] = [‘WingExt’, ‘Flap’, ‘WingInt’]

Value 0 is special, and applied to all points that do not belong to an input family. So if you give a marker 0 to a family, then you will get all boundary conditions in the ouput base.

• only_profiles: bool, default= False

Only the profiles selected by profile_points will be calculated.

• profile_points: list(list(float)), default= []

Profiles will be plotted according to the surface points. If you give the coordinates of a point in the mesh, the boundary layer profile issued from the nearest point to the surface will be computed.

• bl_parameters: dict, default= {}

Dictionary with keys in [“2D_offset_vector”, “reference_Mach_number”, “reference_velocity”, “reference_pressure”, “reference_density”, “reference_temperature”, “Prandtl_number”, “gas_constant_gamma”, “gas_constant_R”]

• offsurf: int, default= -1

Distance from the wall. To get values on a surface lying at this distance from the wall.

• eps_Mis: float, default= 0.01

Threshold on the isentropic Mach number for edge_crit = 8.

• variables: list(str), default= [“density”, “x_velocity”, “y_velocity”, “z_velocity”, “pressure”]

Names of input primitive variables (density, velocity components, static pressure).

• nb_zone_passingthrough: int, default= 50

Maximum number of zone a vector could pass through. When the maximum is reached, the process will stop to propagate the vector (will be cut before the end).

• keep_interpolation_data: bool, default= False

Allow intermediary results (interpolation coefficients, C objects to remain instantiated in the treatment, in order not to be recomputed in a further call of the same treatement. Warning, the mesh must remain the same.

• instant_selection: int or str, default= 0

Instant name (or index) to execute in the base.

• change_input_data_type: bool, default= False

C routines need to have a specific type for the provided variables. Mostly float64 for float, and int32 for int. And C ordered tables. If the provided base does not have the correct type and ordering, a retyped/ordered copy is performed. To save memory, this copy could replace the input variable. This will avoid to keep both tables instantiated at the same time, but may increase the memory used in the initial base (as new copies will mostly use a bigger type size). This copy will only apply to coordinates and retrieved variables.

## Preconditions¶

Zones must be unstructured.

Input flow field values must be dimensionalized. The standard variables named density, x_velocity, y_velocity, z_velocity and pressure must be given in instants.

If edge_crit >= 8, then the variables total_temperature and Mach_number must also be given.

## Postconditions¶

Three output bases.

If you give the coordinates of a point in the mesh, it will compute the boundary layer profile issued from the nearest point to the surface.

A profile is a line issued from a point (mesh point) on a surface that is perpendicular to this surface. The points along the boundary layer profile are computed from the volume grid points. A point along the boundary layer profile is the result of the intersection between the normal and a face of a volume element.

## Example¶

import antares
myt = antares.Treatment('Bl')
myt['base'] = base
myt['coordinates'] = ['points_xc', 'points_yc', 'points_zc']
myt['families'] = {'Slat': 2, 'Flap': 4, 'Main': 3}
myt['mesh_type'] = '2.5D'
myt['max_bl_pts'] = 50
myt['smoothing_cycles'] = 0
myt['only_profiles'] = False
myt['bl_type'] = 'aerodynamic' # 'aerodynamic' 'thermal'
myt['offsurf'] = 0.5
myt['profile_points'] = [[1.0 , 1.0 , 0.0],]
t['bl_parameters'] = {'2D_offset_vector': 2, 'gas_constant_R': 287.0,
'reference_Mach_number': 0.1715, 'reference_velocity': 36.4249853252,
'reference_pressure': 101325.0, 'reference_density': 3.14466310931,
'reference_temperature': 112.269190122, 'gas_constant_gamma': 1.4,
'Prandtl_number': 0.72, }
res_base, profbase, offsurfbase = t.execute()


## Main functions¶

class antares.treatment.codespecific.boundarylayer.TreatmentBl.TreatmentBl

This class is used to analyze boundary layers.

The treatment is divided in many parts. Some of them are developed in C language to get CPU performance. C routines are made available in python with the module ‘ctypes’.

1. read results and user parameters

2. build a volume element connectivity structure

3. compute the surface normals

4. interpolate the flow field values to the surface normals

5. compute the boundary layer edge see file calculate_boundary_layer_edge.c

6. define the local (s, n) coordinate system

7. compute the boundary layer values in the (s, n) coordinate system see file calculate_boundary_layer.c

8. modify the boundary layer edge

9. define the local (s, n) coordinate system

10. compute the boundary layer values in the (s, n) coordinate system

static asvoid(arr)

View the array as dtype np.void (bytes).

Based on https://stackoverflow.com/a/16973510/190597(Jaime, 2013-06) The items along the last axis are viewed as one value. This allows comparisons to be performed which treat entire rows as one value.

compute_constant_data(base)

Compute normal vectors and interpolation coefficients

compute_interzonetopo(base, cstzone, zone_loc_bnds)

Create additional boundary to consider elements coming from other zone when doing normal vector calculation.

execute()

Execute the treatment.

Returns

the base containing the results

Return type

Base

initialization()

Set the constant part.

Only needed once (both global then cstzone)

merge_prism_layer(nb_zone_passingthrough)

Merge prism layer layer between additional vector and main one.

### Notations:¶

Indices $$\displaystyle \delta$$ and $$\displaystyle \infty$$ are used respectively for the boundary layer edge and for the free stream conditions.

‘s’ denotes the coordinate of the streamwise direction, and ‘n’ the coordinate of the normal to the streamwise direction. ‘s’ and ‘n’ form the (s, n) streamline-oriented coordinate system. This coordinate system is curvilinear, orthogonal and attached to the surface.

### Boundary layer profile (normal to the wall):¶

If asked, the profiles are output in a formated ascii file containing:

• a first section (header) with free stream conditions (given as parameters):

 Ma Mach number $$\displaystyle M_\infty$$ P Pressure $$\displaystyle p_\infty$$ Rho Density $$\displaystyle \rho_\infty$$ T Temperature V Velocity modulus $$\displaystyle V_\infty$$ Cp Pressure coefficient Ptot Total pressure Pdyn Dynamic pressure
• a section that is repeated as long as there are profiles:

 Profile number requested Requested coordinates of the profile origin real True coordinates of the profile origin. Profiles are computed at surface grid points Normal vector Normal unit vector at the profile origin

then come the integral values:

 Del Boundary layer thickness Dels Displacement thickness [s-dir] Deln Displacement thickness [n-dir] Delt Thermal boundary layer thickness approximation Thtss Momentum thickness [s-dir] Thtnn Momentum thickness [n-dir] H12 Incompressible shape factor

then the wall values:

 FirstLayer height a P Pressure cp Pressure coefficient cfsI, cfnI Skin friction coefficients [s, n dirs], $$\displaystyle V_\infty$$ cfrI Resulting skin friction coefficient, $$\displaystyle V_\infty$$ cfsD, cfnD Skin friction coefficients [s, n dirs], $$\displaystyle V_\delta$$ cfrD Resulting skin friction coefficient, $$\displaystyle V_\delta$$ Tw Temperature Rhow Density muew Dynamic viscosity nuew Kinematic viscosity Taus Shear stress component [s-dir] Taun Shear stress component [n-dir] Taur Resulting shear stress

then the edge values:

 Prism height Total height of the prism layer Vdel Velocity modulus Ma Mach number ptot Total pressure cp Pressure coefficient

finally, come different values along the profile:

 S Distance to the wall S/S_edge adimensionalized distance to the wall Vx, Vy, Vz Cartesian velocity components Vres/V_inf Velocity modulus normalized by the free stream velocity modulus Vs/Ve s normalized velocity component Vn/Ve n normalized velocity component beta Skew angle log10(Y+) logarithm of the nondimensional wall distance Vr/V Velocity modulus normalized by the shear stress velocity modulus Vr/V*_theo Velocity modulus normalized by the shear stress velocity modulus based on a flow over a flat plate (theory) Rho Density cp_static Static pressure coefficient cp_dyn Dynamic pressure coefficient Ptot Total pressure PtotLoss Total pressure loss T Temperature Ttot Total temperature Trec to be defined Vi Isentropic velocity modulus Vi-Vr Isentropic velocity modulus - velocity modulus MassFlow Mass flow tke Turbulent Kinetic Energy mu_t/mu_l ratio of viscosities mu_t Eddy viscosity Mach Mach number Mach_is Isentropic Mach Number

### Boundary layer surface:¶

A surface contains in the mesh can be divided in many parts. A marker is a set of triangular and quadrilateral faces.

All the surfaces of a computation may be defined as families (or markers).

It can return a 2D field on the surface (or selected markers) with the following boundary layer quantities:

 X, Y, Z Cartesian coordinates of the surface point VecNx, VecNy, VecNz Cartesian coordinates of the surface normal vector VxSurf, VySurf, VzSurf Cartesian velocity components close to the wall [first grid layer] for wall streamlines VxEdge, VyEdge, VzEdge Cartesian velocity components at the boundary layer edge for inviscid streamlines Vs, Vn s, n velocity components close to the wall beta_w Profile skew angle between wall and boundary layer edge in the s, n coordinate system Y+ Nondimensional wall distance delta Boundary layer thickness, $$\displaystyle \delta$$ deltas Displacement thickness [s-dir] based on the velocity at the boundary layer edge, $$\displaystyle \delta^{*}_s$$ deltan Displacement thickness [n-dir] based on the velocity at the boundary layer edge, $$\displaystyle \delta^{*}_n$$ deltat Thermal boundary layer thickness, $$\displaystyle \delta_T$$ (approximation) deltaes Kinetic energy thickness [s-dir] deltaen Kinetic energy thickness [n-dir]; thetass Momentum thickness [s-dir], $$\displaystyle \theta_{ss}$$ thetann Momentum thickness [n-dir], $$\displaystyle \theta_{nn}$$ H12 Incompressible shape factor. A shape factor is used in boudary layer flow to determine the nature of the flow. The higher the value of H, the stronger the adverse pressure gradient. A large shape factor is an indicator of a boundary layer near separation. A high adverse pressure gradient can greatly reduce the Reynolds number at which transition into turbulence may occur. Conventionally, H = 2.59 (Blasius boundary layer) is typical of laminar flows, while H = 1.3 - 1.4 is typical of turbulent flows. Ptot Total pressure cp Pressure coefficient Ma_edge Mach number at the boundary layer edge CfsInf, CfnInf Skin friction coefficients [s, n dirs], $$\displaystyle V_\infty$$ CfrInf Resulting skin friction coefficient, $$\displaystyle V_\infty$$ CfsDel, CfnDel Skin friction coefficient [s, n dirs], $$\displaystyle V_\delta$$ CfrDel Resulting skin friction coefficient, $$\displaystyle V_\delta$$ maxEddyVisc Maximum eddy viscosity along the profile T_wall Wall temperature T_edge Temperature at the boundary layer edge T_rec Temperature at the boundary layer edge, TO BE DETAILED profTag Profile tag pointMeshIndex Index of the point in the array of grid points heightPrismLayers Total height of the prism layer nbPrismLayers Number of prism layers BlThickOutPrismThick, (delta - heightPrismLayers)/delta*100 firstLayerHeight Height of the first prism layer ptsExceeded 1 if the number of profile points is close to the maximum number of boundary layer points given as input, NOT YET AVAILABLE profFound 100 tells that the boundary layer edge is found Lam/Turb Transition Laminar/Turbulent nbBlPts Number of points in the boundary layer

The mathematical definitions are given below:

 Displacement thicknesses (based on the velocity at the boundary layer edge): $$\displaystyle \delta^{*}_s = \int_0^\delta (1-\frac{\rho}{\rho_\delta}\frac{V_s}{V_\delta}) ds$$ (deltas), $$\displaystyle \delta^{*}_n = -\int_0^\delta (1-\frac{\rho}{\rho_\delta}\frac{V_n}{V_\delta}) ds$$ (deltan) Momentum thicknesses: $$\displaystyle \theta_{ss} = \int_0^\delta \frac{\rho}{\rho_\delta}\frac{V_s}{V_\delta}(1-\frac{V_s}{V_\delta}) ds$$ (thetass), $$\displaystyle \theta_{nn} = \int_0^\delta \frac{\rho}{\rho_\delta}\frac{V_n}{V_\delta}(1-\frac{V_n}{V_\delta}) ds$$ (thetann) Kinetic energy thicknesses: $$\displaystyle \theta_{es} = \int_0^\delta \frac{\rho}{\rho_\delta}\frac{V_s}{V_\delta}(1-\frac{V_s^2}{V_\delta^2}) ds$$ (deltaes), $$\displaystyle \theta_{en} = \int_0^\delta \frac{\rho}{\rho_\delta}\frac{V_n}{V_\delta}(1-\frac{V_n^2}{V_\delta^2}) ds$$ (deltaen) Shape factors: $$\displaystyle H12 = \frac{\delta^{*}_s}{\theta_{ss}}$$, $$\displaystyle H22 = \frac{\delta^{*}_n}{\theta_{nn}}$$ Pressure coefficients: $$\displaystyle Cp = \frac{p - p_{\infty}}{\frac{1}{2}\rho_{\infty} V_{\infty}^2}$$ (cp_static), $$\displaystyle Kp = \frac{p_{tot} - p_{\infty}}{\frac{1}{2}\rho_{\infty} V_{\infty}^2}$$ (cp_dyn) Skin friction coefficients: $$\displaystyle Cf_{\delta} = (\mu_w + \mu_t) \frac{{\frac{\partial V}{\partial w}}_{|_w}}{\frac{1}{2}\rho_w V_\delta^2}$$, $$\displaystyle Cf_{\infty} = (\mu_w + \mu_t) \frac{{\frac{\partial V}{\partial w}}_{|_w}}{\frac{1}{2}\rho_\infty V_\infty^2}$$ $$\displaystyle Y^+ = Y \frac{V^\star}{\nu}$$ Velocity normalized by the shear stress velocity: $$\displaystyle V^+ = \frac{V}{V^\star}$$ Shear stress at the wall: $$\displaystyle \tau_w = \mu_w \frac{\partial u}{\partial y}$$ Shear stress velocity (friction velocity) at the wall: $$\displaystyle V^\star = \sqrt{\frac{\tau_w}{\rho_w}}$$ Total pressure loss: $$\displaystyle \frac{\Delta P_t}{P_t} = \frac{P_t - {P_t}_\infty}{{P_t}_\infty}$$ (PtotLoss) Skew angles: $$\displaystyle \beta = \gamma - \gamma_\delta$$, $$\displaystyle \beta_w = \gamma_w - \gamma_\delta$$

Users must keep in mind the following restrictions when analyzing their results.

## Regions:¶

There are two types of regions where an exact boundary layer edge can not be found:

1. The velocity component Vs becomes very small or zero:
1. Stagnation point regions

At a stagnation point, Vs=0. Then, there is no boundary layer thickness. Near stagnation points, it is also difficult to find a boundary layer thickness accurately.

2. Trailing edge regions

In thick trailing edge regions, Vs becomes very small.

2. The surface normals do not leave the boundary layer:
1. Intersection regions (wing-body, body-tail, etc)
1. The surface normal does not leak from the boundary layer.

2. The surface normal passes the boundary layer from another component first.

3. The surface normal hits another surface.

No accurate boundary layer edge can be found in these cases.

## Thickness:¶

The best situation is when the boundary layer edge lies in the prism layer. In this case, the difference between the true boundary layer edge and the grid point laying on the profile that is upper (or lower) than this true edge is small.

If the boundary layer edge lies outside the prism layer, then the diffence may be more important. Of course, that depends on the mesh quality, but, in general, mesh cell ratio increase quite fast outside the prism layer. In this case, the computed boundary layer thickness can be overestimated, and the step between two successive points on the profile can be large. Then, it is necessary to check the displacement (or momentum) thickness to check the lower pointsand the upper point (respective to the boundary layer edge).

## Boundary layer edge:¶

The edge of the boundary layer is detected with the following multistage algorithm.

A first guess of the boundary layer edge is made with four different methods (see input parameter edge_crit) (function calculate_boundary_layer_edge()):

1. the first derivative of the boundary layer profile velocity is smaller than the 1% and the step between one derivative and the previous one is smaller than the 0.01%. (function edge_velocity_slope())

2. $$\displaystyle \delta = \epsilon y_{max}$$ with $$\displaystyle y_{max}$$ is the wall distance for which the diagnostic function $$\displaystyle y^a |\frac{\partial u}{\partial y}|^b$$ is maximum. The first 10 points are omitted from this maximum search. The constants for a turbulent boundary layer are evaluated from Coles velocity profiles: $$\displaystyle a=1, b=1, \epsilon = 1.936$$. The constants for a laminar boundary layer come from quasi-similar compressible solutions including heat transfer effects: $$\displaystyle a=3.9, b=1, \epsilon = 1.294$$. The reference is [STOCKHAASE].

1. the relative difference between the isentropic velocity and the velocity of the profile is smaller than a threshold.

$$\displaystyle V_r = \sqrt{V_x^2 + V_y^2 + V_z^2}$$ and $$\displaystyle V_i = V_{\infty} \sqrt{\frac{2}{(\gamma-1)*M_{\infty}^2} *(1-(0.5*\frac{2*(P-P_{\infty})}{\rho_{\infty}*V_{\infty}^2}*\gamma*M_{\infty}^2 + 1)^{(\gamma-1)/\gamma}) + 1}$$ and

The edge is the first point giving $$\displaystyle \frac{| V_r - V_i|}{Vr} < \epsilon_{V_i}$$ with $$\displaystyle \epsilon_{V_i} = 2.5 \, 10^{-3}$$

If none, then the edge is given by $$\displaystyle \min{\frac{| V_r - V_i|}{Vr}}$$

1. hypothesis: the isentropic Mach number is constant in the boundary layer. Comparing with the Mach number gives the edge. $$\displaystyle P_{t_{is}} = P_{ref} (T_t/T_{ref})^{\gamma/(\gamma-1)}$$, $$\displaystyle T_{t_{is}} = (P_{t_{is}}/P)^{(\gamma-1)/\gamma}$$, $$\displaystyle M_{is} = \sqrt{2 |T_{t_{is}}-1| / (\gamma-1)}$$

$$\displaystyle \delta M = M_{is} - M$$, $$\displaystyle \delta M_{min} = \min(\delta M)$$, $$\displaystyle M_{is_{cor}} = M_{is} - \delta M_{min}$$.

The edge is the first point giving $$\displaystyle | M - M_{is_{cor}}| / M_{is_{cor}} < \epsilon_{M_{is}}$$ with $$\displaystyle \epsilon_{M_{is}} = 0.01$$ by default.

The boundary layer edge is at the position where the first derivative of “V+” becomes approximately zero.

Start defining a range for the edge search. This range, stored in “points”, will be the 160% of the first guess of the boundary layer edge.

The first step is to search the first maximum of the first derivative of “V+” and then, starting from this maximum, search where the first derivative is smaller than two. Once this first edge is found, the next maximum is compared to the first one. If it is bigger, search again the boundary layer edge starting from this new maximum.

When the final edge is found, perform two checks. If the thickness of the boundary layer is greater than two, move the edge back to the last maximun or minimum of the first derivate of “V+” before the distance to the wall is two.

The second check is done using the isentropic velocity. Search for the point where the isentropic velocity meets the profile velocity but only until three points below the current boundary layer edge. If the point is found, average it with the current edge to get the final and corrected boundary layer edge.

## Interpolation:¶

The interpolation of the CFD values to the normal (profile) is done with a polynomial expansion of the form $$\displaystyle G^i = a + b x + c y$$, which represents a linear variation of G in both x and y directions within a triangular element. The three constants a, b, and c are computed with the three triangle points. more details can be found in [CHUNG].

Note

Sutherland law for molecular viscosity: $$\displaystyle \mu = \frac{C_1 \sqrt{T_w}}{1+\frac{S}{T_w}}$$ with $$\displaystyle C_1 = 1.46 10^{-6} kg.m^{-1}.s^{-1}.K^{-1/2}$$ and $$\displaystyle S = 112. K$$ (functions Skinfric() of calculate_boundary_layer.c and bl_profile(), clear_values() of calculate_boundary_layer_edge.c)

Note

Function $$\displaystyle U^+$$ versus $$\displaystyle y^+$$:

linear part: $$\displaystyle U^+ = y^+$$ for $$\displaystyle y^+ < 11.63$$

log part: $$\displaystyle U^+ = 5.75 * \log(y^+) + 5.5$$ for $$\displaystyle y^+ > 11.63$$

Note

if not given, T_wall is computed with the equation of state of perfect gas: $$\displaystyle T_w = \frac{\gamma}{(\gamma-1)C_p} \frac{P}{\rho}$$ where $$\displaystyle C_p = 1,003.41 J.kg^{-1}.K$$ is here the specific heat at constant pressure (at T = 250 K), and $$\displaystyle \gamma = \frac{C_p}{C_v}$$ given as input, with $$\displaystyle C_v$$ the specific heat at constant volume

Note

The boundary layer thickness $$\displaystyle (\delta)$$ is the average of a numerical value (distance of the boundary layer edge from the wall) and a correlation value coming from the shape factor (function Thickness() of calculate_boundary_layer.c).

$$\displaystyle \delta = \frac{\delta_{sf}+\delta_{num}}{2}$$

$$\displaystyle \delta_{sf} = \delta^{*I} * (\frac{\theta_s^{I}}{\delta_s^{*I}} H_{1}+1)$$ where the superscript I means that the incompressible formula is used and with $$\displaystyle H_{1}$$ the mass flow shape factor computed with a correlation of Green [GREEN].

Note

$$\displaystyle \delta_T$$ (deltat) is the E. Polhausen approximation of the thermal boundary layer thickness for Prandtl number greater than 0.6

$$\displaystyle \delta_T = \frac{\delta}{{Pr}^{1/3}}$$

Note

Numerical values that are not available as inputs are computed with:

$$\displaystyle M_a = \frac{V}{\sqrt(\gamma*R*T)}$$

$$\displaystyle P_{tot} = (1+.2*M_a*M_a)^{3.5} * P$$ (so $$\displaystyle \gamma = 1.4$$)

GREEN

PROCUREMENT EXECUTIVE MINISTRY OF DEFENCE AERONAUTICAL RESEARCH COUNCIL REPORTS AND MEMORANDA no 3791 Prediction of Turbulent Boundary Layers and Wakes in Compressible Flow by a Lag-Entrainment Method By J. E. GREEN, D. J. WEEKS AND J. W. F. BRODMAN, Aerodynamics Dept., R.A.E, Farnborough

STOCKHAASE

Hans W. Stock and Werner Haase. “Feasibility Study of e^N Transition Prediction in Navier-Stokes Methods for Airfoils”, AIAA Journal, Vol. 37, No. 10 (1999), pp. 1187-1196. https://doi.org/10.2514/2.612

CHUNG

Computational Fluid Dynamics. T. J. Chung. Cambridge University Press, 2002. pp 273–275

## Example¶

import os
if not os.path.isdir('OUTPUT'):
os.makedirs('OUTPUT')

import antares

r['filename'] = os.path.join('..', 'data', 'BL', 'EUROLIFT_2D', 'grid1_igs_Local_Sweep.nc')

r['base'] = base
r['filename'] = os.path.join('..', 'data', 'BL', 'EUROLIFT_2D', 'sol_Local_Sweep_ETW1.pval.30000')

t = antares.Treatment('Bl')
t['base'] = base
t['coordinates'] = ['points_xc', 'points_yc', 'points_zc']
t['families'] = {'Slat_2': 2, 'Flap_4': 4, 'Main_3': 3}
t['mesh_type'] = '2.5D'
t['max_bl_pts'] = 50
t['smoothing_cycles'] = 0
t['only_profiles'] = False
t['bl_type'] = 'thermal' # 'aerodynamic' 'thermal'
t['bl_type'] = 'aerodynamic' # 'aerodynamic' 'thermal'
t['offsurf'] = 0.5
t['profile_points'] = [[1.0 , 1.0 , 0.0],
[2.0 , 1.0 , 0.0],
[3.0 , 1.0 , 0.0],
[5.0 , 1.0 , 0.0],
[7.0 , 1.0 , 0.0],
[9.0 , 1.0 , 0.0],
[10.0, 1.0,  0.0],
[15.0, 1.0,  0.0]
]
t['bl_parameters'] = {'2D_offset_vector': 2, 'gas_constant_R': 287.0,
'reference_Mach_number': 0.1715, 'reference_velocity': 36.4249853252,
'reference_pressure': 101325.0, 'reference_density': 3.14466310931,
'reference_temperature': 112.269190122, 'gas_constant_gamma': 1.4,
'Prandtl_number': 0.72, }
res_base, profbase, offsurfbase = t.execute()

w = antares.Writer('bin_tp')
w['base'] = res_base
w['filename'] = os.path.join('OUTPUT', 'ex_blsurf.plt')
w.dump()

if offsurfbase:
w = antares.Writer('bin_tp')
w['base'] = offsurfbase
w['filename'] = os.path.join('OUTPUT', 'ex_offblsurf.plt')
w.dump()