Simulation of laser phase transformation using pde thermal model toolbox

Dear Matlabers,

I try to simulate the temperature time course at the orbitrary position during laser hardening using rectangular laser. I choose pde thermal model toolbox as tool. The most important thing is to simulate the austenitization during heating and no austenitization during cooling (which is theoritically using temperature gradient as judgement) However, it doesnt works correctly. I hope someone it can fix this problem. Most importantly, it can be used to simulate different customization during heating and cooling for phase transformation and melting process. Thank you for your comments!

clear
clear global variable
close all
dbstop if error
tic % start time
%%
% global variables
global v distance A
A=0.36; %Absorption[-], 0.36@(1030 nm 20°C grinded surface AISI4140)[BIAS]
%A=3.62e-5*state.u+0.3494; %increase from 0.36 at 20°C to 0.41 at 1400°C
v=20e-3; %[m/s] scanning speed
distance=4e-3; %[m] moving distance 10 mm
t=distance/v; % [s] time for scanning speed of 20 mm/s passing length of 10 mm 
Tr=20;%Tr room temperature in degC
W=7E-3/2;  %half of geometry width [m] x-achse  original 2*W=100 mm
D=8E-3/2;  %half of geometry depth [m] y-achse  original 2*D=10 mm
H=4E-3;    %geometry height  [m] z-achse   orginal H=13 mm 
%% Model building boundary heat flux
% gm=multicuboid(2*W,2*D,H); %multicuboid(W,D,H); [m]
gm=multicuboid(2*W,2*D,[H-0.5e-3 0.1e-3 0.1e-3 0.1e-3 0.1e-3 0.1e-3],'ZOffset',[0 H-0.5e-3 H-0.4e-3  H-0.3e-3  H-0.2e-3  H-0.1e-3]); %multicuboid(W,D,H); [m]
thermalmodel = createpde('thermal','transient'); %transient thermal model done
thermalmodel.Geometry=gm; 
pdegplot(thermalmodel,'FaceLabel','on','FaceAlpha',0.5,'CellLabels','on','FaceAlpha',0.5);%face label and transparency of 3d 
set(gca,'xlim',[-W W],'ylim',[-D D],'zlim',[0 H]);
xlabel 'x in m'
ylabel 'y in m'
zlabel 'z in m'
title('Fig.1: Geometry with Face Labels')
msh=generateMesh(thermalmodel,'Hmax',0.5e-3,'Hmin',0.1e-3,'Hgrad',2,'GeometricOrder','quadratic')  %meshing [m]; should use stacked cuboid
figure;
pdemesh(thermalmodel)
set(gca,'xlim',[-W W],'ylim',[-D D],'zlim',[0 H])
xlabel 'x in m'
ylabel 'y in m'
zlabel 'z in m'
title('Fig.2: Meshing of workpiece')
% pdemesh(thermalmodel,'NodeLabels','on')
%% 
%showing heat flux
% x=-W:W/50:W;  %[m]
% y=-D:D/50:D;  %[m]
% t1=1*t;
% [X,Y]=meshgrid(x,y);
% location.x = X;
% location.y = Y;
% state.time=t1;
% F = flux1(location,state);   %with time-dependent heat source
% figure
% surf(x,y,F)
% d=colorbar;
% d.Label.String = 'laser intensity in W/m^2';
% xlabel 'x in m'
% ylabel 'y in m'
% zlabel 'intensity in W/m^2'
% title('Fig.3: Laser intensity distribution of rectangular beam') 
%%
%boundary condition: make sure all the surface has heat flux, in order to
%make calculation of each time step converge
thermalBC(thermalmodel,'Face',27,'HeatFlux',@flux1);%Convectioncoeffi: 1000 [W/m/C] forced gas,Tamb=20°C; emissivity=absorptance (assumed as gray body)
thermalBC(thermalmodel,'Face',27,'ConvectionCoefficient',1000,'AmbientTemperature',Tr);
fac=1:gm.NumFaces;
fac([2 7 12 17 22])=[]; %faces inside the object, which should be removed
thermalmodel.StefanBoltzmannConstant = 5.670373E-8;
thermalBC(thermalmodel,'Face',fac,'Emissivity',A,'AmbientTemperature',Tr);
% thermalBC(thermalmodel,'Face',fac,'Temperature',Tr);
% initial condition
thermalIC(thermalmodel,Tr); % [degC] initial temperature 20°C
%%
%time step
device.SolverOptions.ReportStatistics='on';
dt=1e-3; %1 ms [s] for 20mm/s    
% TLIST=[0 logspace(-4,-3,20) 2e-3:dt:t];
TLIST=0:dt:t;
%%
% thermal properties@[1998 Thomas] [Edelstahlwerke]
% k=45;
% cp=470; %[J/kg/K] @20°C
% rho=7720; %kg/m^3 Assumed T independent
thermalProperties(thermalmodel,'ThermalConductivity',@kFunc, ...
                                   'MassDensity',@rhoFunc, ...
                                   'SpecificHeat',@cpFunc)   %thermal properties
%% method 2 direct calculation
R=solve(thermalmodel,TLIST); 
Tmin=20; %[Kelvin], initial temperature 20°C/293:15 K
Tmax=max(R.Temperature(:,end)); 
%%
%temperature and heating/cooling rate at surface center F2 and in depth as a function of time CHECKED 05052020
[gradTx,gradTy,gradTz] = evaluateTemperatureGradient(R,[0;0;H],1:numel(TLIST));
u0= interpolateTemperature(R,[0;0;H],1:numel(TLIST)); 
u1= interpolateTemperature(R,[0;0;H-0.1e-3],1:numel(TLIST)); 
u2= interpolateTemperature(R,[0;0;H-0.2e-3],1:numel(TLIST)); 
u3= interpolateTemperature(R,[0;0;H-0.3e-3],1:numel(TLIST)); 
u4= interpolateTemperature(R,[0;0;H-0.4e-3],1:numel(TLIST)); 
u5= interpolateTemperature(R,[0;0;H-0.5e-3],1:numel(TLIST));
figure; 
plot(TLIST,u0);
hold on
plot(TLIST,u1);
plot(TLIST,u2);
plot(TLIST,u3);
plot(TLIST,u4);
plot(TLIST,u5);
grid on
legend('surface','100 µm depth','200 µm depth','300 µm depth','400 µm depth','500 µm depth')
title 'Fig.8: transient temperature at surface center F2'
xlabel 'time in second'
ylabel 'temperature in degC'
%% integral transformation time   rewrite use recycle
% t2=tlist(find(u2)>ac1);
% U2=u1(find(u2)>ac1);
% I2=trapz(t2,U2); %integral transition time in 100 µm depth
% 
%% please ignore
% data=[TLIST',u0',u1',u2',u3',u4',u5'];
% [m,n]=size(data);
% data_cell=mat2cell(data,ones(m,1),ones(n,1));
% %title={'time','surface Kelvin','100µm in degC','200µm in degC','300µm in degC'};
% result=(data_cell);
% %result=[title; data_cell];
% s=xlswrite('20mms_fp_0.36x2670w_a1.5mm_M0.75-0.1mm_multilayer_dtvar_cphkSFB.xlsx', result);
%%
toc %stop time
disp(['running time: ',num2str(toc)])
%%          
%--------------------------------------------------------------------------
%function
function [QFlux] = flux1(location,state)  % laser input with line laser CHECKED 05052020
    global v distance A
    P=2670; %[W],laser input power, change fa
    a=1.5E-3/2; %[m],half of laser width, 1/e,15072020 based on fit form of measured intensity distribution using 600 µm fiber in focal position
    b=10E-3/2; %[m],half of laser length, assumed 1/e 
       %QFlux=A*P/pi/a/b*exp(-(location.x+10E-3-v*state.time).^2/a^2-(location.y).^2/b^2); % (-b<=y<=b) gaussian source
        QFlux=A*P/2/sqrt(pi)/a/b*exp(-(location.x+0.5*distance-v*state.time).^2/a^2);  % (-b<=y<=b) rectangular gaussian source [Arata 1978] start position 10 mm
        %QFlux=A*P/2/sqrt(pi)/a/b*exp(-(location.x).^2/a^2);  % (-b<=y<=b) rectangular gaussian source [Arata 1978] start position 10 mm
end
% step function
function [s, S] = sFunc(location,state)
C=2.84e20; % 1/s
na=1.525; % [-]
dH=4.185;%enthalpy of activation: 4.4 eV for FP; 4.0 eV for martensite, 4.185 eV from [Mio05]
k=8.62e-5; %bolzmann constant [eV/K]
temp=state.u+273;%state.u in degC, temp in Kelvin
%temp=state.u;
bT=C*exp(-dH./(k*temp));
const=bT.*state.time;
s=1-exp(-const.^na);
S=(C*dH*na*state.time.*exp(-(C*state.time.*exp(-dH./(k*temp))).^na).*exp(-dH./(k*temp)).*(C*state.time.*exp(-dH./(k*temp))).^(na - 1))./(k*temp.^2);%d(fa)/du
%FAt=C*na*exp(-(C*t.*exp(-dH./(k*u))).^na).*exp(-dH./(k*u)).*(C*t.*exp(-dH./(k*u))).^(na - 1);%d(fa)/dt
end
% thermal conductivity
function k1 = kFunc(location,state)
[s,~]=sFunc(location,state);
kp=-5e-5*state.u.^2+0.0199*state.u+36.149; %W/(mK), state.u in degC,FP Friedhelm SFB136 07102020
%kp=1E-5*state.u.^2-0.0278*state.u+46.558;%W/(mK); temperature in degC [2007 Lakhkar] 07102020
% kp=-3E-05*state.u.^2 - 0.0012*state.u + 39.262;%FP temperature in degC [2005 Miokovic]
ka=0.0066*state.u+19.51;%W/(mK), stateu in degC,Austenite Friedhelm SFB136 07102020
%ka=30;%W/(mK); Austenite temperature in degC [2007 Lakhkar]
% ka=20;%W/(mK); Austenite temperature in degC [2005 Miokovic]
k1 =(kp.*(1-s)+ka.*s).*double(state.ux<1e3)+ka.*double(state.ux>1e3);%W/(mK), stateu in degC,FP to austenite, Friedhelm SFB136
end
% heat capacity with latent heat
function cp = cpFunc(location,state)  
[s,S]=sFunc(location,state);
cpp=-8.6654E-09*state.u.^4 + 1.2813E-05*state.u.^3 - 5.3581E-03*state.u.^2 + 1.0587*state.u + 4.2794E+02;% [J/kg/K]FP temperature in degC Friedhelm; SFB136 07102020
%cpp=-1.43e-8*state.u.^4+3.99e-5*state.u.^3-0.0395*state.u.^2+16.844*state.u-2109;%[J/kg/K]FP; temperature in Kelvin [2007 Lakhkar]
%cpm=cpp;%[J/kg/K]martensite temperature in degC Friedhelm; SFB136 07102020
cpa=0.1134*state.u+504.94;%[J/kg/K]austenite temperature in degC Friedhelm; SFB136 07102020
%cpa=607;%[J/kg/K]austenite; temperature in Kelvin [2007 Lakhkar]
Hp=1000*(-8.5234E-09*state.u.^4 + 1.2898E-05*state.u.^3 - 5.7136E-03*state.u.^2 + 1.2636*state.u - 3.6727E+01);% [J/kg]FP temperature in degC Friedhelm; SFB136 07102020
Ha=1000*(0.556*state.u+66.391); % [J/kg]austenite temperature in degC Friedhelm; SFB136 07102020
cppa=(cpp.*(1-s)+cpa.*s+(Ha-Hp).*S).*double(state.ux<1e3)+cpa.*double(state.ux>1e3);%seems endless calculation
%cppa=(1-s).*cpp+s.*cpa+10e6*S;% [J/kg/K]apparent heat capacity FP to austenite temperature in Kelvin [1987 Civan]
%cpma=(1-s).*cpm+s.*cpa+(Ha-Hm).*S;% [J/kg/K]apparent heat capacity martensite to austenite temperature in Kelvin [1987 Civan]
cp=cppa;
end
% density
function rho = rhoFunc(location,state)
[s,~]=sFunc(location,state);
rhop=-7E-05*state.u.^2 - 0.2819*state.u + 7840.3;%[kg/m^3]FP temperature in degC Friedhelm; SFB136 07102020
rhoa=-0.5694*state.u+8079.5;% [kg/m^3]Austenite temperature in degC Friedhelm; SFB136 07102020
rho=(rhop.*(1-s)+rhoa.*s).*double(state.ux<1e3)+rhoa.*double(state.ux>1e3);%seems endless calculation
end
%