www.gusucode.com > pde 案例源码 matlab代码程序 > pde/heatedBlockWithSlot.m
%% Solving a Heat Transfer Problem With Temperature-Dependent Properties % This example shows an idealized thermal analysis of a rectangular, metal % block with a rectangular cavity in the center. One of the purposes of this % example is to show how temperature-dependent thermal conductivity can % be modeled in PDE Toolbox. % % The partial differential equation for transient conduction heat transfer % is: % % $$ % \rho C_p \frac{\partial T}{\partial t} -\nabla \cdot (k \nabla T) = f % $$ % % where $T$ is the temperature, $\rho$ is the material density, $C_p$ % is the specific heat, and $k$ is the thermal conductivity. $f$ is the % heat generated inside the body which is zero in this example. % % Copyright 2015-2016 The MathWorks, Inc. %% Create a PDE Model with a single dependent variable numberOfPDE = 1; pdem = createpde(numberOfPDE); %% Geometry % To construct the geometry we create one rectangle the size of the block % and a second rectangle the size of the slot. 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]'; %% % Subtract the second rectangle from the first to create the block with a % slot. g = decsg(gdm,'R1-R2',['R1'; 'R2']'); %% Convert the decsg format into a geometry object % The geometry is appended to the PDEModel geometryFromEdges(pdem,g); %% % Plot the geometry with edge labels displayed. The edge labels will be % used below in the function for defining boundary conditions. figure pdegplot(pdem,'EdgeLabels','on'); axis([-.9 .9 -.9 .9]); title 'Block Geometry With Edge Labels Displayed' %% PDE Coefficients and Boundary Conditions % Boundary conditions are defined below. For the steady state cases, we set % the temperature on the left edge to 100 degrees. On the right edge, there % is a prescribed heat flux out of the block. The top and bottom edges as % well as the edges inside the cavity are all insulated, that is, no heat % is transferred across these edges. uRight = applyBoundaryCondition(pdem,'neumann','Edge',1,'g',-10); uLeft = applyBoundaryCondition(pdem,'dirichlet','Edge',6,'u',100); %% Coefficient Definition % We will first consider the case where the material properties are simple % scalars, all equal one. We will later consider a case where the % thermal conductivity is a function of temperature. rho = 1; % density Cp = 1; % specific heat k = 1; % thermal conductivity c = k; a = 0; f = 0; d = rho*Cp; specifyCoefficients(pdem,'m',0,'d',0,'c',c,'a',a,'f',f); %% Mesh % Call |generateMesh| to create a mesh with elements no larger than % about 0.2. hmax = .2; % element size msh=generateMesh(pdem,'Hmax',hmax); figure pdeplot(pdem); axis equal title 'Block With Finite Element Mesh Displayed' %% Steady State Solution % First calculate the steady-state solution. % We want to compare this with the solution obtained below from a % transient analysis. R = solvepde(pdem); u = R.NodalSolution; figure pdeplot(pdem,'XYData',u,'Contour','on','ColorMap','hot'); axis equal title 'Temperature, Steady State Solution' %% PDE Coefficients and Boundary Conditions % Boundary conditions are defined below. % In the transient cases, the temperature % on the left edge is zero at time=0 and ramps to 100 degrees over .5 % seconds. On the right edge, there is a prescribed heat flux out of the % block. The top and bottom edges as well as the edges inside the cavity % are all insulated, i.e. no heat is transferred across these edges. % applyBoundaryCondition(pdem,'dirichlet','Edge',6,'u',@transientBCHeatedBlock); % Specify coefficients with non-zero d-coefficient to define the % time-dependent PDE. % specifyCoefficients(pdem,'m',0,'d',d,'c',c,'a',a,'f',f); %% Transient Solution % Perform a transient analysis from zero to five seconds. % The solution will be saved every .1 seconds so that plots of % the results as functions of time can be created. tlist = 0:.1:5; setInitialConditions(pdem, 0); R = solvepde(pdem,tlist); u = R.NodalSolution; %% % Two plots are useful in understanding the results from this transient % analysis. The first is a plot of the temperature at the final time. % The second is a plot of the temperature at a specific point in the block, % in this case near the center of the right edge, as a function of time. % To identify a node near the center of the right edge, it is convenient % to define this short utility function. % getClosestNode = @(p,x,y) min((p(1,:) - x).^2 + (p(2,:) - y).^2); % Call this function to get a node near the center of the right edge. [~,nid] = getClosestNode( msh.Nodes, .5, 0 ); %% % The two plots are shown side-by-side in the figure below. h = figure; h.Position = [1 1 2 1].*h.Position; subplot(1,2,1); axis equal pdeplot(pdem,'XYData',u(:,end),'Contour','on','ColorMap','hot'); axis equal title 'Temperature, Final Time, Transient Solution' subplot(1,2,2); axis equal plot(tlist, u(nid,:)); grid on title 'Temperature at Right Edge as a Function of Time'; xlabel 'Time, seconds' ylabel 'Temperature, degrees-Celsius' %% % As can be seen, the temperature distribution at this time is very similar % to that obtained from the steady-state solution above. At the right edge, % for times less than about one-half second, the temperature is less than % zero. This is because heat is leaving the block faster than it is arriving % from the left edge. At times greater than about three seconds, the % temperature has essentially reached steady-state. %% Steady State Solution (Temperature-Dependent Thermal Conductivity) % It is not uncommon for material properties to be functions of the % dependent variables. In this case we will assume that the thermal % conductivity is a simple linear function of temperature. The most convenient way % to define this function is by using a string expression: c = @(loc,state) 0.3+0.003*state.u; %% % In this case, the variable u is the temperature. % We will assume that the density and specific heat are not functions of % temperature but, if they were temperature-dependent, these coefficients % could also be defined as functions. specifyCoefficients(pdem,'m',0,'d',0,'c',c,'a',a,'f',f); % %Because the c-coefficient is a function of the dependent variable, we now %have a nonlinear PDE. Update the boundary condition on the left edge to be %a fixed temperature of 100 for steady state solution. Other edges would %continue to have boundary conditions as specified previously. applyBoundaryCondition(pdem,'dirichlet','Edge',6,'u',100); % %Providing a good initial guess to nonlinear solver yields faster %convergence. Hence, specify the solution values corresponding to final %time-step of transient solution, with constant c-coefficient, as initial %conditions to nonlinear solver. The time-dependent results are passed to %setInitialConditions and it selects the final time-step by default. If a %prior time-step were required it could be specified by an additional %argument that specified the time index. % setInitialConditions(pdem, R); R = solvepde(pdem); u = R.NodalSolution; figure pdeplot(pdem,'XYData',u,'Contour','on','ColorMap','hot'); axis equal title 'Temperature, Steady State Solution (Nonlinear)' %% % As can be seen, compared with the constant-conductivity case, the % temperature on the right-hand edge is lower. This is due to the lower % conductivity in regions with lower temperature. %% Transient Solution (Temperature-Dependent Thermal Conductivity) % We now perform a transient analysis with the temperature-dependent % conductivity. Update the boundary condition on the left edge to have % time-dependent boundary conditions. % applyBoundaryCondition(pdem,'dirichlet','Edge',6,'u',@transientBCHeatedBlock); % Specify coefficients with non-zero d-coefficient and temperature % dependent c-coefficients to define the nonlinear time-dependent PDE. % specifyCoefficients(pdem,'m',0,'d',d,'c',c,'a',a,'f',f); % Just as for the linear case, the timespan is from zero to five seconds setInitialConditions(pdem, 0); R = solvepde(pdem,tlist); u = R.NodalSolution; %% % Again we will plot the temperature at the final time and the temperature % at the right edge as a function of time. These two plots are shown % side-by-side in the figure below. % h = figure; h.Position = [1 1 2 1].*h.Position; subplot(1,2,1); axis equal pdeplot(pdem,'XYData',u(:,end),'Contour','on','ColorMap','hot'); axis equal title 'Temperature, Final Time, Transient Solution (Nonlinear)' subplot(1,2,2); axis equal plot(tlist(1:size(u,2)), u(nid,:)); grid on title 'Temperature at Right Edge as a Function of Time (Nonlinear)'; xlabel 'Time, seconds' ylabel 'Temperature, degrees-Celsius' %% % The plot of temperature at the final time is only slightly different % from the comparable plot from the linear analysis; temperature at the right % edge is slightly lower than the linear case. The plot of temperature as a % function of time is considerably different from the linear case. Because % of the lower conductivity at lower temperatures, the heat takes longer % to reach the right edge of the block. In the linear case, the temperature % is essentially constant at around three seconds but for this nonlinear % case, the temperature curve is just beginning to flatten at five seconds. displayEndOfDemoMessage(mfilename)