This preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
Unformatted text preview: 4/24/10 EP 471 — Engineering Problem SolvingII
Exercise 19: Solution of Parabolic PDEs Using Matlab’s pdepe Utility You might recall from Exercise 16‘1that we classiﬁed second order PDEs as elliptic,
parabolic or hyperbolic depending on the relative size of the coefﬁcients in: 2' 2 2
aa—ZH} a Z +c—a——+dg+eg+fz
6x2 8x6y 8y2 ~ 8x If we restrict ourselves to one spatial dimension, transient heat conduction problems fall
into this class with b = c = O and with the spatial variable y replaced with time If. You ’ may recall from that exercise that an equation is classiﬁed as parabolic if b2 = 4ac. This 
is the case with b = c = 0. ' The governing equation in transient heat conduction with no internal heat generation and
constant properties is:  s—,——=V2,T \‘ 05:1. (1). where T 18 the temperature, p is the density, cp is the heat capacity, and k 13 the thermal
conductivity of the solid. The spatial dependence ls embodied 1n the Laplacian operator,
and its speciﬁc form varies depending on the coordinate system used 1n the problem. Matlab has a utility for solving systems of parabolic partial differential equations of the
form: ' 624 at; _,,, 5 ,,, Eu , 5a
a a —_ —: — a 9 —— + a 9 —
C(xt u, ,ax) 6t x 6x06 f(xt u, 6x» S(xt u, 36x ) Here the “c” is a capacity, “f” is a ﬂux ands“ ” is a source term. In its simplest form, the
ﬂux might be deﬁned as f: 614/ 6x. In general heat conduction problems with properties that are not constant, it may be written as f: k(T)8T / 6x. Them“ ” mparameter is used to choose between Cartesian (m= 0), cylindrical (m— — 1) or spherical (m= 2) geometries.
When m takes 011 one of these three integer values, it becomes consistent with the
appropriate form of the Laplacian in Equation (1). (It should be noted that, in the case of '
cylindrical or spherical problems, it is also implied that these are symmetric problems, as
there is no dependence on angular variables.) Also, although our classiﬁcation of
equations (elliptic, parabolic, hyperbolic) 13 based on a linear form for the equation, the
Matlab form is intended to permit nonlinear dependencies in c, f and s. i ' Download trans_heat_cond_mat1ab.m from the eCOW2 web page along with its
associated functions: pde_’trans_heat.m, pde_trans_bc.m and pde_trans_ic.m. In this 4/24/10 problem we are beginning with a slab With temperature T o = 0 throughout the slab. No
heat source is present. The left boundary condition is that the temperature remains at
zero. The right boundary condition is that the temperature is raised linearly with time.
This problem also has an analytical solution that can be compared to the numerical solution. The key statements in the main script are: x = linspace(0,l,lOl); t = [O l 3 10 30 100 300 1000 3000];' m = 0; soln_pde = pdepe (m, 'pde_trans_heat', 'pde_trans_ic' , 'pde_trans#~bc' ,x,t) ; T he“ ” array speciﬁes the mesh density, “1‘” provides the times at which the solution is
to be computed “m” is zero because it’s a slab (Cartesian) geometry and the ﬁrst, second
and third functions mprovide the speciﬁc functional form of the PDE (pde_ trans _heat), the
initial conditions (pde_ trans _ic) and the boundary conditions (pde trans b_c). The
pdepe utility submits all this to a solver whose inner workings are invisible to us. The output, soln_pde, is athree dimensional array, soln _pde(1,j k), where the ﬁrst index
indicates time, the second space and the third is the component. (When the system
consists of a single entity, like temperature in the present problem, It is just 1.) According to Matlab documentation, “The coupling of the partial derivatives with respect to time is restricted to multiplication
by a diagonal matrix C(x,l‘,u, au / 626). The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic
equation and otherwise to a parabolic equation. There must be at least one parabolic _ equation” This seems odd in that a component equation with c = 0 Would eliminate time, producing
a PDE in just one spatial variable, which would make it.. .an ordinary differential
equation (which I guess we could also call “elliptic, ” but that seems like cheating.) It
looks to me like this utility is most useful for systems of parabolic PDEs. We should just leave it at that. For transient heat cenduction, the more general form of the energy balance pr‘oduces: Since p, 617 and k are all material properties, they are generally functions of, among other
things, temperature, so having the ability to specify these as functions of the variable
could be quite useﬁrl for realistic problems. When they are treated as constants, we pull
the k out of the differential operator and divide the whole equation by k, producing the l/a in front of the 8T / at term in Equation (1). . 4/24/10 In the problem you’ve downloaded, the source term is zero and the material properties
are all constants. The pde_trans_heat.m script consists entirelyof specifying the form of
c, f and s in Matlab’s default form (highlighted equation on page one): function [c,f,s] = pde_trans_heat(x,t,T,DTDx)
global alpha; c — l/alpha; f DTDX; s 0; We also have to specify initial and boundary conditions. In the present problem, we
begin with temperature equal to zero everywhere, so the initial condition function is just: function T_O = pde_trans_ic(x)
T_O = 0; Boundary conditions are trickier, just as they were in the bvp4 c solver. Matlab is again
expectingthe boundary conditions to be speciﬁed in homogeneous form: “At the,
boundary x = a (left) and x = [7 (right), the solution components satisfy a boundary
condition of the form: ‘ p(x,r,u)+q(x,t)f(x,t,u,:—u) =0
X The matrix q(x,t) is a diagonal matrix with elements that are either identically zero or I never zero. Note that boundary conditions are expressed in terms of the ﬂux f rather than
the gradient au / 6x. Also, of the two coefﬁcients, only p Can depend on u. The ﬁle
pde_trans_bc.m consists of: I ' function [pl,ql,pr,qr] = pde_trans_bc(xl,Tl,xr,Tr,t)
global T_dot; pl‘= Tl; ql = 0; . pr = Tr  T_dot*t; 2
qr = O; In this case we are specifying a temperature equal to zero at all times on the left boundary
and a linearly increasing temperature with time on the right boundary. Since neither of
these involve the ﬂux, values of g on both the left and the right Sides are set to zero. On
the left side, we simply set pl equal to Tl. Given the assumed form of the equation, this
is equivalent to setting the temperature on the left boundary to zero. On the right side,
we set pr equal to the difference between the temperature on the right side and the '
ramped temperature. Again, given the assumed form for the boundary condition, this is
equivalent to setting the temperature on the right boundary equal to a linearly increasing function of time. 4/24/10 For this problem, there is also an analytical solution of the form: i(—.1)"1;9mjn—Msin(m/z) x12 12 7: "=1 n 71' a T(x,z‘) = 12% + 2M" The rest of the script involves evaluating the temperature proﬁle at various times and
comparing the results to those generated by the pdepe utility. Exercise: Realistic modiﬁcations to this problem involve convectiVe heat transfer
conditions to a Surrounding medium and the inclusion of a heat generation source.
Suppose we begin from a uniform temperature of 100 °C and have a constant heat source
' of 10000 W/m3. Suppose further that the surrounding medium is at 20 °C and heat
transfer coefﬁcients on the left and right boundaries are 1000 and 2000 W/mZK
respectively. Modify the scripts you’ve downloaded to incorporate these changes.
Although there isn’t an analytical solution for the entire transient, the solution should
approach an‘asymptotic steadystate for which there is an analytical solution. ...
View
Full Document
 Fall '09

Click to edit the document details