Renewable energy engineer: MQL recipes

Recipes

Yearly yield of a rooftop PV array
Annual energy of a wind turbine with a Weibull wind
Battery bank of an off-grid cabin
Levelized cost of electricity of a solar farm
Daily energy and ramp of a building load curve
Sun angles and plane of array irradiance on a tilted PV panel
ABCD constants and Ferranti effect of a 400 kV line
Size a counterflow heat exchanger with the effectiveness-NTU method
DC power flow of a 4-bus grid with an LU factorisation
Power factor correction of a workshop

executed on MentDB

Yearly yield of a rooftop PV array

From the monthly in-plane irradiation and a performance ratio, compute the yearly energy of a 6 kWp array, then the DC power of a hot summer hour with the cell temperature.

Pipeline

01 · irradiationSum of the monthsstat sum→02 · yieldE = P H PRmath calc→03 · thermalNOCT cell temperaturemath calc

Script

admin

#1. 6 kWp rooftop array: monthly in-plane irradiation in kWh/m2, performance ratio 0.80;
-> "[H]" "[55, 78, 120, 150, 172, 180, 192, 170, 135, 95, 58, 46]";
-> "[P]" 6;
-> "[Hy]" (stat sum [H]);
#2. Yearly energy E = P * H / 1 kW/m2 * PR;
-> "[E]" (math calc "[P]*[Hy]*0.80");
#3. Hot summer hour: 950 W/m2, 33 degC ambient, NOCT 45 degC, power coefficient -0.35 %/K;
-> "[Tc]" (math calc "33+(45-20)/800*950");
-> "[Pdc]" (math calc "[P]*950/1000*(1-0.0035*([Tc]-25))");
concat "yearly yield = " (math round [E] 0) " kWh, specific yield = " (math round (math calc "[E]/[P]") 0) " kWh/kWp, best month = " (math round (math calc (concat (stat max [H]) "*[P]*0.80")) 0) " kWh, cell at " (math round [Tc] 1) " degC gives " (math round [Pdc] 2) " kW DC"

Result

mentdb

yearly yield = 6965 kWh, specific yield = 1161 kWh/kWp, best month = 922 kWh, cell at 62.7 degC gives 4.95 kW DC

Go further

Add the inverter clipping on hourly data
Compare the orientations with a table of irradiations

executed on MentDB

Annual energy of a wind turbine with a Weibull wind

Find the rated wind speed of a 3.6 MW turbine, integrate its power curve over a Weibull wind distribution and get the annual energy and the capacity factor.

Pipeline

01 · ratedRated wind speednum solve→02 · integratePower x densitynum integrate→03 · rated partWeibull probabilitiesstat dist_cdf→04 · aepAEP and capacity factormath calc

Script

admin

#1. 3.6 MW turbine, rotor 120 m, Cp 0.42, air 1.20 kg/m3, cut-in 3 and cut-out 25 m/s;
-> "[A]" (math calc "pi*120^2/4");
-> "[Pv]" (concat "0.5*1.20*" [A] "*0.42*v^3/1e6");
#2. Rated wind speed where the power reaches 3.6 MW;
-> "[vr]" (num solve (concat [Pv] "-3.6") "v" 5 20);
#3. Weibull wind k 2.1, c 8.2 m/s: integrate the power curve times the density;
-> "[pdf]" "(2.1/8.2)*(v/8.2)^1.1*exp(-(v/8.2)^2.1)";
-> "[part]" (num integrate (concat [Pv] "*" [pdf]) "v" 3 [vr]);
-> "[rated]" (math calc (concat "3.6*(" (stat dist_cdf "weibull" "[2.1, 8.2]" 25) "-" (stat dist_cdf "weibull" "[2.1, 8.2]" [vr]) ")"));
#4. Annual energy production and capacity factor;
-> "[aep]" (math calc "8760*([part]+[rated])");
concat "rated wind = " (math round [vr] 2) " m/s, mean power = " (math round (math calc "[part]+[rated]") 3) " MW, AEP = " (math round [aep] 0) " MWh, capacity factor = " (math round (math calc "100*[aep]/(3.6*8760)") 1) " %"

Result

mentdb

rated wind = 10.81 m/s, mean power = 1.409 MW, AEP = 12347 MWh, capacity factor = 39.2 %

Go further

Correct the wind to hub height with the power law
Fit the Weibull parameters on measured winds with stat fit_model

executed on MentDB

Battery bank of an off-grid cabin

Size the battery of an off-grid cabin for three days of autonomy with the depth of discharge, the efficiency and the cold derating, in kWh, Ah and modules.

Pipeline

01 · needNominal capacitymath calc→02 · convertkWh to Wh, then Ahunit convert→03 · modulesWhole modulesmath ceil

Script

admin

#1. Off-grid cabin: 6.5 kWh per day, 3 days of autonomy, LFP depth of discharge 80 %, efficiency 92 %, cold derating 95 %;
-> "[need]" (math calc "6.5*3/(0.80*0.92*0.95)");
#2. Capacity of the 48 V bank in Ah;
-> "[Ah]" (math calc (concat (unit convert [need] "kWh" "Wh") "/48"));
#3. Number of 51.2 V 100 Ah modules of 5.12 kWh, and C-rate for a 4 kW peak load;
-> "[n]" (math ceil (math calc "[need]/5.12"));
-> "[C]" (math calc "4/([n]*5.12)");
concat "nominal capacity = " (math round [need] 2) " kWh = " (math round [Ah] 0) " Ah at 48 V, " [n] " modules = " (math calc "[n]*5.12") " kWh, peak C-rate = " (math round [C] 3) " C"

Result

mentdb

nominal capacity = 27.89 kWh = 581 Ah at 48 V, 6 modules = 30.72 kWh, peak C-rate = 0.13 C

Go further

Size the PV array that recharges it in one day
Check the cycle life with the yearly throughput

executed on MentDB

Levelized cost of electricity of a solar farm

Discount the costs and the degrading production of a 10 MWp solar farm over 25 years to get its LCOE, and check it with the capital recovery factor.

Pipeline

01 · costsDiscounted opexnum sum→02 · energyDiscounted energynum sum→03 · lcoeCost / energymath calc

Script

admin

#1. 10 MWp solar farm: capex 7.5 Meur, opex 120 keur per year, 15 GWh first year, degradation 0.5 % per year, WACC 6 %, 25 years;
-> "[r]" 0.06;
#2. Discounted costs and discounted energy;
-> "[cost]" (math calc (concat "7.5e6+" (num sum "120000/(1.06)^t" "t" 1 25)));
-> "[energy]" (num sum "15000*(1-0.005)^(t-1)/(1.06)^t" "t" 1 25);
#3. Levelized cost of electricity in eur/MWh;
-> "[lcoe]" (math calc "[cost]/[energy]");
#4. Annuity check with the capital recovery factor, no degradation;
-> "[crf]" (math calc "[r]*(1+[r])^25/((1+[r])^25-1)");
-> "[simple]" (math calc "(7.5e6*[crf]+120000)/15000");
concat "LCOE = " (math round [lcoe] 2) " eur/MWh = " (math round (math calc "[lcoe]/10") 2) " c/kWh, discounted energy = " (math round [energy] 0) " MWh, CRF = " (math round [crf] 5) ", annuity LCOE without degradation = " (math round [simple] 2) " eur/MWh"

Result

mentdb

LCOE = 49.28 eur/MWh = 4.93 c/kWh, discounted energy = 183333 MWh, CRF = 0.07823, annuity LCOE without degradation = 47.11 eur/MWh

Go further

Find the WACC that gives a 45 eur/MWh tariff with num solve
Add the inverter replacement in year 12

executed on MentDB

Daily energy and ramp of a building load curve

A building load is measured every 2 hours. The script integrates a monotone interpolation of the curve to get the daily energy, finds the fastest ramp on a 6-minute grid to size a battery inverter, and draws the load curve through the measures.

Pipeline

01 · energyDaily energy in kWhnum interpolate_integral→02 · rampFastest rampnum interpolate_derivative→03 · drawLoad curveplot interpolate

Script

admin

#1. Load of an office building in kW, measured every 2 h;
-> "[h]" "[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]";
-> "[P]" "[42, 40, 41, 55, 160, 185, 170, 182, 178, 120, 70, 48, 42]";
#2. Daily energy: exact integral of the monotone interpolant, kWh;
-> "[E]" (num interpolate_integral [h] [P] 0 24 "pchip");
#3. Fastest ramp in kW/h on a 6 min grid, to size the battery inverter;
-> "[ramp]" (stat max (num interpolate_derivative [h] [P] (matrix col (num table "t" "t" 0 24 0.1) 1) "pchip"));
#4. Load curve through the measures;
plot interpolate [h] [P] "pchip"

Result

mentdb

The chart is shown in the Plot tab of the editor.

Go further

Return [E] and [ramp] with concat to read the numbers (about 2.58 MWh and 71 kW/h).
Compare the modes linear, cubic and pchip: cubic overshoots near the morning step.

executed on MentDB

Sun angles and plane of array irradiance on a tilted PV panel

For a PV plant tilted 30 deg south on the summer solstice at 10 h solar time, the script gets the declination from the day of the year, the sun elevation and azimuth with atan2, the angle of incidence on the panel and the plane of array irradiance with the isotropic sky model.

Pipeline

01 · dayDay of the yeardate day_of_year→02 · sunElevationmath calc→03 · azimuthAzimuth from the southmath atan2→04 · incidenceAngle of incidencemath acos

Script

admin

#1. PV plant at 43.6 N, panels tilted 30 deg facing south, 10 h solar time on 2026-06-21;
-> "[lat]" (math deg_to_rad 43.6); -> "[beta]" (math deg_to_rad 30);
-> "[N]" (date day_of_year "2026-06-21");
-> "[dec]" (math calc "23.45*pi/180*sin(2*pi*(284+[N])/365)");
-> "[w]" (math deg_to_rad (math calc "15*(10-12)"));
#2. Sun elevation and azimuth from the south, west positive, with atan2;
-> "[alt]" (math calc "asin(sin([lat])*sin([dec])+cos([lat])*cos([dec])*cos([w]))");
-> "[az]" (math atan2 (math calc "sin([w])") (math calc "cos([w])*sin([lat])-tan([dec])*cos([lat])"));
#3. Angle of incidence on the tilted panel;
-> "[cos]" (math calc "sin([alt])*cos([beta])+cos([alt])*sin([beta])*cos([az])");
#4. Plane of array irradiance: DNI 820 W/m2, DHI 120 W/m2, ground albedo 0.2;
-> "[ghi]" (math calc "820*sin([alt])+120");
-> "[poa]" (math calc "820*[cos]+120*(1+cos([beta]))/2+0.2*[ghi]*(1-cos([beta]))/2");
concat "day " [N] ", declination = " (math round (math rad_to_deg [dec]) 2) " deg, sun elevation = " (math round (math rad_to_deg [alt]) 2) " deg, azimuth = " (math round (math rad_to_deg [az]) 2) " deg from south, incidence = " (math round (math rad_to_deg (math acos [cos])) 2) " deg, GHI = " (math round [ghi] 0) " W/m2, POA = " (math round [poa] 0) " W/m2"

Result

mentdb

day 172, declination = 23.45 deg, sun elevation = 58.19 deg, azimuth = -60.48 deg from south, incidence = 30.03 deg, GHI = 817 W/m2, POA = 833 W/m2

Go further

Loop on the hours of the day and integrate the POA with num integrate to get the daily irradiation.
Find the tilt that maximises the irradiance of a given day with num maximize.

executed on MentDB

ABCD constants and Ferranti effect of a 400 kV line

For a 300 km, 400 kV overhead line, the script computes in complex numbers the propagation constant and the characteristic impedance, the exact ABCD constants with the complex hyperbolic cosine and sine, the voltage rise at no load (Ferranti effect), the surge impedance loading and the sending end voltage at 500 MW.

Pipeline

01 · gammaPropagation constantcomplex sqrt→02 · aA = cosh gamma lcomplex cosh→03 · bB = Zc sinh gamma lcomplex sinh→04 · loadVs = A Vr + B Ircomplex from_polar→05 · readModulus and anglecomplex abs

Script

admin

#1. 400 kV line of 300 km, per km: z = 0.03 + 0.35i ohm, y = 4.4e-6i S;
-> "[z]" "0.03+0.35i"; -> "[y]" "0.0000044i"; -> "[l]" 300;
#2. Propagation constant and characteristic impedance;
-> "[gl]" (complex mul (complex sqrt (complex mul [z] [y])) [l]);
-> "[Zc]" (complex sqrt (complex div [z] [y]));
#3. ABCD constants of the long line: A = cosh gamma l, B = Zc sinh gamma l;
-> "[A]" (complex cosh [gl]);
-> "[B]" (complex mul [Zc] (complex sinh [gl]));
#4. Ferranti effect at no load, surge impedance loading and the sending end voltage at full load;
-> "[Vr]" (math calc "400e3/sqrt(3)");
-> "[Ir]" (complex from_polar (math calc "500e6/(sqrt(3)*400e3*0.95)") (math calc "-acos(0.95)"));
-> "[Vs]" (complex add (complex mul [A] [Vr]) (complex mul [B] [Ir]));
concat "gamma l = " [gl] ", Zc = " [Zc] " ohm, A = " [A] [_n_] "no load rise = " (math round (math calc (concat "100*(1/" (complex abs [A]) "-1)")) 2) " %, SIL = " (math round (math calc (concat "400^2/" (complex abs [Zc]))) 0) " MW, sending voltage at 500 MW pf 0.95 = " (math round (math calc (concat (complex abs [Vs]) "*sqrt(3)/1000")) 1) " kV, load angle = " (math round (math rad_to_deg (complex arg [Vs])) 2) " deg"

Result

mentdb

gamma l = 0.015940688761504986+0.3726313265926406i, Zc = 282.2964595398792-12.076279364776505i ohm, A = 0.9314909267241552+0.005803731455354173i
no load rise = 7.35 %, SIL = 566 MW, sending voltage at 500 MW pf 0.95 = 444.1 kV, load angle = 16.63 deg

Go further

Size the shunt reactors that cancel the Ferranti rise at no load.
Compare with the nominal pi model: A = 1 + ZY/2.

executed on MentDB

Size a counterflow heat exchanger with the effectiveness-NTU method

A heat recovery loop needs an effectiveness of 0.80. The script solves the counterflow effectiveness-NTU relation for the NTU, gets the exchange area, the duty and the outlet temperatures, and checks the duty with the log mean temperature difference.

Pipeline

01 · ratioCapacity ratiomath min→02 · ntuNTU for epsilon = 0.8num solve→03 · dutyDuty and outletsmath calc→04 · checkLMTDmath log

Script

admin

#1. Counterflow heat exchanger of a heat recovery loop: hot side 2.5 kW/K, cold side 1.8 kW/K, U = 850 W/m2/K;
-> "[Ch]" 2500; -> "[Cc]" 1800; -> "[U]" 850;
-> "[Cmin]" (math min [Ch] [Cc]); -> "[Cr]" (math calc (concat [Cmin] "/" (math max [Ch] [Cc])));
#2. NTU needed for an effectiveness of 0.80, epsilon = 1 - e^-NTU 1-Cr / 1 - Cr e^-NTU 1-Cr;
-> "[NTU]" (num solve (concat "(1-exp(-x*(1-" [Cr] ")))/(1-" [Cr] "*exp(-x*(1-" [Cr] ")))-0.80") "x" 0.1 20);
-> "[A]" (math calc "[NTU]*[Cmin]/[U]");
#3. Duty and outlet temperatures for inlets at 90 C and 15 C;
-> "[q]" (math calc "0.80*[Cmin]*(90-15)");
-> "[Tho]" (math calc "90-[q]/[Ch]"); -> "[Tco]" (math calc "15+[q]/[Cc]");
#4. Check with the log mean temperature difference;
-> "[lmtd]" (math calc (concat "((90-[Tco])-([Tho]-15))/" (math log (math calc "(90-[Tco])/([Tho]-15)"))));
concat "Cr = " (math round [Cr] 3) ", NTU = " (math round [NTU] 4) ", area = " (math round [A] 2) " m2, duty = " (math round (math calc "[q]/1000") 1) " kW, hot out = " (math round [Tho] 2) " C, cold out = " (math round [Tco] 2) " C, LMTD = " (math round [lmtd] 3) " K, U A LMTD = " (math round (math calc "[U]*[A]*[lmtd]/1000") 1) " kW"

Result

mentdb

Cr = 0.72, NTU = 2.6836, area = 5.68 m2, duty = 108 kW, hot out = 46.8 C, cold out = 75 C, LMTD = 22.358 K, U A LMTD = 108 kW

Go further

Plot the effectiveness against the NTU for counterflow and parallel flow with plot function.
Add the fouling resistances to U and see the area grow.

executed on MentDB

DC power flow of a 4-bus grid with an LU factorisation

The DC load flow solves B theta = P on the buses without the slack. The script factorises the susceptance matrix with LU and partial pivoting, does the forward and back substitutions, gets the bus angles, the determinant from the diagonal of U and the active power flows of the five lines.

Pipeline

01 · factorP B = L Umatrix lu→02 · forwardL y = P bmatrix solve→03 · backU theta = ymatrix get→04 · flowsLine flows in MWmath calc

Script

admin

#1. DC power flow of a 4-bus grid, bus 1 is the slack, reactances in pu: 1-2 0.1, 1-3 0.2, 2-3 0.25, 2-4 0.2, 3-4 0.1;
#2. Reduced susceptance matrix of buses 2, 3, 4 and injections in pu on 100 MVA;
-> "[B]" "[[19, -4, -5], [-4, 19, -10], [-5, -10, 15]]";
-> "[P]" "[0.5, -1.2, -0.8]";
#3. LU factorisation with partial pivoting P B = L U, then forward and back substitution;
json load "rx_p4e_lu" (matrix lu [B]);
-> "[y]" (matrix solve (json select "rx_p4e_lu" "/L") (matrix mul (json select "rx_p4e_lu" "/P") [P]));
-> "[th]" (matrix solve (json select "rx_p4e_lu" "/U") [y]);
-> "[detU]" (math calc (concat (matrix get (json select "rx_p4e_lu" "/U") 1 1) "*" (matrix get (json select "rx_p4e_lu" "/U") 2 2) "*" (matrix get (json select "rx_p4e_lu" "/U") 3 3)));
json unload "rx_p4e_lu";
#4. Angles and line flows P = theta i - theta j / x in MW;
json load "rx_p4e_th" (concat "{\"t\": " [th] "}");
-> "[t2]" (json select "rx_p4e_th" "/t/[0]"); -> "[t3]" (json select "rx_p4e_th" "/t/[1]"); -> "[t4]" (json select "rx_p4e_th" "/t/[2]");
json unload "rx_p4e_th";
concat "angles = " (math round (math rad_to_deg [t2]) 3) ", " (math round (math rad_to_deg [t3]) 3) ", " (math round (math rad_to_deg [t4]) 3) " deg, det B = " (math round [detU] 4) [_n_] "flows MW: 1-2 " (math round (math calc "100*(0-[t2])/0.1") 1) ", 1-3 " (math round (math calc "100*(0-[t3])/0.2") 1) ", 2-3 " (math round (math calc "100*([t2]-[t3])/0.25") 1) ", 2-4 " (math round (math calc "100*([t2]-[t4])/0.2") 1) ", 3-4 " (math round (math calc "100*([t3]-[t4])/0.1") 1)

Result

mentdb

angles = -3.521, -10.146, -10.994 deg, det B = 2400
flows MW: 1-2 61.5, 1-3 88.5, 2-3 46.3, 2-4 65.2, 3-4 14.8

Go further

Reuse the same LU for several injection scenarios: only the substitutions change.
Get the PTDF matrix from matrix inverse to see how a new generator loads each line.

executed on MentDB

Power factor correction of a workshop

A workshop draws 400 kW at a power factor of 0.78. The script writes the apparent power as a complex number, sizes the capacitor bank for a power factor of 0.95 with the tangents of the angles, gets the capacitance per phase of a delta bank and the reduction of the line current and of the Joule losses.

Pipeline

01 · anglephi = acos pfmath acos→02 · powerS = P + jQcomplex from_polar→03 · bankQc = P (tan phi1 - tan phi2)math tan→04 · afterNew S and pfcomplex sub

Script

admin

#1. Workshop on a 400 V 50 Hz network: 400 kW of motors at power factor 0.78 lagging, target 0.95;
-> "[P]" 400; -> "[pf1]" 0.78; -> "[pf2]" 0.95;
-> "[f1]" (math acos [pf1]); -> "[f2]" (math acos [pf2]);
#2. Apparent power as a complex number S = P + jQ;
-> "[S1]" (complex from_polar (math calc "[P]/[pf1]") [f1]);
#3. Capacitor bank Qc = P tan phi1 - tan phi2;
-> "[Qc]" (math calc (concat "[P]*(" (math tan [f1]) "-" (math tan [f2]) ")"));
-> "[S2]" (complex sub [S1] (concat [Qc] "i"));
#4. Capacitance per phase of a delta bank and line current saving;
-> "[C]" (math calc "[Qc]*1000/(3*2*pi*50*400^2)*1e6");
concat "S before = " [S1] " kVA, Qc = " (math round [Qc] 1) " kvar, S after = " (math round (complex abs [S2]) 1) " kVA at pf " (math round (math cos (complex arg [S2])) 4) ", C per phase = " (math round [C] 1) " microF, line current " (math round (math calc (concat (complex abs [S1]) "/(sqrt(3)*0.4)")) 0) " A to " (math round (math calc (concat (complex abs [S2]) "/(sqrt(3)*0.4)")) 0) " A, Joule losses cut by " (math round (math calc (concat "100*(1-(" (complex abs [S2]) "/" (complex abs [S1]) ")^2)")) 1) " %"

Result

mentdb

S before = 400+320.9125712238362i kVA, Qc = 189.4 kvar, S after = 421.1 kVA at pf 0.95, C per phase = 1256.3 microF, line current 740 A to 608 A, Joule losses cut by 32.6 %

Go further

Choose standard steps of 25 kvar with math ceil and check that the pf never goes leading at light load.
Compute the payback from the reactive energy tariff.