Этот пример показывает, как решить уравнение тепла с зависящей от температуры теплопроводностью. Пример показывает идеализированный тепловой анализ прямоугольного блока с прямоугольной полостью в центре.
Дифференциальное уравнение с частными производными для переходной теплопередачи является:
где - температура, - плотность материала, является удельной теплотой, и - теплопроводность. - тепло, генерируемое внутри тела, которое в этом примере равняется нулю.
Создайте статическую тепловую модель.
thermalmodelS = createpde('thermal','steadystate');
Создайте 2-D геометрию путем рисования на один прямоугольник размера блока и на второй прямоугольник размера паза.
r1 = [3 4 -.5 .5 .5 -.5 -.8 -.8 .8 .8]; r2 = [3 4 -.05 .05 .05 -.05 -.4 -.4 .4 .4]; gdm = [r1; r2]';
Вычесть второй прямоугольник из первого, чтобы создать блок с пазом.
g = decsg(gdm,'R1-R2',['R1'; 'R2']');
Преобразуйте формат decsg в геометрический объект. Включите геометрию в модель.
geometryFromEdges(thermalmodelS,g);
Постройте график геометрии с отображенными метками ребер. Метки ребра будут использоваться ниже в функции определения граничных условий.
figure pdegplot(thermalmodelS,'EdgeLabels','on'); axis([-.9 .9 -.9 .9]); title 'Block Geometry With Edge Labels Displayed'
Установите температуру на левый край равной 100 степеней. На правом ребре есть предписанный тепловой поток из блока. Верхние и нижние ребра и ребра внутри полости все изолированы, то есть тепло не передается через эти ребра.
thermalBC(thermalmodelS,'Edge',6,'Temperature',100); thermalBC(thermalmodelS,'Edge',1,'HeatFlux',-10);
Задайте теплопроводность материала. Сначала рассмотрим постоянную теплопроводность, например, равную единице. Позже рассмотрим случай, когда теплопроводность является функцией от температуры.
thermalProperties(thermalmodelS,'ThermalConductivity',1);Создайте mesh с элементами не более 0,2.
generateMesh(thermalmodelS,'Hmax',0.2); figure pdeplot(thermalmodelS); axis equal title 'Block With Finite Element Mesh Displayed'
Вычислите установившееся решение.
R = solve(thermalmodelS); T = R.Temperature; figure pdeplot(thermalmodelS,'XYData',T,'Contour','on','ColorMap','hot'); axis equal title 'Temperature, Steady State Solution'
Создайте переходную тепловую модель и включите геометрию.
thermalmodelT = createpde('thermal','transient'); r1 = [3 4 -.5 .5 .5 -.5 -.8 -.8 .8 .8]; r2 = [3 4 -.05 .05 .05 -.05 -.4 -.4 .4 .4]; gdm = [r1; r2]'; g = decsg(gdm,'R1-R2',['R1'; 'R2']'); geometryFromEdges(thermalmodelT,g);
Задайте теплопроводность, массовую плотность и удельное тепло материала.
thermalProperties(thermalmodelT,'ThermalConductivity',1,... 'MassDensity',1,... 'SpecificHeat',1);
Задайте граничные условия. В переходных случаях температура на левом краю равна нулю в момент времени = 0 и наклоняется на 100 степени в течение 5 секунд. Вы можете найти функцию helper transientBCHeatedBlock под matlab/R20XXx/examples/pde/main.
thermalBC(thermalmodelT,'Edge',6,'Temperature',@transientBCHeatedBlock);
На правом ребре есть предписанный тепловой поток из блока.
thermalBC(thermalmodelT,'Edge',1,'HeatFlux',-10);
Верхние и нижние ребра, а также ребра внутри полости все изолированы, то есть тепло не передается через эти ребра.
Создайте mesh с элементами не более 0,2.
msh = generateMesh(thermalmodelT,'Hmax',0.2); figure pdeplot(thermalmodelT); axis equal title 'Block With Finite Element Mesh Displayed'
Вычислите переходное решение. Выполните временной анализ от нуля до пяти секунд. Тулбокс сохраняет решение каждые 1 секунду, чтобы можно было создать графики результатов как функций времени.
tlist = 0:.1:5; thermalIC(thermalmodelT,0); R = solve(thermalmodelT,tlist); T = R.Temperature;
Два графика полезны для понимания результатов этого переходного анализа. Первый - график температуры в окончательное время. Второй - график температуры в определенной точке блока, в этом случае около центра правого ребра, как функции времени. Чтобы идентифицировать узел вблизи центра правого ребра, удобно задать эту короткую служебную функцию.
getClosestNode = @(p,x,y) min((p(1,:) - x).^2 + (p(2,:) - y).^2);
Вызовите эту функцию, чтобы получить узел около центра правого ребра.
[~,nid] = getClosestNode( msh.Nodes, .5, 0 );
Эти два графика показаны один за другим на рисунке ниже. Распределение температуры в это время очень похоже на распределение, полученное из установившегося решения выше. На правом ребре, в разы меньше, чем примерно на половину секунды, температура меньше нуля. Это связано с тем, что тепло выходит из блока быстрее, чем оно приходит с левого края. В случаях, превышающих примерно три секунды, температура по существу достигла установившегося состояния.
h = figure; h.Position = [1 1 2 1].*h.Position; subplot(1,2,1); axis equal pdeplot(thermalmodelT,'XYData',T(:,end),'Contour','on','ColorMap','hot'); axis equal title 'Temperature, Final Time, Transient Solution' subplot(1,2,2); axis equal plot(tlist, T(nid,:)); grid on title 'Temperature at Right Edge as a Function of Time'; xlabel 'Time, seconds' ylabel 'Temperature, degrees-Celsius'
Часто, когда свойства материала являются функциями зависимых переменных. Например, предположим, что теплопроводность является простой линейной функцией от температуры:
k = @(~,state) 0.3+0.003*state.u;
В этом случае переменная u является температурой. В данном примере примите, что плотность и удельное тепло не являются функциями температуры.
thermalProperties(thermalmodelS,'ThermalConductivity',k);Вычислите установившееся решение. По сравнению со случаем постоянной проводимости температура на правом ребре ниже. Это связано с более низкой проводимостью в областях с более низкой температурой.
R = solve(thermalmodelS); T = R.Temperature; figure pdeplot(thermalmodelS,'XYData',T,'Contour','on','ColorMap','hot'); axis equal title 'Temperature, Steady State Solution'
Теперь выполните переходный анализ с температурно-зависимой проводимостью.
thermalProperties(thermalmodelT,'ThermalConductivity',k,... 'MassDensity',1,... 'SpecificHeat',1);
Используйте тот же интервал времени tlist = 0:.1:5 что касается линейного случая.
thermalIC(thermalmodelT,0); R = solve(thermalmodelT,tlist); T = R.Temperature;
Постройте график температуры на конечном временном шаге и температуры на правом ребре как функцию времени. График температуры на заключительном временном шаге лишь немного отличается от сопоставимого графика от линейного анализа: температура на правом ребре немного ниже линейного случая. График температуры как функции времени значительно отличается от линейного случая. Из-за меньшей проводимости при более низких температурах, тепло занимает больше времени, чтобы достигнуть правого ребра блока. В линейном случае температура по существу постоянна на уровне около трех секунд, но для этого нелинейного случая кривая температуры только начинает сглаживаться через пять секунд.
h = figure; h.Position = [1 1 2 1].*h.Position; subplot(1,2,1); axis equal pdeplot(thermalmodelT,'XYData',T(:,end),'Contour','on','ColorMap','hot'); axis equal title 'Temperature, Final Time, Transient Solution' subplot(1,2,2); axis equal plot(tlist(1:size(T,2)), T(nid,:)); grid on title 'Temperature at Right Edge as a Function of Time (Nonlinear)'; xlabel 'Time, seconds' ylabel 'Temperature, degrees-Celsius'
