How can i implemente this PDE boundary conditions?

I wrote this code using pdepe to implement this system (I wrote the constant terms like a, b, c, d, and for simplicity)

function[c,f,s]=heatsempccl(z,t,u,dudz)
global rhoc cp kc h eps q sigma Tinf volf volm rhom Hu mc n
c=[1;1];
f=[(kc/(rhoc*cp));0].*dudz;
A1=75000;       %[1/s]
E1=74690;       %[J/mol]
A2=21666.667;    %[1/s]
E2=58370;       %[J/mol]
R=8.314462;      %[J/(mol*K)]
k1=A1*exp(-E1/(R*u(1)));
k2=A2*exp(-E2/(R*u(1)));
dalfadt=(k1+k2*u(2)^mc)*(1+u(2))^n;
s=[((volm*rhom*Hu)/rhoc*cp)*dalfadt;dalfadt];
end

How can i write the boundary conditions d(alpha)/dz=0 at both ends?

[Torsten]

pr(2)=pl(2)=0

qr(2)=ql(2)=1

in the function routine where you define the boundary conditions.

[Claudio D'Amato]

function[pl,ql,pr,qr]=heatbcccl(zl,ul,zr,ur,t)
global rhoc cp kc h eps q sigma Tinf
pl=[q-h*(ul(1)-Tinf)-eps*sigma*(ul(1)^4-Tinf^4);0];
ql=[rhoc*cp;0];
pr=[0;0];
qr=[1;0];
end

I had tried to write it, but if i run the code, Matlab gives me this error message:

"Error using pdepe (line 293)

Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving

spatial derivative."

[Torsten]

You set ql(2)=qr(2)=0 instead of 1.

[Claudio D'Amato]

If i set ql(2)=qr(2)=1 he send me this message

"Warning: Failure at t=8.624748e+00. Unable to meet integration tolerances without reducing the step size below the

smallest value allowed (2.842171e-14) at time t."

My doubt in this case is this: if f is f=[(kc/(rhoc*cp));0].*dudz, this boundary conditions p(x,t,u)+q(x,t)*f(x,t,u,du/dx)=0 in this case are not 0=0?

[Bill Greene]

What is the mathematical definition for the boundary conditions on T? Your code for these looks very strange to me. Also, I suggest posting a complete listing of your current code.

[Torsten]

This doesn't matter since f is not used for the integration of alpha.

Did you look at the results up to t=8,624748 ?

[Claudio D'Amato]

This is a mathematical model of the cure of a composite material. The first equation is the Fourier heat equation, the second equation represents the cure kinetics of the composite (alpha is the degree of cure). The boundary conditions on T represent the heat flux balance on the surfaces. On the left side there is a balance regarding the heat flux supplied, the convective one, the radiative one and the conductive one, on the right side there is the condition of thermal insulation. I look the results up to t=8,624748 and they seem to me to be physically correct results. Here i post the whole code (sorry for the notes in Italian, I'm Italian).

global rhoc cp kc h eps q sigma Tinf volf volm rhom Hu mc n
rhoc=1452.696;     %densità del composito   [kg/m^3]
cp=1125.767;       %calore specifico del composito   [J/kg*K]
kc=0.354;       %conducibilità termica del composito   [W/m*K]
h=15;        %coefficiente convettivo [(W/m^2)*K]
eps=0.95;      %emissività del composito
q=5000;        %flusso di calore fornito [W/m^2]
Tinf=293;      %temperatura aria  [K]
sigma=5.670367e-8;         %costante di stefann-boltzmann[W/(m^2*K^4)]
volf=0.62;                %frazione in volume fibre
volm=1-volf;              %frazione in volume matrice
rhom=1050;                %%densità della matrice   [kg/m^3]
Hu=446000;                     %calore di reazione   [J/kg]
mc=1.216;                      %costante catalitica
n=1.284;                      %costante catalitica
zmesh=linspace(0,0.003,80);
tspan=linspace(0,1000,4000);
m=0;
options=odeset('RelTol',1e-10,'AbsTol',1e-10);
    u=pdepe(m,@heatsempccl,@heaticccl,@heatbcccl,zmesh,tspan,options);
    
function[c,f,s]=heatsempccl(z,t,u,dudz)
global rhoc cp kc h eps q sigma Tinf volf volm rhom Hu mc n
f=[(kc/(rhoc*cp));0].*dudz;
A1=75000;       %[1/s]
E1=74690;       %[J/mol]
A2=21666.667;    %[1/s]
E2=58370;       %[J/mol]
R=8.314462;      %[J/(mol*K)]
k1=A1*exp(-E1/(R*u(1)));
k2=A2*exp(-E2/(R*u(1)));
dalfadt=(k1+k2*u(2)^mc)*(1+u(2))^n;
c=[1;1/dalfadt];
s=[((volm*rhom*Hu)/rhoc*cp)*dalfadt;1];
end

function u0=heaticccl(z)
u0=[293;0];     
end

function[pl,ql,pr,qr]=heatbcccl(zl,ul,zr,ur,t)
global rhoc cp kc h eps q sigma Tinf
pl=[q-h*(ul(1)-Tinf)-eps*sigma*(ul(1)^4-Tinf^4);0];
ql=[rhoc*cp;1];
pr=[0;0];
qr=[1;1];
end

[Torsten]

Why do you use c(2)=1/dalphadt and s(2)=1 instead of c(2)=1 and s(2)=dalphadt ?

And in the definition of dalphadt, you use (1+u(2))^n instead of (1-u(2))^n.

Taking larger tolerances may also improve convergence.

[Claudio D'Amato]

perfect, all solved. thank you so much for your availability and patience!!!