pdepe with discontinuos domain
1 view (last 30 days)
Show older comments
Hello All,
I am trying to implement a set of pdes (1-D, time dependent) having parameters with different values in different domains (for example, 0<x<17, 17<x<34). The interfaces of domain are total 6 in number. The code is working fine (without errors).
Problem: The obtained solution is not smooth at the given interface, although a lot of points of domain are defined near the interface.
The code is as follows. Any kind of help is appreciated.
Thanks and regards,
-Meenal
function abc
clear all; clc;
global Temperature L tend cp k rho wb qm cb Ta rho_b qs
L=0.093;
tend=3;
Ta=37;
cp=[3800; 3700; 3800; 3700; 3800; 2300; 4300];
k=[0.5; 0.5; 0.5; 0.5; 0.5; 1.16; 0.34];
rho=[1052; 1052; 1052; 1052; 1052; 1052; 1000];
wb=[0.00833; 0.0005; 0.00833; 0.0005; 0.00833; 0.0003; 0.00033];
qm=[25000; 25000; 25000; 25000; 25000; 3700; 3600];
cb=3800;
rho_b=1052;
qs=0;
m=0;
x = [linspace(0,17,40) linspace(17.01,34,40) linspace(34.01,51,40) linspace(51.01,68,40) linspace(68.01,85,40) linspace(85.01,89,10) linspace(89.01,93,10)];
t=linspace(0,tend,41);
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
Temperature = sol(:,:,1);
figure
surf(x,t,Temperature)
title('Numerical solution with nonuniform mesh')
xlabel('Distance x')
ylabel('Time t')
zlabel('Solution Temperature')
figure
plot(x,Temperature(end,:))
xlabel('Distance x')
ylabel('Solution Temperature')
title('Solution profiles at several times')
function [c,f,s] = pdex2pde(x,t,T,DTDx)
global cp k rho wb qm cb Ta rho_b qs
if x <= 17
c = rho(1).*cp(1);
f = k(1).*(DTDx);
s = wb(1).*cb*rho_b.*(Ta-T)+qm(1)+qs;
elseif (x>17 && x<=34)
c = rho(2).*cp(2);
f = k(2).*(DTDx);
s = wb(2).*cb*rho_b.*(Ta-T)+qm(2)+qs;
elseif (x>34 && x<=51)
c = rho(3).*cp(3);
f = k(3).*(DTDx);
s = wb(3).*cb*rho_b.*(Ta-T)+qm(3)+qs;
elseif (x>51 && x<=68)
c = rho(4).*cp(4);
f = k(4).*(DTDx);
s = wb(4).*cb*rho_b.*(Ta-T)+qm(4)+qs;
elseif (x>68 && x<=85)
c = rho(5).*cp(5);
f = k(5).*(DTDx);
s = wb(5).*cb*rho_b.*(Ta-T)+qm(5)+qs;
elseif (x>85 && x<=89)
c = rho(6).*cp(6);
f = k(6).*(DTDx);
s = wb(6).*cb*rho_b.*(Ta-T)+qm(6)+qs;
elseif (x>89 && x<=93)
c = rho(7).*cp(7);
f = k(7).*(DTDx);
s = wb(7).*cb*rho_b.*(Ta-T)+qm(7)+qs;
end
function T0 = pdex2ic(x)
T0 =309;
function [pl,ql,pr,qr] = pdex2bc(xl,Tl,xr,Tr,t)
pl = 0;
ql = 1;
pr = Tr - 312;
qr = 0;
0 Comments
Answers (0)
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!