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

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

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

Физические свойства и геометрия этой проблемы описаны в Сингхе, джайне, и 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 для геометрии. Выберите размер 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