- Are you looking into the equation1) rho * Cp * (du/dt) = k * (d^2u/dt^2) + s? and 2) deltaT/deltat = (s*L) / (h * A)?
- If we are rewriting the equation 2) like(Potf - Poti) / (tf - ti) = (s * L) / (h * A)then with the available value of Potf, Poti, tf, ti, h and A, you can find the value of s right?
I need to know how to put the s variable of pdefun (from pdepe solver) as a function of x and time
2 次查看(过去 30 天)
显示 更早的评论
JOAN PERE PONSETI
2017-1-6
Hi, I need help to solve a proble I'm in for a long ago. The thing is that I can't find any doc explaining how to set the variables [c,f,s] of pdefun (in particular the 's') as a function of the point of the 'xmesh'. It is to solve the heat transfer and difussion equation that is rho*cp*du/dt-k*d/dx(du/dx)=g.
from mathworks blogs and answers, I know that c=rho*cp (which is constant with time and x) f=k (also constant), but g is power generation, and that is not constant with time neither with x, the power is generated the first time values and just in a given set of points of xmesh.
Here I paste a part of the code where I state the pdefun:
function[c,f,s] = pdef(x, u, t, DuDx)
c = rho*cp;
f = k*DuDx;
i=1;
while(i<=length(x))
if(x(i)<=percent*radi)
s = (1/(h*pi*radi)^2)).*interp1(t40,Pot40,t,'linear','extrap');
else
s = 0;
end
i=i+1;
end
end
Please, if anyone could tell me if this is well developed. The s variable I mean.
8 个评论
Nikhil Sreekumar
2017-1-9
Hi Joan,
I had some doubts could you please clarify?
Thanks
Nikhil
JOAN PERE PONSETI
2017-1-10
Hi, yes, my doubt is regarding to the equation you call 1), but variable u in my case is temperature T. The thing is that, matlab says that s is a function of x,u,t,dudx and I just want it to put s(x,t) because the source term depends of the position within the radius (x) and the time t, with the power. My source term is the power density generation, which is P/volume and the power depends on time and the volume depends on x. this is why my source term depends on x and time.
thanks!
JP
Torsten
2017-1-11
The formulation
"s is a function of x,u,t,dudx"
means that s may depend on all of these variables, but of course it is possible that it depends on x and t only.
Best wishes
Torsten.
JOAN PERE PONSETI
2017-1-11
Yes, I know, but my problem is that I don't know exactly how the anonymous function of pdefun works, I mean, if the pdepe calls for the pdefun each time in 'tspan' and each x in 'xmesh', if so, pdefun returns 's' as a single value (double) but don't know how to give the function the vectors of tspan and xmesh and then return that single value.
Thanks!
JOAN PERE PONSETI
2017-1-11
The thing is that there is a velue for s in each time and each x, then, for a given x point there will be never the same value of s because it also changes on time. So, if I have a tspan of 20 intervals, and and xmesh of 50 points, s should be a matrix (50x20) which will have in each cell the value of the source term at each time and each x.
Torsten
2017-1-11
编辑:Torsten
2017-1-11
Why do you think you have to supply tspan and xmesh in "pdefun" ?
From the documentation of pdepe:
pdefun computes the terms c, f, and s (Equation 1-3). It has the form
[c,f,s] = pdefun(x,t,u,dudx)
The input arguments are scalars x and t and vectors u and dudx that approximate the solution u and its partial derivative with respect to x, respectively. c, f, and s are column vectors. c stores the diagonal elements of the matrix c (Equation 1-3).
So your code should read:
function[c,f,s] = pdef(x, u, t, DuDx)
c = rho*cp;
f = k*DuDx;
if x<=percent*radi
s = (1/(h*pi*radi)^2))*interp1(t40,Pot40,t,'linear','extrap');
else
s = 0;
end
end
Of course, the variables percent, radi, h, t40 and Pot40 must be visible in "pdef".
Best wishes
Torsten.
JOAN PERE PONSETI
2017-1-11
Hi again, I implemented this code to see if now it works but it gives me this error message back.
Not enough input arguments.
Error in CLDissipation1D/pdef (line 375) if(x==xmesh(1))
Error in pdepe (line 246) [c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});
Error in CLDissipation1D (line 340) sol = pdepe(m, @pdef, @pdeic, @pdebc, xmesh, t40);
function[c,f,s] = pdef(x,t,DuDx,S,xmesh,tspan)
c=rho*cp;
f=k*DuDx;
% find out the volume we need the previous point too
i=1;
encontrado=0;
while(encontrado==0)
if(x==xmesh(1))
vol=inf;
encontrado=1;
elseif(x==xmesh(i))
vol=(pi*(x^2-xmesh(i-1)^2)*h);
encontrado=1;
end
i=i+1;
end
% repeat the same loop to find the power at time t
i=1;
encontrado=0;
while(encontrado==0)
if(t==tspan(i))
pot=Pot40(i);
end
i=i+1;
end
% Finally we define the variable s as the ratio of power/volum
s=pot/(vol*pin);
% ========================================================================= % defines the Initial Conditions
end
Thanks a lot again!
采纳的回答
Torsten
2017-1-11
1. pdepe calls pdef for arbitrary times, not only those specified in tspan.
2. sol = pdepe(m, @(x,t,DuDx)pdef(x,t,DuDx,S,xmesh,tspan), @pdeic, @pdebc);
3. Comparisons for equality like if x==xmesh(1) probably won't work because of rounding errors.
Best wishes
Torsten.
9 个评论
JOAN PERE PONSETI
2017-1-13
Ok, I have tried this, but matlab returns this error.
Not enough input arguments.
Error in CLDissipation1D>@(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh) (line 296) sol = pdepe(m, @(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh), @pdeic, @pdebc, xmesh, t40);
Error in pdepe (line 246) [c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});
Error in CLDissipation1D (line 296) sol = pdepe(m, @(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh), @pdeic, @pdebc, xmesh, t40);
sol = pdepe(m, @(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh), @pdeic, @pdebc, xmesh, t40);
function[c,f,s] = pdef(x,t,DuDx,vol,xmesh)
c=rho*cp;
f=k*DuDx;
s=interp1(t40,Pot40,t,'linear','extrap')/(pin.*interp1(xmesh,vol,x,'linear','extrap'));
% ========================================================================= % defines the Initial Conditions
end
Thanks a lot!
JOAN PERE PONSETI
2017-1-13
And now, this. Here the problem is in pdeic??
Error using CLDissipation1D/pdeic Too many input arguments.
Error in pdepe (line 229) temp = feval(ic,xmesh(1),varargin{:});
Error in CLDissipation1D (line 188) sol = pdepe(m, @(vol)pdef,@pdeic, @pdebc, xmesh, t40,[], vol);
sol = pdepe(m, @(vol)pdef,@pdeic, @pdebc, xmesh, t40,[], vol);
u=sol(:,:,1);
%========================================================================= % defines the properties of the PDE function
function[c,f,s] = pdef(x,t,u,DuDx,vol)
c=rho*cp;
f=k*DuDx;
s=interp1(t40,Pot40,t,'linear','extrap')/(pin.*interp1(xmesh,vol,x,'spline','extrap'));
% ========================================================================= % defines the Initial Conditions
end
function[u0] = pdeic(x)
u0 = tinitial;
end % ========================================================================= % defines the Boundary Conditions
function[pl,ql,pr,qr] = pdebc(xl, ul, xr, ur, t)
%for cylindrical pl and ql are 0, because of symmetry. No flux is present
%at the center of the sphere.
pl = 0;
ql = 0;
pr = ur-interp1(t40,Temp40,t,'linear','extrap');
qr = 0;
end
Torsten
2017-1-16
Take a look at this page on how to pass extra parameters by the anonymous function method:
https://de.mathworks.com/help/optim/ug/passing-extra-parameters.html
Best wishes
Torsten.
JOAN PERE PONSETI
2017-1-21
Thanks a lot! I've solved this by stating the variables as global. But now, I have another problem with the solution I get. The thing is that, there is a power dissipation at the center, and I have the initial conditions and the boundary ones. Logic says that if there is power dissipation inside, the temperature within the core must be always higher than the one outside (at the same time for both temperatures). But what I get is the graph below, and furthermore, the blue-crossed line of the graph doesn't start at the 300 K (aprox) which is the initial condition.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/185519/image.png)
I give you the code I have used and let's see if you can figure out what's the error.
Thanks!
Joan Pere
It looks like the differential is not working properly.. don't know ehat happens
t40=times{var};
Temp40=temperatures{var};
Pot40=potencies{var};
% Assignem els vectors a variables que es mostraran en el Workspace
assignin('base','t40',t40);
assignin('base','Temp40',Temp40);
assignin('base','Pot40',Pot40)
% figure
% plot(t40,Pot40);
%
% figure
% plot(t40,Temp40);
% capacitor properties
k = 1.5; % W/(m K) 30 Alumina %aluminum 237
rho = 3000; % alumina density 4000 Aluminum 2300
cp = 900; % J/K/kg Aluminum 910 Alumina 880
%capacitor geometry
h=0.01;
radi=0.005;
percent=0.25;
tinitial = Temp40(1); % intial temp throughout the tile
m = 1; % assume cylindrical symmetry
xmesh = linspace(0, radi, 200); % 100 punts a discretitzar
assignin('base','xmesh',xmesh)
global vol
for i=1:length(xmesh)
if(i==1)
vol(i)=pi*h*(radi)^2;
else
vol(i)=inf;
end
end
assignin('base','vol',vol)
sol = pdepe(m, @pdef,@pdeic, @pdebc, xmesh, t40);
T1=sol(:,:,1);
assignin('base','T1',T1)
% A surface plot is often a good way to study a solution.
figure;
subplot(2,1,1)
surf(xmesh,t40,T1);
title(sprintf('SOLUCIÓ NUMÈRICA AMB 200 PUNTS.\nArxiu %s.', VCC40));
xlabel('Distància al centre (m)');
ylabel('Temps (segons)');
% A solution profile can also be illuminating.
subplot(2,1,2)
plot(t40,T1(:,1),'x');
hold on
plot(t40,T1(:,end),'+');
hold on
plot(t40,Temp40,'o');
legend('Temperatura al centre','Temperatura externa','Condicions Incials','Location','best');
title(sprintf('PERFIL DE TEMPERATURES EN 2 PUNTS (x=0 i x=R).\nArxiu %s.', VCC40));
xlabel('Temps (segons)');
ylabel('Temperatura (Kelvin)');
%========================================================================== % defines the properties of the PDE function
function[c,f,s] = pdef(x,u,t,DuDx)
c=rho*cp;
f=k*DuDx;
s=(1/interp1(xmesh,vol,x,'linear','extrap')).*interp1(t40,Pot40,t,'linear','extrap');
end
% ========================================================================= % defines the Initial Conditions
function[u0] = pdeic(x)
u0 = tinitial;
end
% ========================================================================= % defines the Boundary Conditions
function[pl,ql,pr,qr] = pdebc(xl, ul, xr, ur, t)
%for cylindrical pl and ql are 0, because of symmetry. No flux is present at the center of the sphere.
pl = 0;
ql = 0;
pr = ur-interp1(t40,Temp40,t,'linear','extrap');
qr = 0;
end
The volume is needed because my source term is the power density (Power/Volume) then the power is a different one each time (that's why I have to interpolate) and the volume I state that the volume dissipation is a quarter of the radius (percent=0.25)
Torsten
2017-1-23
1. In the boundary part (@pdebc), pl=ql=0 can't be correct.
2. I don't see where "percent" is used in the calculation of the source term.
3. I don't understand your variable "vol" because you assign "inf" to all the inner grid points.
Best wishes
Torsten.
JOAN PERE PONSETI
2017-2-24
1. pl and ql must be zero because there is no te flux into that point, which is the center of the cilindre.
2. the percent variable let's say is irrelevant here.
3. I see now that the volume variable code is not complete, I might have deleted some code lines. But the thing is that, i need the inner points to have a given value (power/volume, both of them are scalars) but i need the outter ones to be 0, thats why I set de volume to inf, since it divides and makes the value of the source term 0 (power always have a positive value).
Thanks!
JOAN PERE PONSETI
2017-2-24
编辑:JOAN PERE PONSETI
2017-2-24
Hello again, I need some help in another code, similar to this one but now using the solvepde function. When I have to assign the coefficients. I guess that when I write the following code:
sourcefun= @(region,state) interp1(t40,Pot40,state.t,'linear','extrap')/(volumeD*count);
specifyCoefficients(pdem,'m',0,'d',rho*cp,'c',k,'a',0,'f', sourcefun);
Matlab returs an error that says the following:
Error using pde.CoefficientAssignment/checkCoefFcnHdlArgCounts (line 531)
Function handle specifying a coefficient must accept two
input arguments and return one output argument.
As if the function doesn't return any value, is it parmited to use the state and region variables in a function like interp1??
Thanks a lot!!
Joan Pere
JOAN PERE PONSETI
2017-2-24
Another thing, when I write the next code % Specify Coefficients
f = @(region,state) interp1(t40,Pot40,state.time,'linear','extrap');
specifyCoefficients(pdem,'m',0,'d',rho*cp,'c',k,'a',0,'f', f);
and I set a breakpoint at f, it always says that state structure is 0 at all its elements but it doesn't return any error until the fifth or sixth step... don't know why...
thanks!
JClarcq
2018-2-19
Hi Joan,
I have the same problem than you where my heatsource is a term temperature dependent via an interp1 of a material properties.
Have you been able to solve this inputs argument issue? Thanks,
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)