System of seven ODEs solve using solve_ivp or implement RK4

Question

I'm trying solve a system of coupled ordinary differential equations, formed by 7 ODEs in python, using solve_ivp or either implement a fuction for RK4.

The general physical problem is as follows: Cooling of photovoltaic modules with heat exchanger coupling to the module. In this way, the module generates electrical energy and thermal energy.

I have a polynomial function, G(t) = 9.8385e-13*t^4 - 1.82918e-8*t^3 + 5.991355e-05*t^2 + 2.312059e-1*t + 25, which works for an approximate range of 0 < t < 9000, which represents solar radiation as a function of time of day. This function was obtained through a "polyfit" applied to real data (file upload here. Its a CSV - https://files.fm/u/9y4evkf6c).

This function is used as input for the ODEs, which represent an electrical and a thermal system as a function of time. To solve the electrical model, I created some scripts that solve the diode equation for the photovoltaic module in question, and the output of this script is the photovoltaic power (called in the PPV thermal model) generated as a function of the module temperature and radiation. This script works great and solves part of my problem. My difficulty lies in solving the equations of the thermal model, which receives as input parameters G(t) and PPV.

The equations result in this system:

System of EDOS

Labels:

Tvidro = Tglass = T1

Tcel = Tpv = T2

Ttedlar = T3

Tabs = Tabsorber = T4

Ttubo = Ttube = T5

Tfsai = Tfluid_out = T6

Tiso = Tinsulation = T7

Using method/function for RK4, the complete code is like this (you can go direct to part "#DEFINE MODEL EQUATIONS - ODES)" :

import numpy as np
import matplotlib.pyplot as plt
import csv
from numpy.polynomial.polynomial import polyval

############################################################

with open('directory of data called teste_dados_radiacao',"r") as i:
    rawdata = list(csv.reader(i, delimiter = ";"))

exampledata = np.array(rawdata[1:], dtype=float)

xdata = exampledata[:,0]
ydata = exampledata[:,1]

curve = np.array(np.polyfit(xdata, ydata, 4))

rev_curve = np.array(list(reversed(curve)), dtype=float)
print(rev_curve)

#G_ajustado = polyval(xdata, rev_curve)

""" plt.plot(xdata, ydata, label = "dados experimentais")
plt.plot(xdata, model, label = "model")
plt.legend()
plt.show() """

#############################################################
#CONSTANTS
Tamb = 25 #°C #ambient temperatura
SIGMA = 5.67e-8 #W/m2K4
E_VIDRO = 0.90 #between 0.85 e 0.83 #nasrin2017 0.04
VENTO = 2 #m/s
T_GROUND = Tamb + 2 #°C
T_CEU = 0.00552*Tamb**1.5
Vf = 1 #m/s
Do = 10e-3 #m
Di = 8e-3 #m
NS = 6*10 #number of cells
T_F_ENT = 20 #°C

#INPUTS
Tcel = 25
Tv = 25
Tiso = 30
Av = 1.638*0.982
ALPHA_VIDRO = 0.9
L_VIDRO = 3e-3 #m
RHO_VIDRO = 2500 #kg/m3
M_VIDRO = Av*L_VIDRO*RHO_VIDRO #kg
CP_VIDRO = 500 #j/kgK
K_VIDRO = 2 #W/mK
TAU_VIDRO = 0.95

Pac = 0.85
H_CELL = 0.156 #m
A_CELL = NS*H_CELL**2
ALPHA_CELL = 0.9
L_CEL = 3e-3
RHO_CEL = 2330
M_CEL = A_CELL*L_CEL*RHO_CEL #kg - estimated
CP_CEL = 900 #J/kgK
K_CEL = 140 #W/mK
BETA_T = 0.43/100 # %/°C
N_ELE_REF = 0.1368 #13.68%
N_ELE = N_ELE_REF*(1 - BETA_T*(Tcel - 25)) #273 + 25   - tcel  kelvin

A_tedlar = Av
L_TEDLAR = 0.33e-3
RHO_TEDLAR = 1500
M_TEDLAR = Av*L_TEDLAR*RHO_TEDLAR
CP_TEDLAR = 1090 #1090 OU 2090
K_TEDLAR = 0.35
ALPHA_TEDLAR = 0.34 #doc nasa ou zero

#parameters
RHO_ABS = 2700
A_ABS = Av
CP_ABS =900
L_ABS = 3e-3 #mm
M_ABS = A_ABS*RHO_ABS*L_ABS
K_ABS = 300
A_ABS_TUBO = 10*1.60*0.01+0.154*9*0.01
A_ABS_ISO = Av-A_ABS_TUBO

RHO_TUBO = 2700
CP_TUBO = 900
N_TUBOS = 10
L_TUBO = N_TUBOS*1.6
M_TUBO = RHO_TUBO*L_TUBO*(3.1415/4)*(Do**2 - Di**2)
K_TUBO = 300
A_TUBO_F = 0.387 #pi*Di*(L*10 VOLTAS + R(156MM)*9)
A_TUBO_ISO = 0.484 #pi*Do*(L*10 VOLTAS + R(156MM)*9)

A_ISO = Av
RHO_ISO = 50
L_ISO = 40e-3
M_ISO = A_ISO*RHO_ISO*L_ISO
CP_ISO = 670
K_ISO = 0.0375
E_ISO = 0.75 #ESTIMATED

RHO_FLUIDO = 997
M_FLUIDO = L_TUBO*(3.1415/4)*Di**2*RHO_FLUIDO
CP_FLUIDO = 4186 #j/kgK
MI_FLUIDO = 0.890e-3 #Pa*s ou N/m2 * s
K_FLUIDO = 0.607
M_PONTO = 0.05 #kg/s ou 0.5 kg/m3

#DIMENSIONLESS
Pr = CP_FLUIDO*MI_FLUIDO/K_FLUIDO #water 25°C
Re = RHO_FLUIDO*Vf*Di/MI_FLUIDO

if (Re<=2300):
    Nuf = 4.364
else:
    Nuf = 0.023*(Re**0.8)*(Pr*0.4)*Re

#COEFFICIENTS
h_rad_vidro_ceu = SIGMA*E_VIDRO*(Tv**2 - T_CEU)*(Tv + T_CEU)
h_conv_vidro_amb = 2.8 + 3*VENTO
h_conv_tubo_fluido = 0.5*30#Nuf
h_cond_vidro_cel = 1/((L_VIDRO/K_VIDRO) + (L_CEL/K_CEL))
h_cond_cel_tedlar = 1/((L_TEDLAR/K_TEDLAR) + (L_CEL/K_CEL))
h_cond_tedlar_abs = 1/((L_TEDLAR/K_TEDLAR) + (L_ABS/K_ABS))
h_cond_abs_tubo = 1/((L_TUBO/K_TUBO) + (L_ABS/K_ABS))
h_cond_abs_iso = 1/((L_ISO/K_ISO) + (L_ABS/K_ABS))
h_cond_tubo_iso = 1/((L_ISO/K_ISO) + (L_TUBO/K_TUBO))
h_conv_iso_amb = h_conv_vidro_amb
h_rad_iso_ground = SIGMA*E_ISO*(Tiso**2 - T_GROUND**2)*(Tiso + T_GROUND)

#GROUPS
A1 = (1/(M_VIDRO*CP_VIDRO))*(ALPHA_VIDRO*Av)#*G(t)) G_ajustado = polyval(dt,rev_curve)
A2 = (1/(M_VIDRO*CP_VIDRO))*(Av*(h_rad_vidro_ceu + h_conv_vidro_amb + h_cond_vidro_cel))
A3 = (1/(M_VIDRO*CP_VIDRO))*Av*h_cond_vidro_cel
A4 = (1/(M_VIDRO*CP_VIDRO))*Av*(h_conv_vidro_amb + h_rad_vidro_ceu)

A5 = (1/(M_CEL*CP_CEL))*(Pac*A_CELL*TAU_VIDRO*ALPHA_CELL) #*G(t)
A6 = -1*A5*N_ELE #*G(t)
A7 = (1/(M_CEL*CP_CEL))*A_CELL*h_cond_vidro_cel
A8 = (1/(M_CEL*CP_CEL))*A_CELL*(h_cond_vidro_cel + h_cond_cel_tedlar)
A9 = (1/(M_CEL*CP_CEL))*A_CELL*h_cond_cel_tedlar

A10 = (1/(M_TEDLAR*CP_TEDLAR))*A_tedlar*(1 - Pac)*TAU_VIDRO*ALPHA_TEDLAR#G(t)
A11 = (1/(M_TEDLAR*CP_TEDLAR))*A_tedlar*(h_cond_cel_tedlar + h_cond_tedlar_abs)
A12 = (1/(M_TEDLAR*CP_TEDLAR))*A_tedlar*h_cond_cel_tedlar
A13 = (1/(M_TEDLAR*CP_TEDLAR))*A_tedlar*h_cond_tedlar_abs

A14 = (1/(M_ABS*CP_ABS))*A_ABS*h_cond_tedlar_abs
A15 = (1/(M_ABS*CP_ABS))*(A_ABS*h_cond_tedlar_abs + A_ABS_TUBO*h_cond_abs_tubo + A_ABS_ISO*h_cond_abs_iso)
A16 = (1/(M_ABS*CP_ABS))*A_ABS_TUBO*h_cond_abs_tubo
A17 = (1/(M_ABS*CP_ABS))*A_ABS_ISO*h_cond_abs_iso

A18 = (1/(M_TUBO*CP_TUBO))*A_ABS_TUBO*h_cond_abs_tubo
A19 = (1/(M_TUBO*CP_TUBO))*(A_ABS_TUBO*h_cond_abs_tubo + A_TUBO_F*h_conv_tubo_fluido + A_TUBO_ISO*h_cond_tubo_iso)
A20 = (1/(M_TUBO*CP_TUBO))*A_TUBO_F*h_conv_tubo_fluido*0.5
A21 = (1/(M_TUBO*CP_TUBO))*A_TUBO_ISO*h_cond_tubo_iso

A22 = (1/(M_FLUIDO*CP_FLUIDO))*A_TUBO_F*h_conv_tubo_fluido
A23 = (1/(M_FLUIDO*CP_FLUIDO))*(A_TUBO_F*h_conv_tubo_fluido*0.5 + M_PONTO*CP_FLUIDO)
A24 = (1/(M_FLUIDO*CP_FLUIDO))*(T_F_ENT*(M_PONTO*CP_FLUIDO - h_conv_tubo_fluido*A_TUBO_F*0.5))

A25 = (1/(M_ISO*CP_ISO))*A_ABS_ISO*h_cond_abs_iso
A26 = (1/(M_ISO*CP_ISO))*(A_ABS_ISO*h_cond_abs_iso + A_TUBO_ISO*h_cond_tubo_iso + A_ISO*h_conv_iso_amb + A_ISO*h_rad_iso_ground)
A27 = (1/(M_ISO*CP_ISO))*A_TUBO_ISO*h_cond_tubo_iso
A28 = (1/(M_ISO*CP_ISO))*A_ISO*(h_conv_iso_amb*Tamb + h_rad_iso_ground*T_GROUND)

#DEFINE MODEL EQUATIONS - ODES -  (GLASS, PV CELL, TEDLAR, ABSORBER, TUBE, FLUID, INSULATION)   #  dT1dt = A1*G_ajustado - A2*x[0] + A3*x[1] + A4 # dT2dt = A5*G_ajustado - A6*G_ajustado + A7*x[0] - A8*x[1] + A9*x[2]# dT3dt = A10*G_ajustado - A11*x[2] + A12*x[1] +A13*x[3]

def SysEdo(x, k):#tv-x[0] tcel-x[1] ttedlar-x[2] tabs-x[3] ttubo-x[4] tiso-x[5] tfs-x[6] 
    dT1dt = A1*polyval(k,rev_curve) - A2*x[0] + A3*x[1] + A4
    dT2dt = A5*polyval(k,rev_curve) - A6*polyval(k,rev_curve) + A7*x[0] - A8*x[1] + A9*x[2]
    dT3dt = A10*polyval(k,rev_curve) - A11*x[2] + A12*x[1] +A13*x[3]
    dT4dt = A14*x[2] - A15*x[3] + A16*x[4] + A17*x[5]
    dT5dt = A18*x[3] - A19*x[4] + A20*x[6] + A20*T_F_ENT + A21*x[5]
    dT6dt = A22*x[4] - A23*x[6] + A24
    dT7dt = A25*x[3] - A26*x[5] + A27*x[4] + A28

    Tdot = np.array([dT1dt, dT2dt, dT3dt, dT4dt, dT5dt, dT6dt, dT7dt])

    return Tdot

#RungeKutta4
def RK4(f, x0, t0, tf, dt):
    
    t = np.arange(t0, tf, dt) #time vector
    nt = t.size #lenght of time vector
    nx = x0.size #length of state variables?
    x = np.zeros((nx,nt)) #initialize 2D vector
    
    x[:,0] = x0 #initial conditions

    
    #RK4 constants
    for k in range(nt-1):

        k1 = dt*f(t[k], x[:,k],k)
        k2 = dt*f(t[k] + dt/2, x[:,k] + k1/2, k)
        k3 = dt*f(t[k] + dt/2, x[:,k] + k2/2, k)
        k4 = dt*f(t[k] + dt, x[:,k] + k3, k)
        dx = (k1 + 2*k2 + 2*k2 + k4)/6

        x[:,k+1] = x[:,k] + dx

    return x,t

#Define problems
f = lambda t, x, k : SysEdo(x, k)

#initial state - t0 is initial time - tf is final time - dt is time step
x0 = np.array([30, 30, 30, 30, 30, 30, 30])
t0 = 0
tf = 1000
dt = 1

#EDO SOLVE
x, t = RK4(f, x0, t0, tf, dt)

plt.figure()
plt.plot(t, x[0], '-', label='Tvidro')
"""
plt.plot(t, x[1], '-', label='Tpv')
plt.plot(t, x[2], '-', label='Ttedlar')
plt.plot(t, x[3], '-', label='Tabs')
plt.plot(t, x[4], '-', label='Tiso')
plt.plot(t, x[5], '-', label='Ttubo')
plt.plot(t, x[6], '-', label='Tfsai')"""
plt.title('Gráfico')
plt.legend(['Tvidro', 'Tpv', 'Ttedlar', 'Tabs', 'Tiso', 'Ttubo', 'Tfsai'], shadow=False)
plt.xlabel('t (s)')
plt.ylabel('Temperatura (°C)') 
plt.xlim(0,20)
plt.ylim(0,150) 
plt.grid('on')
plt.show()

Thank you in advance, I am also open to completely start the implementation from scratch if there is a better way to do this with python or matlab.

score 0 · Answer 1 · answered Apr 27 '22 at 05:14

0

You can just replace

x, t = RK4(f, x0, t0, tf, dt)

with

t = arange(t0,tf+0.5*dt,dt)
res = solve_ivp(f,(t0,tf),x0,t_eval=t,args=(k,), method="DOP853", atol=1e-6,rtol=1e-8)
x = res.y[0]

Adapt the last 3 parameters to your liking.

answered Apr 27 '22 at 05:14

Lutz Lehmann

25,219
2
22
51

Dear @Lutz Lehmann, thanks for your reply. I changed the code like you suggest, but "args=(k,)" appears as warning - "k is not defined". K is my **for** iterator. I think i've to delete line 193 **"def RK4(f,x0,t0,tf,dt):** until line 221: **#EDO SOLVE**. But, when i do this, **f and k appears like "not defined".** Some suggestion? – Elias Apr 28 '22 at 00:12
Yes, but that just highlights a problem in your code. The `k` as argument should be the polynomial (coefficients array) that get evaluated with `polyval`, probably `G`. In RK4 the position is filled with the time index, which almost never makes sense. – Lutz Lehmann Apr 28 '22 at 04:32

System of seven ODEs solve using solve_ivp or implement RK4

1 Answers1