[skeptik]

Hi all. I try to make a script which could produce geometry and mesh for different airfoils in SALOME.

First, I made a script using analytical curves and got big gaps between curves (on image).

uncert_snap.jpg

In SALOME I used smth like this:

Yc1="300*(15.957/6*(pow(t/300,3)-3*0.2025*pow(t/300,2)+pow(0.2025,2)*(3-0.2025)*t/300))"
Yc2="300*(15.957/6*pow(0.2025,3)*(1-t/300))"
     
Xup1="300*( t/300-0.6*(0.2969*sqrt(t/300)-0.1260*t/300-0.3516*pow(t/300,2)+0.2843*pow(t/300,3)-0.1015*pow(t/300,4))*sin(atan(15.957/6*(3*pow(t/300,2)-6*0.2025*t/300+pow(0.2025,2)*(3-0.2025)))) )"
Yup1="300*( 15.957/6*(pow(t/300,3)-3*0.2025*pow(t/300,2)+pow(0.2025,2)*(3-0.2025)*t/300)+(0.6)*(0.2969*sqrt(t/300)-0.1260*t/300-0.3516*pow(t/300,2)+0.2843*pow(t/300,3)-0.1015*pow(t/300,4))*cos(atan(15.957/6*(3*pow(t/300,2)-6*0.2025*t/300+pow(0.2025,3)*(3-0.2025)))) )"

Xup2="300*( t/300-0.6*(0.2969*sqrt(t/300)-0.1260*t/300-0.3516*pow(t/300,2)+0.2843*pow(t/300,3)-0.1015*pow(t/300,4))*sin(atan(-15.957/6*pow(0.2025,3))) )"
Yup2="300*( 15.957/6*pow(0.2025,3)*(1-t/300)+(0.6)*(0.2969*sqrt(t/300)-0.1260*t/300-0.3516*pow(t/300,2)+0.2843*pow(t/300,3)-0.1015*pow(t/300,4))*cos(atan(-15.957/6*pow(0.2025,3)) ) )"

Xlo1="300*( t/300+0.6*(0.2969*sqrt(t/300)-0.1260*t/300-0.3516*pow(t/300,2)+0.2843*pow(t/300,3)-0.1015*pow(t/300,4))*sin(atan(15.957/6*(3*pow(t/300,2)-6*0.2025*t/300+pow(0.2025,2)*(3-0.2025)))) )"
Ylo1="300*( 15.957/6*(pow(t/300,3)-3*0.2025*pow(t/300,2)+pow(0.2025,2)*(3-0.2025)*t/300)-(0.6)*(0.2969*sqrt(t/300)-0.1260*t/300-0.3516*pow(t/300,2)+0.2843*pow(t/300,3)-0.1015*pow(t/300,4))*cos(atan(15.957/6*(3*pow(t/300,2)-6*0.2025*t/300+pow(0.2025,3)*(3-0.2025)))) )"

Xlo2="300*( t/300+0.6*(0.2969*sqrt(t/300)-0.1260*t/300-0.3516*pow(t/300,2)+0.2843*pow(t/300,3)-0.1015*pow(t/300,4))*sin(atan(-15.957/6*pow(0.2025,3))) )"
Ylo2="300*( 15.957/6*pow(0.2025,3)*(1-t/300)-(0.12/0.2)*(0.2969*sqrt(t/300)-0.1260*t/300-0.3516*pow(t/300,2)+0.2843*pow(t/300,3)-0.1015*pow(t/300,4))*cos(atan(-15.957/6*pow(0.2025,3))) )"

geompy.init_geom(theStudy)

Curve_1 = geompy.MakeCurveParametric(Xup1, Yup1, "0", 0, p*300, 50, geompy.GEOM.Interpolation, True)
Curve_2 = geompy.MakeCurveParametric(Xup2, Yup2, "0", p*300, 300, 50, geompy.GEOM.Interpolation, True)

Curve_3 = geompy.MakeCurveParametric(Xlo1, Ylo1, "0", 0, p*300, 50, geompy.GEOM.Interpolation, True)
Curve_4 = geompy.MakeCurveParametric(Xlo2, Ylo2, "0", p*300+10, 300, 50, geompy.GEOM.Interpolation, True)

Curve_5 = geompy.MakeCurveParametric("t", Yc1, "0", 0, p*300, 50, geompy.GEOM.Interpolation, True)
Curve_6 = geompy.MakeCurveParametric("t", Yc2, "0", p*300, 300, 50, geompy.GEOM.Interpolation, True)

this is really bad code and I had no results cause this code produce gaps in airfoil.

Next, I read topic from http://www.salome-platform.org/forum...11/thread_2877 , read enough about Python scripting and upgrade my script

#define profile parameters
#chord in millimeters
c = 300
#max thickness
t=0.12
#location of the maximum camber
p=0.15
#define constants NACA 5-digits
a1=0.2969
a2=-0.1260
a3=-0.3516
a4=0.2843
a5=-0.1015
k1=15.957
m=0.2025
r=0

def Yt(x_arg):
    r=t/0.2*(a1*math.sqrt(x_arg)+a2*x_arg+a3*math.pow(x_arg,2)+a4*math.pow(x_arg,3)+a5*math.pow(x_arg,4))
    return r

def Yc(x_arg):
    if x_arg<p:
        r=k1/6*(math.pow(x_arg,3)-3*m*math.pow(x_arg,2)+math.pow(m,2)*(3-m)*x_arg)
    else:
        r=k1/6*math.pow(m,3)*(1-x_arg)
    return r
    
def Xu(x_arg):
    if x_arg<p:
        r=x_arg-Yt(x_arg)*math.sin(math.atan(k1/6*(3*math.pow(x_arg,2)-6*m*x_arg+math.pow(m,2)*(3-m))))
    else:
        r=x_arg-Yt(x_arg)*math.sin(math.atan(-k1/6*math.pow(m,3)))
    return r
    
def Yu(x_arg):
    if x_arg<p:
        r=Yc(x_arg)+Yt(x_arg)*math.cos(math.atan(k1/6*(3*math.pow(x_arg,2)-6*m*x_arg+math.pow(m,2)*(3-m))))
    else:
        r=Yc(x_arg)+Yt(x_arg)*math.cos(math.atan(-k1/6*math.pow(m,3)))
    return r

def Xlo(x_arg):
    if x_arg<p:
        r=x_arg+Yt(x_arg)*math.sin(math.atan(k1/6*(3*math.pow(x_arg,2)-6*m*x_arg+math.pow(m,2)*(3-m))))
    else:
        r=x_arg+Yt(x_arg)*math.sin(math.atan(-k1/6*math.pow(m,3)))
    return r
    
def Ylo(x_arg):
    if x_arg<p:
        r=Yc(x_arg)-Yt(x_arg)*math.cos(math.atan(k1/6*(3*math.pow(x_arg,2)-6*m*x_arg+math.pow(m,2)*(3-m))))
    else:
        r=Yc(x_arg)-Yt(x_arg)*math.cos(math.atan(-k1/6*math.pow(m,3)))
    return r    
    
nPts = 250
vId=0
step=1.0/(nPts-1)
x = 0
z = 0
ptList = []

#In my script i use linear refinement to the tip of airfoil. 
dStart = 2.0/math.pow(nPts,2) 
dStep=0
dX=0

print "Creating %s points"%nPts
print "Step %f"%step

for i in range(nPts):
        v = geompy.MakeVertex(c*dX, c*Yc(dX), z)
        xp = c*dX
        yp = c*Yc(dX)
        vId+=1
        #print "Creating x/c point %f"%x
        #print c*x, c*Yc(x)
        geompy.addToStudy(v, "Vertex_%d"%(vId) )
        x+=step
        dStep+=1
        dX=dX+dStep*dStart
        #print "Different step %f"%dX
        ptList.append(v)
#print "Done"

polyline = geompy.MakePolyline(ptList)
#interpol = geompy.MakeInterpol(ptList)
geompy.addToStudy(polyline, "polyline" )
#geompy.addToStudy(interpol, "interpol" )
print "Camber line from points list... done"

del ptList
del i

print "Creating upper line"
x=0
ptList = []

dStart = 2.0/math.pow(nPts,2)
dStep=0
dX=0

for i in range(nPts):
        v = geompy.MakeVertex(c*Xu(dX), c*Yu(dX), z)
        vId+=1
        print c*Xu(x), c*Yu(x)
        geompy.addToStudy(v, "Vertex_%d"%(vId) )
        x+=step
        dStep+=1
        dX=dX+dStep*dStart
        ptList.append(v)
#print "Done"

polyline_u = geompy.MakePolyline(ptList)
#interpol_u = geompy.MakeInterpol(ptList)
geompy.addToStudy(polyline_u, "polyline_u" )
#geompy.addToStudy(interpol_u, "interpol_u" )
print "Upper line from points list... done"

del ptList
del i

print "Creating lower line"
x=0
ptList = []

dStart = 2.0/math.pow(nPts,2)
dStep=0
dX=0

for i in range(nPts):
        v = geompy.MakeVertex(c*Xlo(dX), c*Ylo(dX), z)
        vId+=1
        print c*Xlo(x), c*Ylo(x)
        geompy.addToStudy(v, "Vertex_%d"%(vId) )
        x+=step
        dStep+=1
        dX=dX+dStep*dStart
        ptList.append(v)
#print "Done"

polyline_lo = geompy.MakePolyline(ptList)
#interpol_lo = geompy.MakeInterpol(ptList)
geompy.addToStudy(polyline_lo, "polyline_lo" )
#geompy.addToStudy(interpol_lo, "interpol_lo" )
print "Lower line from points list... done"

del ptList
del i

geompy.init_geom(theStudy)

polyline_u_vertex_500 = geompy.GetSubShape(polyline_u, [500])
polyline_lo_vertex_500 = geompy.GetSubShape(polyline_lo, [500])
Line_1 = geompy.MakeLineTwoPnt(polyline_u_vertex_500, polyline_lo_vertex_500)
geompy.addToStudy( Line_1, 'Line_1' )

So It's pretty good and i have continious curve.

profile_good.jpg

So i need recommendations about my last script.
May be somebody has an experience mesh generation in SALOME cause next steps will concern on hexa meshing.

I see another way to generate airfoil geometry: read a file with points generated by this java script http://www.ppart.de/aerodynamics/profiles/NACA5.html and read this file in SALOME. May be somebody has a template for this case?

[skeptik]

made a working script for NACA 5-digit profile.

It contains data for NACA 23012, chord length 300 mm.

It creates square domain like in airFoil2D, tutorial of simpleFoam solver.

Next, SMESH generate tetra mesh (2D). Like in geometry there are also present simple variable parameters (maxTetra, minTetra) for tetra generation in the script under #SMESH.

The script is not ideal, it contains a lot of trash (there is blocking structure for hexa but it's not ready yet), but it works. I checked it)

May be it helps somebody to make something on python or make another airfoil mesh.


snapTetraScript.jpg

NACA23012_tetra_quad_domain.zip