Do you think the following code is good? and what changes can I make to make it better. (It's for my tipe, cpge....) if anyone has any idea of another simulation that maybe didn't cross my mind can you please propose it? thank you

import numpy as np
import matplotlib.pyplot as plt

# Define the parameters of the problem
num_panels = 50  # Number of panels
alpha = 5.0  # Angle of attack in degrees
chord_length = 1.0  # Chord length of the airfoil

# Convert angle of attack to radians
alpha = np.radians(alpha)

# Define the airfoil geometry (symmetric airfoil for simplicity)
theta = np.linspace(0, 2*np.pi, num_panels+1)
x = (1 - np.cos(theta)) / 2 * chord_length
y = np.zeros_like(x)

# Compute the panel angles
beta = np.arctan2(y[1:] - y[:-1], x[1:] - x[:-1])

# Compute the length of each panel
S = np.sqrt((x[1:] - x[:-1])**2 + (y[1:] - y[:-1])**2)

# Compute the influence matrix
A = np.zeros((num_panels, num_panels))
for i in range(num_panels):
    for j in range(num_panels):
        if i == j:
            A[i, j] = np.pi
        else:
            A[i, j] = 2 * np.arctan2(y[j+1] - y[j] - (x[j+1] - x[j]) * np.tan(beta[i]), (x[j+1] - x[j]) + (y[j+1] - y[j]) * np.tan(beta[i]))

# Compute the right-hand side of the system
rhs = -2 * np.pi * np.cos(beta - alpha)

# Solve the system for the panel strengths
gamma = np.linalg.solve(A, rhs)

# Compute the tangential velocity on each panel
Vt = np.zeros(num_panels)
for i in range(num_panels):
    for j in range(num_panels):
        if i != j:
            Vt[i] += gamma[j] / (2 * np.pi) * (np.sin(beta[i] - beta[j]) + (x[j+1] - x[j]) / (y[j+1] - y[j]) * (np.cos(beta[i] - beta[j]) - 1))

# Add the free stream velocity
Vt += np.cos(beta - alpha)

# Compute the pressure coefficient on each panel
Cp = 1 - Vt**2

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(x[:-1], Cp, label='Pressure Coefficient')
plt.xlabel('x')
plt.ylabel('Cp')
plt.title('Pressure Coefficient Distribution over the Airfoil')
plt.grid()
plt.legend()
plt.gca().invert_yaxis()
plt.show()

# Plot the airfoil shape
plt.figure(figsize=(10, 6))
plt.plot(x, y, 'k-', label='Airfoil Shape')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Airfoil Shape')
plt.grid()
plt.legend()
plt.axis('equal')
plt.show()

Hi,

the following line of code:

y = np.zeros_like(x)

creates a list array of all zeros. This in turn causes the following line:

 Vt[i] += gamma[j] / (2 * np.pi) * (np.sin(beta[i] - beta[j]) + (x[j+1] - x[j]) / (y[j+1] - y[j]) * (np.cos(beta[i] - beta[j]) - 1))

to generate a divide by zero exception because y[j+1] is always zero and in the denominator of a term in the expression.

So, to rectify this issue, what exactly are the values that you are expecting to populate the numpy array y? If you can answer that, then you can fix the issue. This is the root of the error for generating the divide by zero exception.

Hint:

Comment out the code after the numpy array assignment y. Add the following print statement after the y numpy array assignment to verify that the contents of the array are as expected. Once having verified that the contents of the array are correct, then uncomment out the rest of the code :

print('The value of y is: ', y)

Try eliminating loops. It will be a good exercise and will result in much faster code. numpy is best at vectorized operations.

Some things that might help you.

a = np.zeros(10)
for i in range(10):
    if i = 2:
        a[i] = 1
    else:
        a[i] = 2

Can be replaced with:

a = np.zeros(10)
rng = np.arange(10)
mask = rng == 2
a[mask] = 1
a[~mask] = 2

I am almost certain that broadcasting and masks will do the trick for you.

However, if it is not enough, explore some of numpy’s functions.

Couple of basics that are sometimes helpful in these cases:

np.select
np.where

Finally, if nothing helps take a look at np.einsum.

1 Like