Error when using if-else statement in MATLAB function block. How can I fix this error described below?

function m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
q_w = Q/(2*pi*rporous_outer*height_porous); %heat input converted to heat flux
if T_w>=T_sat
    m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
else
    m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x); 
    %the power of 0.4 is what triggers the complex number    
end

When I run this code, I get the following error, which is triggered by the line ' m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);

I've checked if any of the values result in a negative number, but none of the terms are negative, nor their result. I've run this separately, and this the same line works fine and I get a positive results without an complex number or any error? Why would this occur, does this error have to do anything with the if-else statement used? Any thoughts would be much appreciated.

I think this issue is triggered by the use of ^0.4, which cannot be changed. Why would an error/complex result occur in simulink and not in MATLAB?

Many thanks

[Sam Chak]

Hi @Delika Weragoda

Could you please provide us with the values of the parameters? This will allow us to test the function 'fcn()' and compute 'm'.

m = fcn(1, 10, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)

m = Inf

function m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
    q_w = Q/(2*pi*rporous_outer*height_porous); %heat input converted to heat flux
    if T_w>=T_sat
        m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
    else
        m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x); 
    end
end

Answer 1

[Sulaymon Eshkabilov]

Have you tried to change ^0.4 with: ((...)^2)^(1/5) or ((...)^2)^(0.2)

Answer 2

[Walter Roberson]

m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x); 
    1    2    345              4 5                        43 456        5       6    7       654    32 3   210 1 0
    

Simulink is calculating the expression

(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))

as being negative. When negative is raised to the fractional power 0.4 it results in a complex value.

My guess is that psi is outside the range -1 to +1

If you want a complex result, then code

m = (k_f*complex(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x); 

If you do not want a complex result, you will need to protect against it, such as

if abs(psi) >= 1
    m = nan;
elseif T_w>=T_sat
    m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
else
    m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x); 
end

[Delika Weragoda]

Hi @Sam Chak/ @Walter Roberson

@Walter Roberson none of the values are negative including psi, I am not sure whether its triggered by Q and T_w, but when I run the same line of code [i.e. m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x)], in the MATLAB function block without an if-else statement and the other condition, I don't get any of these error, works perfectly.

Please see the values of the parameters below,

Q - a signal taken from heat flow sensor block (constant in this case = 10)

T_sat = 307

T_w = taken from temperature sensor block (all positive)

sigma = 0.0002

rho_v = 7.1

h_fg M = 142e3

P_v = 9.15e4

R = 41.55

g = 9.81

beta = 0.002

x = 0.037

nu_novec = 3.2e-7

alpha_novec = 4.12e-8;

k_m = 3160

d_pore = 0.05

psi = 0.46

k_f = 0.075

rporous_outer = 0.0103

height_porous = 0.074

@Sulaymon Eshkabilov thank you for the suggestion, I've tried both does work.

[Sam Chak]

Hi @Delika Weragoda

The value for 'h_fg' and sensor data 'T_w' are not provided.

Q             = 10;
T_sat         = 307;
T_w           = 299:305;    % <-- this one
sigma         = 0.0002;
rho_v         = 7.1;
h_fg          =             % <-- and this one
M             = 142e3;
P_v           = 9.15e4;
R             = 41.55; 
g             = 9.81;
beta          = 0.002;
x             = 0.037;
nu_novec      = 3.2e-7;
alpha_novec   = 4.12e-8;
k_m           = 3160;
d_pore        = 0.05;
psi           = 0.46;
k_f           = 0.075;
rporous_outer = 0.0103;
height_porous = 0.074;

for j = numel(T_w)
    m(j) = fcn(Q, T_sat, T_w(j), sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous);
end
stem(m), grid on

function m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
    q_w = Q/(2*pi*rporous_outer*height_porous); %heat input converted to heat flux
    if T_w>=T_sat
        m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
    else
        m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x); 
    end
end

[Delika Weragoda]

Hi @Sam Chak Sorry, please see the values below.

Q             = 10; %This is taken as an input from the model (see below)
T_sat         = 307;
T_w           = 299:305;    % This is also taken as an input from the model (see below) 
sigma         = 0.0002;
rho_v         = 7.1;
h_fg          = 142e3            
M             = 0.2;
P_v           = 9.15e4;
R             = 41.55; 
g             = 9.81;
beta          = 0.002;
x             = 0.037;
nu_novec      = 3.2e-7;
alpha_novec   = 4.12e-8;
k_m           = 3160;
d_pore        = 0.05;
psi           = 0.46;
k_f           = 0.075;
rporous_outer = 0.0103;
height_porous = 0.074;

Please see the model for your reference

[Sam Chak]

@Delika Weragoda

Could you please verify if the m function returns the expected values? So far no errors.

%% Parameters
Q             = 10;
T_sat         = 307;
sigma         = 0.0002;
rho_v         = 7.1;
h_fg          = 142e3;
M             = 142e3;
P_v           = 9.15e4;
R             = 41.55; 
g             = 9.81;
beta          = 0.002;
x             = 0.037;
nu_novec      = 3.2e-7;
alpha_novec   = 4.12e-8;
k_m           = 3160;
d_pore        = 0.05;
psi           = 0.46;
k_f           = 0.075;
rporous_outer = 0.0103;
height_porous = 0.074;

%% Case 1a: 1 deg lesser than T_sat
T_w           = 306;
m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)

m = 22.4468

%% Case 1b: absolute zero
T_w           = 0;
m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)

m = 22.4468

%% Case 2a: 1 deg greater than T_sat
T_w           = 308;
m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)

m = 1.1852e+05

%% Case 2b: % speed of light
T_w           = 3e8;    
m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)

m = 1.1852e+05

%% The m-function
function m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
    q_w = Q/(2*pi*rporous_outer*height_porous); %heat input converted to heat flux
    if T_w>=T_sat
        m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
    else
        m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x); 
    end
end

[Delika Weragoda]

Hi @Sam Chak

Yes, these are the expected values. I do get them when I run the code in MATLAB, but not in Simulink.

I also tried simply using the the second m function (only) in the function block and it works fine (gets the results as above). However, I wanted to compare this with experimental results so I added an inport block to compare the two results and then the same error pops up, which is very weird.

I am thinking of upgrading MATLAB tomorrow, will let you know how it goes and if the error still persists.

Thanks

Delika

[Delika Weragoda]

@Sam Chak

I've managed to sort the issue, I've replaced a part of the function (in bold) using another variable (and running the replaced part on workspace) and the model works fine.

m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);

m = m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(q^(0.4)))))/(x)

Thanks for your help.

Best regards

Delika

[Sam Chak]

Hi @Delika Weragoda,

Thank you for updating us. It's great to hear that you were able to resolve the issue through your own efforts. You might want to consider posting the full solution in the Answer section below and then accepting your self-discovered solution to close the issue.

This could be beneficial for other Simscape users who may encounter a similar issue in the future.

[Delika Weragoda]

@Sam Chak So, apparently, I've made a mistake in the parenthesis, that is why the model supposedly worked, but I managed to find the issue. I believe this error in complex results in triggered from the heat sensors' input that I take to calculate the value for m, which is dynamically changing.

The variable m, is a function of Q, which is the input taken from the heat flow rate sensor. This is what is causing the error, as this signal will be a power of 0.4. Is there a way to overcome this?

function m = fcn(k_f, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, Q, rporous_outer, height_porous)

q_w = Q/(2*pi*rporous_outer*height_porous); %(heat flux)

    m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);

Please see the model below for your reference

[Sam Chak]

Hi @Delika Weragoda

m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);

I tend to divide long equations and conquer them later. This approach helps 100% prevent human mistakes such as mismatched parentheses. I suggest you rewrite the equation 'm' in multiple terms and then return here for a review.

%% parameters
a   = ...;
b   = ...;
c   = ...;
d   = ...;

%% terms
term1 = ...;
term2 = ...;
term3 = ...;

%% equation
m   = term1 + (term2/term3)^0.4 ...;

Answer 3

[Delika Weragoda]

@Sam Chak

I found the issue. There is a negative heat flow generated due to the thermal mass that I am using in the model. Why would something like this occur? I've tried debugging the model by monitoring the heat flow at different sections of the model. Please see the results below.

Model:

heat transfer subsystem

Heat flow results

The four heat flow sensors shown in the figure above are used to monitor the heat flow at the four locations. The negative values at the first two time steps is what triggers the error.

I would really appreciate if you can give any feedback on how to rectify this please.

Many thanks

Answer 4

[Sam Chak]

Hi @Delika Weragoda

Try using the Abs block to return the absolute value of the input.