Nonlinear Bessel-ish Function Need Help Solving

2 views (last 30 days)
Kyle Marquis on 25 Aug 2020
Commented: David Goodmanson on 25 Aug 2020
Hi there, I'm looking to solve a nonlinear differential equation that is very similar to a bessel function, but with an added T^4 kick. where g is a constant, and which simplifies using product rule to I tried solving it using bvp4c with the boundary conditions before I knew it was similar to a bessel function. (I used this approach for a very similar problem with no issues) where T sub L is just a positive number specified in the code below. However, bvp4c told me there was an issue with a singular jacobian matrix. Code listed below. Is there something I'm doing wrong here, or is this because bessel functions cause problems for this sort of analysis?
I'm wondering if there's a better way to solve this differential equation using modified bessel functions, because I noticed my equations were similar to those in this youtube tutorial for fin convection heat transfer, and this paper, although I cannot find either of these seperate equations sourced anywhere else, so if there's names for those variations of bessel equations, I would appreciate if you could let me know how to search them up. Regardless, would it be possible to use a shooting algorithm to guess and check values for T, since my rightmost term is T^3 away from being a normal bessel problem? How would I best implement this using Matlab's bessel functions?
Thanks
% Pass variables into functions. Need global variables, although I know
% this is bad practice, I do not currently know how to work around this
global Tb;
global Tr;
global g1;
% Constants
k = 1300; % [W/m/K]
Sigma_Boltzmann = 5.670374419*10^(-8); % [W/(m^2*K^4)] Stefan Boltzmann constant
a = 0.0255/1000;
emissivity = 1;
b = 0.0001; %[m]
Tb = 300; %[K]
% Parameters
L = 0.4; %[m] length of radiator
Tr = 157.74801710; %[K] estimated temp at x0=0, want to iterate this until you get dT/dx0 = 0 at x0=0
N = 40; % Number of points in mesh
g1 = Sigma_Boltzmann*emissivity*L/a/k; % NOTE NO 2 IN FRONT BECAUSE ASSUME RADIATES ONE DIRECTION ONLY
% Boundary Value Problem
x0mesh = linspace(0,L,N);
solinit = bvpinit(x0mesh, @guess);
sol = bvp4c(@bvpfcn, @bcfcn, solinit);
x0vec = sol.x;
Tvec = sol.y(1,:);
Tpvec = sol.yp(1,:);
Tppvec = sol.yp(2,:);
format long
dTdx0atL = Tpvec(end)
qdotstart = k*a*b*Tpvec(1)
figure(1),clf
plot(x0vec, Tvec, '-o')
xlim([0,L])
ylim([0 Tb])
xlabel('x0 distance (m)'), ylabel('T temperature of radiator (K)')
title('Temperature vs Distance Plot')
set(1,'Position',[200 200 800 500])
% FUNCTIONS
function dTdx0 = bvpfcn(x0,T) % equation to solve
global g1;
dTdx0 = zeros(2,1);
dTdx0 = [T(2);+g1*T(1).^4/x0 - T(2)/x0];
end
%--------------------------------
function res = bcfcn(TaR,TaL) % boundary conditions at x0=0 and x0=L respectively
global Tb;
global Tr;
res = [TaR(1)-Tr; TaL(1)-Tb];
end
%--------------------------------
function g = guess(x0) % initial guess for y and y'
g = [-250*x0+300;5*x0];
end
%--------------------------------
David Goodmanson on 25 Aug 2020
Hi Kyle,
I'm afraid the nonlinearity is going to do in any kind of useful relation to bessel functions. Case in point, the function (9*g*x)^(-1/3) is a solution, although not one that satisfies the boundary conditions.

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!