| skeptik |
January 24, 2013 05:43 |
NACA 5-digits airfoil geometry and mesh
2 Attachment(s)
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).
Attachment 18428
In SALOME I used smth like this:
Code:
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
Code:
#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.
Attachment 18429
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? |
