edan0314 commented 2 months ago

Comment:

I'm working with GMSH to mesh a rocket nozzle. My teacher who is helping me with my thesis passed me this code, and the following error message is showing:

Exception: Surface 1 cannot be meshed using the transfinite algorithm (non-matching number of nodes on opposite sides 1524 != 1473 or 99 != 99)

Error: Surface 1 cannot be meshed using the transfinite algorithm (non-matching number of nodes on opposite sides 1524 != 1473 or 99 != 99)

I don't have any idea what is going on! I'm just beginning to work with GMSH, and with python things get much more abstract. If anyone knows how to solve it, I'll be eternally grateful!

The real GMSH code begins at line 190, before that, it is just defining the nozzle geometry, and it is all correct. I checked by plotting. ` import gmsh import sys import numpy as np import matplotlib.pyplot as plt

def neg(lst): return [ -i for i in lst ]

def getS(n,r): if r!=1: s=(r-1)/(r**n-1) elif r==1: s=1/n return s

def getPoints(start,stop,n,r): s=getS(n,r) x=np.zeros(n+1) x[0]=start x[1:]=start+(stop-start)snp.cumsum(np.array([r**i for i in range(n)])) return x

def chamber(Rc,step_count,Lc,progression): x = getPoints(0,Lc,step_count,progression)#np.linspace(0,Lc,step_count) y = Rc*np.ones(len(x)) return [x,y]

def chamber_exit_circ(Rc,Rt,Lc,theta_t_i,step_count): chamber_e_l = Lc angles = np.linspace(90, 90 - theta_t_i, step_count+1) Xc = chamber_e_l Yc = Rc -1.5Rt x = 1.5Rtnp.cos(anglesnp.pi/180) + Xc
y = 1.5Rtnp.sin(angles*np.pi/180) + Yc return [x,y]

def throat_inlet(Rt,theta_t_i,step_count,progression): throat_i_l = 1.5Rtnp.sin(theta_t_inp.pi/180) angles = getPoints(270 - theta_t_i, 270, step_count+1,progression) Xc = throat_i_l Yc = 2.5Rt; x = 1.5Rtnp.cos(anglesnp.pi/180) + Xc y = 1.5Rtnp.sin(anglesnp.pi/180) + Yc return [x,y]

def throat_exit(Rt,theta_N,step_count): throat_e_l = 0.382Rtnp.sin(theta_Nnp.pi/180) angles = np.linspace(270 , 270 + theta_N ,step_count+1) Xc = throat_e_l Yc = 1.382Rt x = 0.382Rtnp.cos(anglesnp.pi/180) + Xc y = 0.382Rtnp.sin(anglesnp.pi/180) + Yc return[x,y]

def bell_curve(Rt,Re,theta_N,theta_e,Ln_ratio,step_count,progression):

bezler quadratic curve equation %

expansion_ratio = (Re/Rt)**2
Nx = 0.382*Rt*np.cos((theta_N - 90)*np.pi/180)
Ny = 1.382*Rt + 0.382*Rt*np.sin((theta_N - 90)*np.pi/180)
Ex = Ln_ratio*( (np.sqrt(expansion_ratio) -1)*Rt + 1.5*Rt*(1/np.cos(15*np.pi/180) -1))/np.tan(15*np.pi/180) # Ln
Ey = Re

m1 = np.tan(theta_N*np.pi/180)
m2 = np.tan(theta_e*np.pi/180)

C1 = Ny - m1*Nx
C2 = Ey - m2*Ex

Qx = (C2 - C1)/(m1-m2)
Qy = (m1*C2 -m2*C1)/(m1-m2)

x = np.zeros(step_count+1)
y = np.zeros(step_count+1)
i=0

for t in getPoints(0,1,step_count,progression):#np.linspace(0,1,step_count):
    x[i] = (1-t)**2*Nx + 2*(1-t)*t*Qx + t**2*Ex 
    y[i] = (1-t)**2*Ny + 2*(1-t)*t*Qy + t**2*Ey 
    i=i+1

return [x,y]

Rocket Dimensions

Rc = 71 # [mm] Rt = 35.5 # [mm] Re = 275 # [mm]

Lc = 100 # [mm] Ln_ratio = 0.8; # Ratio of conical nozzle

theta_t_i = 30 # [deg] theta_N = 30 # [deg] theta_e = 7 # [deg]

step_count_c_e = 40 #20; step_count_t_i = 60 step_count_t_e = 60 #40 step_count_r = 100 #80 progression_i=1/1.1 progression_bell=1.01 progression_r=1/1.05

throat exit section

[x_t_e,y_t_e]=throat_exit(Rt,theta_N,step_count_t_e)

throat inlet section

deltaN_t_i=x_t_e[1]-x_t_e[0] throat_i_l = 1.5Rtnp.sin(theta_t_inp.pi/180) delta0_t_i=deltaN_t_iprogression_i-throat_i_l *(progression_i-1) step_count_t_i=int(np.ceil(np.log(deltaN_t_i/delta0_t_i)/np.log(progression_i))) [x_t_i,y_t_i]=throat_inlet(Rt,theta_t_i,step_count_t_i,progression_i)

chamber exit section

[x_c_e,y_c_e]=chamber_exit_circ(Rc,Rt,Lc,theta_t_i,step_count_c_e) delta0_c_e=x_c_e[1]-x_c_e[0]

chamber exit to throat inlet line

deltaN_ce_ti=x_t_i[1]-x_t_i[0] L_c_e_ti=np.abs(np.tan((90 + theta_t_i)np.pi/180)(y_c_e[-1]-y_t_i[0])) delta0_ce_ti=x_c_e[-1]-x_c_e[-2] progression_ce_ti=(L_c_e_ti-delta0_ce_ti)/(L_c_e_ti-deltaN_ce_ti) step_count_ce_ti=int(np.ceil(np.log(deltaN_ce_ti/delta0_ce_ti)/np.log(progression_ce_ti))) y_ce_ti = getPoints(y_c_e[-1],y_t_i[0],step_count_ce_ti,progression_ce_ti) x_ce_ti = np.tan((90 + theta_t_i)np.pi/180)y_ce_ti

chamber section

step_count_c=int(np.ceil(np.log(Lc(1/progression_i-1)/delta0_c_e+1)/np.log(1/progression_i)))-1 [x_c,y_c]=chamber(Rc,step_count_c,Lc,progression_i) delta0_c=LcgetS(step_count_c,progression_i)

chamber exit section (translate)

x_c_e = x_c_e + x_c[-1] - x_c_e[0]

translate chamber exit to throat inlet, throat inlet and throat exit section sections

x_ce_ti = x_ce_ti + x_c_e[-1] - x_ce_ti[0] x_t_i = x_t_i + x_ce_ti[-1] - x_t_i[0] x_t_e = x_t_e + x_t_i[-1] - x_t_e[0]

bell section

expansion_ratio = (Re/Rt)*2 delta0_t_e=x_t_e[-1]-x_t_e[-2] Ln=Ln_ratio( (np.sqrt(expansion_ratio) -1)Rt + 1.5Rt(1/np.cos(15np.pi/180) -1))/np.tan(15np.pi/180) step_count_nozzle=int(np.ceil(np.log(Ln(progression_bell-1)/delta0_t_e+1)/np.log(progression_bell)))-1 [x_bell,y_bell]=bell_curve(Rt,Re,theta_N,theta_e,Ln_ratio,step_count_nozzle,progression_bell)

translate bell section

x_bell = x_bell + x_t_e[-1] - x_bell[0]

gmsh.initialize() gmsh.model.add("RaoNozzle")

factory=gmsh.model.occ

lc = 1

pos=np.loadtxt('/Users/jsalazar/matlab/RAO-Bell-Nozzle/nozzleContour.csv',delimiter=',', skiprows=1)

x=np.concatenate((x_c,x_c_e[1:],x_ce_ti[1:],x_t_i[1:],x_t_e[1:],x_bell[1:])) y=np.concatenate((y_c,y_c_e[1:],y_ce_ti[1:],y_t_i[1:],y_t_e[1:],y_bell[1:]))

pos=np.concatenate((x,y)) pos=pos.reshape((len(x),2),order='F') npts=len(pos) erase=[]

for i in range(0,npts-1): if np.linalg.norm(pos[i]-pos[i+1])==0: erase.append(i)

pos=np.delete(pos,erase,0) npts=len(pos)

np.savetxt('nozzleContour.csv', pos/1e3,delimiter=',', header='x [m] y [m]')

pAxis=[[]]npts pBell=[[]]npts

lAxis=[[]](npts-1) lBell=[[]](npts-1) lRadial=[[]]*2

for i in range(0,npts): pBell[i]=factory.addPoint(pos[i,0],pos[i,1],0,lc)
pAxis[i]=factory.addPoint(pos[i,0],0,0,lc)

for i in range(0,npts-1): lBell[i]=factory.addLine(pBell[i],pBell[i+1]) lAxis[i]=factory.addLine(pAxis[i],pAxis[i+1])

j=0

for i in [0,npts-1]: print(i)
lRadial[j]=factory.addLine(pAxis[i],pBell[i]) factory.synchronize() gmsh.model.mesh.setTransfiniteCurve(lRadial[j], step_count_r, meshType="Progression", coef=progression_r) j=j+1

lAxis.reverse()
loop=lBell+neg([lRadial[1]])+neg(lAxis)+[lRadial[0]]

cLoop=[] cLoop.append(factory.addCurveLoop(loop)) surface=[] surface.append(factory.addPlaneSurface(cLoop)) factory.synchronize() gmsh.model.mesh.setRecombine(2, surface[0]) factory.synchronize() gmsh.model.mesh.setTransfiniteSurface(surface[0],cornerTags=[pAxis[0],pAxis[-1],pBell[0],pBell[-1]]) factory.synchronize() v=factory.extrude([(2,surface[0])],0,0,1,numElements=[1],recombine=True) factory.synchronize() gmsh.model.addPhysicalGroup(2,[surface[0],v[0][1]], name="frontAndBack") gmsh.model.addPhysicalGroup(2,[v[-1][1]], name="inlet") gmsh.model.addPhysicalGroup(2,[x[1] for x in v[2:npts+1]], name="wall") gmsh.model.addPhysicalGroup(2,[v[npts+1][1]], name="outlet") gmsh.model.addPhysicalGroup(2,[x[1] for x in v[npts+2:2*npts+1]], name="axis") gmsh.model.addPhysicalGroup(3,[v[1][1]], name="internal")

factory.synchronize() gmsh.model.mesh.generate(3) gmsh.write("raoNozzle.msh2")

if '-nopopup' not in sys.argv: gmsh.fltk.run()

gmsh.finalize() `

isuruf commented 2 months ago

This is not a support forum for gmsh. Please use the proper channel.