Теплопроводность в многооменной геометрии с неоднородным тепловым потоком

В этом примере показано, как выполнить анализ 3-D переходной теплопроводности полой сферы, изготовленной из трех различных слоев материала.

Сфера подвергается неравномерному внешнему тепловому потоку.

Физические свойства и геометрия этой задачи описаны в Singh, Jain и Rizwan-uddin (см. Ссылка), который также имеет аналитическое решение для этой задачи. Внутренняя поверхность сферы всегда имеет нулевую температуру. Внешняя полусфера с положительным $$$$значением y имеет неравномерный тепловой поток, определяемый

$$\lambda$$ и $$\lambda$$ - азимутальные и углы возвышения точек в сфере. Изначально температура во всех точках сферы равна нулю.

Создание тепловой модели для нестационарного теплового анализа.

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')

Создайте сетку для геометрии. Выберите размер сетки, который является достаточно грубым, чтобы ускорить решение, но достаточно тонким, чтобы представить геометрию достаточно точно.

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)

Визуализация температурных контуров на поперечном сечении. Сначала определите прямоугольную сетку точек на $$$$плоскости y-z, где $$x\mathrm{=}$$ 0.

[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