В этом примере показано, как выполнить 3-D переходный анализ проводимости тепла полой сферы, сделанной из трех различных слоев материала.
Сфера подвергается неоднородному внешнему потоку тепла.
Физические свойства и геометрия этой проблемы описаны в Сингхе, джайне, и Rizwan-uddin (см. Ссылку), который также имеет аналитическое решение для этой проблемы. Внутренняя поверхность сферы имеет температуру нуля в любом случае. Внешнее полушарие с положительным значение имеет неоднородный поток тепла, заданный
и азимутальные и углы возвышения точек в сфере. Первоначально, температура во всех точках в сфере является нулем.
Создайте тепловую модель для переходного теплового анализа.
thermalmodel = createpde('thermal','transient');
Создайте многослойную сферу с помощью multisphere
функция. Присвойте получившуюся геометрию тепловой модели. Сфера имеет три слоя материала с полым внутренним ядром.
gm = multisphere([1,2,4,6],'Void',[true,false,false,false]);
thermalmodel.Geometry = gm;
Постройте геометрию и покажите метки ячейки и столкнитесь с метками. Используйте FaceAlpha
из 0,25 так, чтобы метки внутренних слоев отобразились.
figure('Position',[10,10,800,400]); subplot(1,2,1) pdegplot(thermalmodel,'FaceAlpha',0.25,'CellLabel','on') title('Geometry with Cell Labels') subplot(1,2,2) pdegplot(thermalmodel,'FaceAlpha',0.25,'FaceLabel','on') title('Geometry with Face Labels')
Сгенерируйте mesh для геометрии. Выберите размер mesh, который достаточно крупен, чтобы ускорить решение, но достаточно прекрасный, чтобы представлять геометрию обоснованно точно.
generateMesh(thermalmodel,'Hmax',1);
Задайте теплопроводность, массовую плотность и удельную теплоемкость для каждого слоя сферы. Свойства материала являются безразмерными значениями, не данными реалистическими свойствами материала.
thermalProperties(thermalmodel,'Cell',1,'ThermalConductivity',1, ... 'MassDensity',1, ... 'SpecificHeat',1); thermalProperties(thermalmodel,'Cell',2,'ThermalConductivity',2, ... 'MassDensity',1, ... 'SpecificHeat',0.5); thermalProperties(thermalmodel,'Cell',3,'ThermalConductivity',4, ... 'MassDensity',1, ... 'SpecificHeat',4/9);
Задайте граничные условия. Самая внутренняя поверхность имеет температуру нуля в любом случае.
thermalBC(thermalmodel,'Face',1,'Temperature',0);
Наружная поверхность сферы имеет внешний поток тепла. Используйте функциональную форму тепловых граничных условий, чтобы задать поток тепла.
function Qflux = externalHeatFlux(region,~)
[phi,theta,~] = cart2sph(region.x,region.y,region.z);
theta = pi/2 - theta; % transform to 0 <= theta <= pi
ids = phi > 0;
Qflux = zeros(size(region.x));
Qflux(ids) = theta(ids).^2.*(pi - theta(ids)).^2.*phi(ids).^2.*(pi - phi(ids)).^2;
end
Постройте поток на поверхности.
[phi,theta,r] = meshgrid(linspace(0,2*pi),linspace(-pi/2,pi/2),6); [x,y,z] = sph2cart(phi,theta,r); region.x = x; region.y = y; region.z = z; flux = externalHeatFlux(region,[]); figure surf(x,y,z,flux,'LineStyle','none') axis equal view(130,10) colorbar xlabel 'x' ylabel 'y' zlabel 'z' title('External Flux')
Включайте это граничное условие в модель.
thermalBC(thermalmodel,'Face',4,'HeatFlux',@externalHeatFlux,'Vectorized','on');
Задайте начальную температуру, чтобы быть нулем во всех точках.
thermalIC(thermalmodel,0);
Задайте вектор такта и решите переходную тепловую задачу.
tlist = [0,2,5:5:50]; R = solve(thermalmodel,tlist);
Чтобы строить контуры в несколько раз, с уровнями контура, являющимися тем же самым для всех графиков, определяют область значений температур в решении. Минимальная температура является нулем, потому что это - граничное условие на внутренней поверхности сферы.
Tmin = 0;
Найдите максимальную температуру из итогового решения такта.
Tmax = max(R.Temperature(:,end));
Постройте контуры в области значений Tmin
к Tmax
в те времена в tlist
.
h = figure; for i = 1:numel(tlist) pdeplot3D(thermalmodel,'ColorMapData',R.Temperature(:,i)) caxis([Tmin,Tmax]) view(130,10) title(['Temperature at Time ' num2str(tlist(i))]); M(i) = getframe; end
Чтобы посмотреть фильм контуров при выполнении этого примера на компьютере, выполните следующую линию:
movie(M,2)
Визуализируйте температурные контуры на поперечном сечении. Во-первых, задайте прямоугольную сетку точек на плоскость, где .
[YG,ZG] = meshgrid(linspace(-6,6,100),linspace(-6,6,100)); XG = zeros(size(YG));
Интерполируйте температуру в узлах решетки. Выполните интерполяцию для нескольких временных шагов, чтобы наблюдать эволюцию температурных контуров.
tIndex = [2,3,5,7,9,11]; varNames = {'Time_index','Time_step'}; index_step = table(tIndex.',tlist(tIndex).','VariableNames',varNames); disp(index_step);
Time_index Time_step __________ _________ 2 2 3 5 5 15 7 25 9 35 11 45
TG = interpolateTemperature(R,XG,YG,ZG,tIndex);
Задайте геометрические сферические слои на поперечном сечении.
t = linspace(0,2*pi); ylayer1 = cos(t); zlayer1 = sin(t); ylayer2 = 2*cos(t); zlayer2 = 2*sin(t); ylayer3 = 4*cos(t); zlayer3 = 4*sin(t); ylayer4 = 6*cos(t); zlayer4 = 6*sin(t);
Постройте контуры в области значений Tmin
к Tmax
для временных шагов, соответствующих индексам времени tIndex
.
figure('Position',[10,10,1000,550]); for i = 1:numel(tIndex) subplot(2,3,i) contour(YG,ZG,reshape(TG(:,i),size(YG)),'ShowText','on') colorbar title(['Temperature at Time ' num2str(tlist(tIndex(i)))]); hold on caxis([Tmin,Tmax]) axis equal % Plot boundaries of spherical layers for reference. plot(ylayer1,zlayer1,'k','LineWidth',1.5) plot(ylayer2,zlayer2,'k','LineWidth',1.5) plot(ylayer3,zlayer3,'k','LineWidth',1.5) plot(ylayer4,zlayer4,'k','LineWidth',1.5) end
[1] Сингх, Suneet, P. K. Джайн и Rizwan-uddin. "Аналитическое Решение для 3D, Неустойчивой Проводимости Тепла в Многоуровневой Сфере". ASME. J. Теплопередача. 138 (10), 2016, стр 101301-101301-11.