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

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

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

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

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

Создайте тепловую модель для переходного теплового анализа.

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 для геометрии. Выберите размер сетки, который является достаточно грубым, чтобы ускорить решение, но достаточно тонким, чтобы представлять геометрию достаточно точно.

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=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