Modified Bessel function second kind

6 次查看（过去 30 天）
Erkan Ozkat
Erkan Ozkat 2016-4-30
评论： David Goodmanson 2025-4-15
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)
% read constants
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
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2025-4-13
Ran in:
@David Goodmanson it's more complicated than that...
syms z
K0 = besselk(sym(0),z)
K0 = 
K1 = besselk(sym(1),z)
K1 = 
subs(K0, z, sym(0))
ans = 
NaN
subs(K1, z, sym(0))
ans = 
NaN
but
besselk(sym(0), sym(0))
ans = 
∞
besselk(sym(1), sym(0))
ans = 
∞
I don't understand at the moment how they can give two different results depending on whether you subs() or not.
David Goodmanson
David Goodmanson 2025-4-15
Hello Walter,
I believe that my comment is correct when the bessel functions are calculated numerically, although that is switching gears from what the OP was asking in terms of syms. The limits
K1(z) --> 1/z and K0(z) --> -ln(z)
both work when z --> 0 from any direction in the complex plane. However, the K0 expression gets to the limit quite slowly. A better approximation along the way is
K0(z) --> -ln(z/2)-eulerC % eulerC = 0.577215664901533...
As to syms considerations, of course there are many questions on this site asking for a symbolic solution to equations so difficult that no such solution is possible. Can't blame syms for those. Other times though, the questions are more along the lines of your observation here, trying to determine why syms is doing something according to its own arbitrary, obscure symsian rules.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

John D'Errico
John D'Errico 2016-5-1
编辑：John D'Errico 2016-5-1
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.
  1 个评论
Walter Roberson
Walter Roberson 2016-5-1
Looking at the definitions of BesselK, it looks like it is definitely undefined at 0, with a (0)^(-1)/0 term

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Bessel functions 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by