# Modified Bessel function second kind

13 views (last 30 days)
Erkan Ozkat on 30 Apr 2016
Commented: Walter Roberson on 1 May 2016
Hi everyone,
I have been working on the analytical solution on keyhole welding. And the first paper in this topic is A.H.F Kaplan 1994. He formulated heat flow inside keyhole in terms of 2nd kind modified Bessel function. However, in my implementation when r approaches zero (r=0), the ratio between Bessel function K1(0)/K0(0) returns NaN. Do you have any idea how to handle this singularity problem? I want to calulate qfront from -r to +r. Related paper is
Kaplan, Alexander. "A model of deep penetration laser welding based on calculation of the keyhole profile." Journal of Physics D: Applied Physics 27.9 (1994): 1805.
My implementation is
if true
model.Material.Density=7860*1e-9; % density (kg/mm^3)
model.Material.Cp=465; % specific heat capacity (J/(kg*K))
model.Material.Lamda=43/1000; % thermal conductivity (W/(mm*K))
model.Material.K=model.Material.Lamda/(model.Material.Density * model.Material.Cp); % thermal diffusivity (mm^2/s)
model.Process.Power=600; % W
model.Process.Speed=50; % mm/s
model.Material.Ta=300; % ambient temperature (K)
model.Material.Tm=1893; % melting temperature (K)
model.Material.Tv=3123; % vaporisation temperature (K)
% function q=getHeatFlowKeyHole(model, abs(r), phi)
Ta=model.Material.Ta; % ambient temperature (K)
Tv=model.Material.Tv; % vaporisation temperature (K)
lamba=model.Material.Lamda;
s=model.Process.Speed;
k=model.Material.K;
Pnumber=s/(2*k);
% Modified Bessel Function Second Kind
K0=besselk(0,Pnumber*r);
K1=besselk(1,Pnumber*r);
qfront=(Tv-Ta) * lamba * Pnumber * (1 + K1 / K0);
end
Thank you Erkan

John D'Errico on 1 May 2016
Edited: John D'Errico on 1 May 2016
1. Are you sure this generates something meaningful for negative r? Those bessel functions return complex results for negative r.
help besselk
besselk Modified Bessel function of the second kind.
K = besselk(NU,Z) is the modified Bessel function of the second kind,
K_nu(Z). The order NU need not be an integer, but must be real.
The argument Z can be complex. The result is real where Z is positive.
If I had to guess, the function is symmetric across zero, since when r is negative, I see a comment in your code about abs( r ), but then you never bother to use abs() in the code.
2. The limit of K1/K0 as r --> 0 from above appears to be +inf. If you try to evaluate those functions exactly at r==0, you have a 0/0 condition, which will be NaN.
Walter Roberson on 1 May 2016
Looking at the definitions of BesselK, it looks like it is definitely undefined at 0, with a (0)^(-1)/0 term