Issue with diverging during runs of own rotor

[Ryde-G]

Hi there,

I am running my own rotor with OpenCOPTER and it is returning nan after around 1300 iterations. I can see that after 1300 iterations that it diverges as it speeds up significantly. Is this a known and/or common problem? If so, what is the cause of the issue so that I can work around it?

The inputs are as follows and I have been using the NACA 0012 airfoil (from the provided airfoil datas)
iterations = 1600
#iterations = 600
wake_history_length = 1*1024

requested_elements = 45

num_rotors = 1
num_blades = 4
R = 7.0104

theta_75 = 4.3*(math.pi/180.0) #pitch of rotor at .75

density = 1.225 #density of air
omega = 34.0 #rad/s
AR = 18.4 #aspect ratio, radius/chord
sos = 340 #speed of sound

V_inf = 20 #The air velocity relative to the body
aoa = 0.3*(math.pi/180.0) #angle of attack
theta_tw_1 = 0*(math.pi/180.0) #twisting angle?

d_psi = 1.0 #rotation angle
dt = d_psi*(math.pi/180.0)/math.fabs(omega)

shed_history_angle = 45.0
shed_history = round(shed_history_angle/d_psi)

d_azimuth = 100*math.pi/num_blades #??

# Root cutout not sure if this is in meters or as a fraction of radius
r_c = 0.2

[Rob-Rau]

Hey thanks for reaching out again!

There are potentially a number of issues causing a divergence. One common issue is having to large or small of a collective value causing either to little thrust generation causing it to go into a vortex ring state, which the method does not like, or too much thrust will also cause problems.

Also, looking at the code you posted, that is a very very large rotor, at a low RPM. I'm not sure how the code will deal with that as I haven't tried something that large before.

[Ryde-G]

Hi! Thanks for getting back so quickly. The thrust value has been around what was expected based on experimental data so hopefully that is not the issue.

We are trying to replicate a more realistic set of parameters for a real helicopter. The example that I am working on has a blade length of 7m and an angular velocity of 34 rad/s. What units is omega in? We assumed rad/s because when we tried putting it in rpm the wake was essentially straight down even in high velocity conditions which I thought must indicate that omega is in rad/s.

[Rob-Rau]

You have the correct units for omega. In most cases, OpenCOPTER is unit agnostic (with the exception of angular rates, which are all rad/s). It just want's consistency in the supplied values.

So, testing this rotor on my machine showed one main problem. You've specified no spanwise twist distribution down the blade. The example.py assumes a linear distribution through the theta_tw_1 variable, but you can give it any distribution you please. When I set theta_tw_1 to -8 (I'm sure other values will work here as well), the simulation will complete without diverging for me. The cause of this is that when the blade was reaching the front of the rotor, it was achieving extremely high loading at the blade tip due to not having any washout. This in turn caused issues in the dynamic inflow model, causing it to diverge.

[Ryde-G]

Thank you so much for your help! Adding the twist made a big difference and the results are a lot closer to the expected value.

By the way, we updated the plotting to be a new 3D plotting routine that uses matplotlib. It is somewhat similar to your vtk output option. Here is the code for it:

def set_axes_equal(ax: plt.Axes):
"""Set 3D plot axes to equal scale.

		Make axes of 3D plot have equal scale so that spheres appear as
		spheres and cubes as cubes.  Required since `ax.axis('equal')`
		and `ax.set_aspect('equal')` don't work on 3D.
		"""
		limits = np.array([
    		ax.get_xlim3d(),
    		ax.get_ylim3d(),
    		ax.get_zlim3d(),
		])
		origin = np.mean(limits, axis=1)
		radius = 0.5 * np.max(np.abs(limits[:, 1] - limits[:, 0]))
		_set_axes_radius(ax, origin, radius)


def _set_axes_radius(ax, origin, radius):
		x, y, z = origin
		ax.set_xlim3d([x - radius, x + radius])
		ax.set_ylim3d([y - radius, y + radius])
		ax.set_zlim3d([z - radius, z + radius])
# Plot
max_dim_1 = -math.inf
min_dim_1 = math.inf
max_dim_2 = -math.inf
min_dim_2 = math.inf
max_dim_3 = -math.inf
min_dim_3 = math.inf
ax = plt.figure().add_subplot(projection='3d')
# Iterating directly on the collections like this is read only
for r_idx, rotor in enumerate(ac_state.rotor_states):
	for b_idx, blade in enumerate(rotor.blade_states):
		x = get_wake_x_component(wake_history.history[0].rotor_wakes[r_idx].tip_vortices[b_idx])
		y = get_wake_y_component(wake_history.history[0].rotor_wakes[r_idx].tip_vortices[b_idx])
		z = get_wake_z_component(wake_history.history[0].rotor_wakes[r_idx].tip_vortices[b_idx])
		max_dim_3 = max(max_dim_2, max(z))
		min_dim_3 = min(min_dim_2, min(z))
		max_dim_2 = max(max_dim_2, max(y))
		min_dim_2 = min(min_dim_2, min(y))
		max_dim_1 = max(max_dim_1, max(x))
		min_dim_1 = min(min_dim_1, min(x))
		ax.plot(x, y, z, linewidth=0.5)
		
		#x_r = math.cos(blade.azimuth)
		#x_b = aircraft.rotors[r_idx].origin[0] + x_r*math.cos(aoa)
		#y_b = aircraft.rotors[r_idx].origin[1]
		#z_b = aircraft.rotors[r_idx].origin[2] - x_r*math.sin(aoa)
		
		#ax.plot([aircraft.rotors[r_idx].origin[0], x_b], [aircraft.rotors[r_idx].origin[1], y_b], [aircraft.rotors[r_idx].origin[2], z_b], "k-", linewidth=1.5)
#end
ax.set_box_aspect([1,1,1])
set_axes_equal(ax)
plt.show()