How to use the PDE Toolbox combined with the script?
14 次查看(过去 30 天)
显示 更早的评论
Hello, is there anybody who can help me with my question? I need to solve a 2D heat diffusion equation PDE, I have the boundary conditions and initial conditions set in an excel file, I also have them in a script. So, as the boundary conditions vary with time but I haven't got any expression for that but the file, I need to import this data from the script or excel file to PDE toolbox.
The thing is the following:
- Boundary Conditions: ExternalTemperature(t) is a vector of a given length containing the external temperatures of a cylinder, so at each time, the bc's change.
- Initial Conditions: IC(0), is the first temperature in that previous vector. ExternalTemperature(1).
There are two zones at the cylinder, and outer one (external crown) which doesn't generates power, so the source term is 0, and an inner one which does generate power, and has a given value for the source term. So, I need to set the coefficients as a function of time (power is a function of time) and space (x-y) since the points outside the inner cylinder doesn't generate power.
Thanks in advance,
Joan Pere!
PS: this is for a final degree project, so any help is gratefully welcomed!
2 个评论
Sam McDonald
2017-3-6
Joan Pere,
It is difficult to address your question without additional information. Could you provide some more detail on where you are getting stuck? Is it loading the data into MATLAB, or using it in the PDE Toolbox to define the boundary conditions? Do you have some example code you are writing or modifying to solve this problem?
In the meantime, if you have not already, take a look at these links from the MathWorks documentation:
采纳的回答
Alan Weiss
2017-3-8
It is possible that you are getting bit by a bug. Sorry. I believe that you will have no problems if you create your circles using a geometry function. For example,
function [x,y] = multicircle(bs,s)
% Create circles centered at (0,0) using four segments per circle
switch nargin
case 0
x = 6*4; % four edge segments per radius
return
case 1
A1 = repmat([0,pi/2,pi,3*pi/2],1,6); % start parameter values
A2 = repmat([pi/2,pi,3*pi/2,2*pi],1,6); % end parameter values
A3 = repmat(1:6,4,1);
A3 = A3(:)'; % region label to left
A4 = repmat([2:6,0],4,1);
A4 = A4(:)'; % region label to right
A = [A1;A2;A3;A4];
x = A(:,bs); % return requested columns
return
case 2
if numel(bs) == 1
bs = bs*ones(size(s));
end
ind = floor((bs-1)/4) + 1; % circle number from 1 through 6
x = ind.*cos(s);
y = ind.*sin(s);
end
I tested this geometry using the following script:
model = createpde();
geometryFromEdges(model,@multicircle);
generateMesh(model);
pdegplot(model)
pdemesh(model)
axis equal
I hope that this helps. Sorry about the bug.
Alan Weiss
MATLAB mathematical toolbox documentation
11 个评论
Alan Weiss
2017-3-22
I think that you are making a serious mistake in not defining subdomains. The mesh will not respect the boundaries that are in your head, and therefore the elements can overlap the regions. I mean your mesh will include points from more than one type of region, no matter how fine your mesh.
Also, I have trouble understanding your code because it seems to be written with the assumption that the regions are in some kind of order. I don't see how you are passing in the vectorlimits data, or t40 or Pot40. I think that this line is erroneous:
if(a(l)>=vectorlimits(j) && a(l)<=vectorlimits(j+1))
You should decide whether you want to write code in a vectorized fashion or in a loop. It is faster to write vectorized code, but until you really know what you are doing, write loops. I assume that vectorlimits is monotone increasing.
for i = 1:L % I prefer upper-case L, l looks like 1
indx = find(vectorlimits - a(i) > 0,1) % find the band
if mod(indx,2) = 1 % do your calculation here
f(i) = ...
else
f(i) = ...
end
end
I hope that this helps,
Alan Weiss
MATLAB mathematical toolbox documentation
Robin
2023-8-10
Hello Alan,
I am having a similar problem. In my case i want to create several circular structures in different locations. My problem is that I don't know how to modify your code. If I'm understanding the code correctly, I just wanted to adapt the last part under case 2. So my idea was to shift the cos or sin function in my desired direction.
x = ind.*cos(s + 1);
y = ind.*sin(s - 1);
But this results in the following error: Error using test_coefficient_def_for_geometry
Meshing failed due to invalid geometry. Each face must have a unique face ID.
Could you elaborate on this?
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Geometry and Mesh 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!