List length problem

Hello everyone, 
I have a problem with my code, when i execute it, it displays the following error message : Cell In[92], line 150
147 wind.value('system_capacity', farm_capacity)
149 # Run simulation
--> 150 wind.execute()
152 # Print results
153 print('Selected wind turbine :', turbine)

Exception: windpower execution error.
fail(wind_turbine_powercurve_powerout, length_equal=wind_turbine_powercurve_windspeeds): 161

However, i've checked the lenght of the mentionned list and these lists has the same lenght :

print(len(wind_speeds))
print(len(turbine_powers))
170
170

I've have tried several options but it stil doesn't working.
Maybe it's just a little error in my code.

Here is my code :

# Import the required modules
import json
import pandas as pd
from pathlib import Path
import PySAM.Windpower as wp
import numpy as np
import seawater as sw
from CoolProp.CoolProp import PropsSI
import scipy.optimize as optimize
from scipy.optimize import fsolve
from scipy.integrate import quad
import matplotlib.pyplot as plt

# --- Electrolyzer Simulation

# Input parameters for the operating conditions of the PEM electrolyzer.
A_electrodes = 0.5  # Electrodes area- m2
D = 0.0001    # Membrane thickness, [m]
N_cells = 310  # Number of cells in series per stack
R_T = 0.167  # Thermal resistance per stack - K/W 
TAU_T = 29  # Thermal time constant per stack - hr
C_T = TAU_T*3600/R_T    #Thermal capacitance
a_1 = 250  # Faraday efficiency - mA/cm
a_2 = 0.97 #Faraday efficiency - 0…1
h_1 = 70  # Convective heat transfer - W/°C
h_2 = 0.035  # Convective heat transfer - W/°C per A
T_ini = 70+273.15 # Initial electrolyzer temperature - °C
P_ely = 30  # Pressure (constant)- [bar]
T_amb = 20+273.15  # Ambient temperature- °C
T_cw_in = 20+273.15  # Cooling water inlet temperature- °C
T_cw_out = 45+273.15  # Cooling water inlet temperature- °C
lambda_a = 14    # Water content at the anode-membrane 
lambda_c = 10    # Water content at the cathode-membrane 
J_ref_a = 170000    # Anode pre-exponential factor, [A/m^2]
J_ref_c = 4600    # Cathode pre-exponential factor, [A/m^2]
E_act_a = 76000    # Activation energy for anode, [J/mol]
E_act_c = 18000    # Activation energy for cathode, [J/mol]
HHV_H2=141804000   # Higher heating value of Hydrogen [J/kg]

# Constants and model parameters
DELT = 3600  # time interval -s
TSTD = 0
P_std = 1.0
TREF = 25.0+273.15  # REFERENCE TEMPERATURE [C]
PREF = 1.0  # REFERENCE PRESSURE [bar]
F = 96485  # FARADAY CONSTANT [C/mol] or [As/mol]
RGAS = 8.3145  # UNIVERSAL GAS CONSTANT [J/K-mol]
n = 2  # NUMBER OF ELECTRONS PER MOLECULE OF H2 [#]
M_H2O = 18.016  # Molecular weight of H2O
RHO_H2O = 1000  # water density
DELTAH0F_H2O = -285800     #enthalpy of formation- J/mol
DELTAH0F_H2 = 0
DELTAH0F_O2 = 0
C_p0_H2O = 75 # molar Heat capacity H2O - JK−1mol−1.
C_p0_H2 = 29  # molar Heat capacity H2 - JK−1mol−1.
C_p0_O2 = 29  # molar Heat capacity O2 - JK−1mol−1.
V_cw = 93 # Cooling water flow rate- Nm3/hr
C_p_H2O = V_cw*RHO_H2O/M_H2O*C_p0_H2O/3600*1000    #Heat capacity of cooling water
c_water = 4200 # [J/kgK]

#solving equations
def lambda_function(x, lambda_a, lambda_c, D, T):
    term1 = 0.5139 * (((lambda_a - lambda_c) / D) * x + lambda_c) - 0.326
    term2 = np.exp(1268 * ((1 / 303) - (1 / T)))
    return 1 / (term1 * term2)

I_L = 1.5e4                                                                        # The current limiting term

def PEM(P):                                                                        # Required power [Watts]
    T = T_ini
    I = P/N_cells/2
    DELTAH = ((T-TREF)*(C_p0_H2 + 0.5*C_p0_O2 - C_p0_H2O) - DELTAH0F_H2O)
    U_rev = 1.229-8.5*10**-4*(T-298)
    U_tn = DELTAH/(n*F)
    Fun_R_ohm = lambda x: 1/((0.5139*(((lambda_a-lambda_c)/D)*x+lambda_c)-0.326)*np.exp(1268*((1/303)-(1/T_ini))))
    R_ohm, _ = quad(Fun_R_ohm, 0, D)
    U_ohm=I/A_electrodes*R_ohm
    J_0_a = J_ref_a*np.exp(-E_act_a/(RGAS*T))
    U_act_a=(RGAS*T/F)*np.arcsinh(I/A_electrodes/(2*J_0_a))
    J_0_c = J_ref_c*np.exp(-E_act_c/(RGAS*T))
    U_act_c=(RGAS*T/F)*np.arcsinh(I/A_electrodes/(2*J_0_c))
    U_cell= U_rev + U_act_a + U_act_c + U_ohm
    Q_dot_gen= N_cells * I * (U_cell - U_tn) 
    Q_dot_loss= (T - T_amb) / R_T
    Q_dot_cw= Q_dot_gen - Q_dot_loss
    m_dot_cw = Q_dot_cw/(4200*20)
    ETA_e = U_tn/ U_cell                                       # Energy efficiency
    P_in=U_cell*I*N_cells                       
    m_dot_H2=1e-3*P_in*ETA_e/HHV_H2*3600                       # Produced Hydrogen mass flow rate [ton/hr]
    m_dot_O2 = 8 * m_dot_H2                                    # Produced Oxygen mass flow rate  [ton/hr]
    m_dot_water = m_dot_H2 + m_dot_O2
    results = [m_dot_H2,m_dot_O2,m_dot_water,ETA_e]
    print('Produced hydrogen mass flow rate : ', round(m_dot_H2, 6), 'ton/hr')
    print('Produced oxygen mass flow rate : ', round(m_dot_O2,6), 'ton/hr')
    print('Water mass flow rater required : ', round(m_dot_water, 6), 'ton/hr')
    print('Electrolyzer efficiency : ', round(ETA_e, 3))

    return (results)

results = PEM (17000000)
m_dot_water = results[2]
m_dot_H2 = results[0]

# --- Wind Turbine Simulation ---

# Create a new windpower model
wind = wp.new()

# Set path for the working directory
working_dir = "C:/Use0-07-2023"

# Set model inputs from JSON generated by SAM
with open(working_dir + '/untitled_windpower.json', 'r') as file:
    sam_data = json.load(file)
    # Loop through each key-value pair
    for k, v in sam_data.items():
        if k != 'number_inputs':
            wind.value(k, v)

# Replace some inputs as needed for this analysis

# Load a wind resource file
wind.value('wind_resource_filename', working_dir + '_pressure_modiv')

# Wind shear is location-specific to adjust wind speed data to turbine hub height
wind.value('wind_resource_shear', 0.14)

# Load wind turbine parameters from JSON
with open(working_dir + '/WMW.json', 'r') as file:
    turbine_data = json.load(file)

# Set up wind farm for one turbine
wind.value('wind_farm_xCoordinates', [0])
wind.value('wind_farm_yCoordinates', [0])

# Loop through turbine list for single turbine
for turbine in turbine_data:
    # Set wind turbine parameters
    wind.value('wind_turbine_rotor_diameter', turbine)
    wind.value('wind_turbine_powercurve_windspeeds', turbine)
    wind.value('wind_turbine_powercurve_powerout', turbine)
    wind.value('wind_turbine_hub_ht', turbine)

    # Set wind farm capacity
    number_of_turbines = len(wind.value('wind_farm_xCoordinates'))
    farm_capacity = number_of_turbines * turbine
    wind.value('system_capacity', farm_capacity)

    # Run simulation
    wind.execute()

    # Print results
    print('Selected wind turbine :', turbine)
    print('Annual produced energy (kWh) = ', round(wind.Outputs.annual_energy, 3))
    print('Capacity factor = ', round(wind.Outputs.capacity_factor, 3))
    hourly_energy = wind.Outputs.gen
    hf = pd.DataFrame(hourly_energy)

        
    # --- RO Desalination Plant Simulation ---

    # Operating Conditions
    T_feed = 25  # Feed water temperature in C
    x_feed = 4  # Feed water salinity in PSU
    m_d = m_dot_water /3.6     # Required desalinated water in kg/s           

    
    # Reverse Osmosis properties
    RR = 0.4  # Membrane recovery ratio
    n = 100  # Total number of membranes
    A_mem = 37  # Active membrane area in m^2
    R = 0.9975  # Membrane rejection coefficient
    Km = 8.03E11  # Membrane permeability resistance in Pa.s/m
    Ds = 1.45E7  # Solute diffusivity in m^2/s
    d = 0.71E-3  # Feed spacer thickness or feed channel thickness in m
    Nch = 59  # Number of feed channels
    Np = 100  # Number of pressure vessels
    Lw = 0.243  # Membrane width in m
    Eta_pump = 0.85  # RO pump efficiency
    Eta_hydroturbine = 0.85  # Hydroturbine efficiency

    # Calculations
    m_feed = m_d / RR  # Feed flow rate in kg/s
    m_rej = m_feed - m_d  # Rejected flow rate in kg/s
    
    P_feed = 101325  # Assuming atmospheric pressure in Pa
    rho_feed = sw.dens0(T_feed, x_feed)  # Seawater density in kg/m^3

    # Assuming that the effect of salinity on seawater viscosity is small,
    # we calculate the viscosity of pure water at the feed temperature and pressure.
    mu = PropsSI('V', 'P', P_feed, 'T', T_feed + 273.15, 'Water')  # Dynamic viscosity of water in Pa.s

    Jw = m_feed / (n * A_mem * rho_feed)  # Volumetric permeate flow rate

    Sc = mu / (rho_feed * Ds)  # Schmidt number
    Re = m_feed / (Nch * Lw * mu * Np)  # Reynolds number
    Kmass = 0.04 * Re**0.75 * Sc**0.33 * Ds / d  # Mass transfer coefficient
    Cw = (np.exp(Jw / Kmass) * x_feed) / (
            np.exp(Jw / Kmass) * (1 - R) + R)  # Membrane wall concentration

    DELTAphi = 805.1E5 * Cw * R  # Transmembrane osmotic pressure
    DELTAP = Jw * Km + DELTAphi  # Transmembrane pressure

    # Alternative correlations for osmotic pressure calculations
    i = 2  # No. of ions, 2 for NaCl
    q = 0.92  # Dissociation and random pairing coefficient
    Mm_nacl = 58.44 / 1000  # Molar mass of NaCl in kg/mol
    DELTAphi_2 = i * q * (x_feed / Mm_nacl) * R * (T_feed + 273.15)  # Transmembrane osmotic pressure - correlation2

    # Power required by the RO pump
    W_pump = DELTAP * m_feed / (Eta_pump * rho_feed)
    print('Wpump = ', (W_pump/1000), 'kW')
    # Recovery power by energy recovery device (ERD)
    W_hydropump = DELTAP * m_rej * Eta_pump / rho_feed
    print('Whydropump = ', (W_hydropump/1000), 'kW')
    # Total RO required power
    W_RO = (W_pump - W_hydropump) / 1000

    # Define simulation duration
    t_hours = 8760

    Annual_energy_required = W_RO * t_hours

    print(f"Total RO required power: {round(W_RO, 3)} kW")
    print(f"Annual RO energy required: {round(Annual_energy_required,3)} kWh")
    print('

---

')
    

hourly_energy = wind.Outputs.gen
hourly_energy_series = pd.Series(hourly_energy)

hf = pd.DataFrame(hourly_energy_series)
#print('Hourly energy (kWh) = ', round(hf, 3))

total_energy_produced = hourly_energy_series.sum()

produced_freshwater_rate = m_dot_water * t_hours / total_energy_produced

total_produced_freshwater = produced_freshwater_rate * total_energy_produced

produced_hydrogen_rate = m_dot_H2 * t_hours / total_energy_produced
total_produced_hydrogen = produced_hydrogen_rate * total_energy_produced

#print('Hourly produced freshwater rate (ton/kWh) = ', round(produced_freshwater_rate,6))
#print('Total produced freshwater (ton) = ', total_produced_freshwater)
 
cumulative_energy = hourly_energy_series.cumsum()
print('Hourly produced energy (kWh) = ', cumulative_energy)

cumulative_freshwater = cumulative_energy * produced_freshwater_rate
print('Hourly produced freshwater (ton) = ', cumulative_freshwater)

cumulative_hydrogen = cumulative_energy * produced_hydrogen_rate
print('Hourly produced hydrogen (ton) = ', cumulative_hydrogen)
plt.plot(cumulative_freshwater, cumulative_energy)
plt.xlabel('Produced Freshwater (ton)')
plt.ylabel('Energy Delivered by Wind Turbine (kWh)')
plt.title('Energy produced by the wind turbine as a function of the freshwated produced')
plt.grid(True)
plt.xlim(0, cumulative_freshwater.max() * 1.1)
plt.ylim(0, cumulative_energy.max() * 1.1)

plt.figure()
plt.plot(cumulative_hydrogen, cumulative_energy, color = 'orange')
plt.xlabel('Produced hydrogen (ton)')
plt.ylabel('Energy Delivered by Wind Turbine (kWh)')
plt.title('Energy produced by the wind turbine as a function of the hydrogen produced')
plt.grid(True)
plt.xlim(0, cumulative_hydrogen.max() * 1.1)
plt.ylim(0, cumulative_energy.max() * 1.1)

plt.show()




I will be very glad if someone could help me.

Good day to you all


 

Hi Martineau,

If you would like me to help troubleshoot your code, please provide a short version of the code with all supporting files, including only the essential part of the code that causes the problem.

Best regards,
Paul.