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!