MATLAB QUESTION:

Modify the code provided to solve the two-dimensional heat conduction equation with Dirichlet boundary conditions (temperature boundary conditions), and show that the implementation is correct.

function [T] = HeatConduction2DTemplate % Solves the 2D heat conduction equation % % d2T/dx2 + d2T/dy2 + egen/k = 1/alpha*(dT/dt) % % % Functions needed to setup: % Line 159: setBC(T,x,y,dx,dy,alpha,egen,k) % Line 189: computeResidual(T,dx,dy,alpha,egen,k) % Line 196: updateSolution(T,R,dt) % % Variables needed to setup % Line 62: Initialize the solution variable T % % Helper functions % plotTemperature(T,x,y,time) Plot the temperature % plotResidual(L2Resid) Plot the iterative residual % [T,S]=Tmms(x,y,alpha,egen,k) Evaluate the manufactured solution

% ========================================================================= % Define Physical properties % ========================================================================= k=171;%w/mk rho=2700;%kg/m3 c=869;%J/kgK alpha=k/(rho*c); %s/K egen=0;

% ========================================================================= % Geometry setup % ========================================================================= x0=-1;%m xL=1;%m nx=11; xvec=linspace(x0,xL,nx); dx=(xL-x0)/(nx-1);

y0=-1;%m yL=1;%m ny=11; yvec=linspace(y0,yL,ny); dy=(yL-y0)/(ny-1);

[x,y]=meshgrid(xvec,yvec);x=x';y=y';

% ========================================================================= % solver parameters % ========================================================================= maxIter=1; cfl=0.5; %max explicit cfl<=0.5 dt=cfl*min(dx,dy)^2/alpha; tolerance=1e-13; % Convergence tolerance, when to exit the time loop when steady state is reached nplot=100; %How often to update plots

% ========================================================================= % Initialize solution % ========================================================================= % Initialize the temperature variable( should be the same size as x and y ) %T=

% ========================================================================= % Time loop % ========================================================================= time=0; %initial time for n=1:maxIter % just a helper variable for short hand indexing i=2:nx-1; j=2:ny-1; % Set up boundary conditions T=setBC(T,x,y,dx,dy,alpha,egen,k); % Compute the iterative residual R=computeResidual(T,dx,dy,alpha,egen,k); % Compute arbitary MMS source term ------------------------------------ % This is for testing to show that you've coded everything correctly % Must be S=0 for actual simulations [~,S] = Tmms(x,y,alpha,egen,k); R=R-S(i,j); % --------------------------------------------------------------------- %---------------------------------------------------------------------- % Compute the L2 norm of the iterative residual. ---------------------- % Remember R=dT/dt so R=0 means you've reached steady state % The norm reduces the residual array to a single value L2Resid(n)=sqrt( sum(sum(R.^2))/numel(R) ); %Normalize Residuals -------------------------------------------------- % This just helps monitor convergence more consistently if n>5 L2Resid(n)=L2Resid(n)/L2Resid(5); end %---------------------------------------------------------------------- % Update Solution % Compute the Temperature at the new time step T = updateSolution(T,R,dt);

%---------------------------------------------------------------------- % Plot solution time=time+dt; if mod(n,nplot)==0 subplot(1,2,1) plotTemperature(T,x,y,time) subplot(1,2,2) plotResidual(L2Resid) pause(0.1) end %----------------------------------------------------------------------

%---------------------------------------------------------------------- % Check if steady state is reached or if solution is diverging if L2Resid(n)1e10 fprintf('Error! Divergence detected ') break end %---------------------------------------------------------------------- end

%---------------------------------------------------------------------- % Plot solution and iterative residual subplot(1,2,1) plotTemperature(T,x,y,time)

subplot(1,2,2) plotResidual(L2Resid) %----------------------------------------------------------------------

% Discretization error can be computed when an exact solution is known % and is the difference between the fully converged numerical solution and % the exact solution. The discretizatione error should decrease with step % size at a rate depending on the order of accuracy of the discretization % scheme. % % The observed order of accuracy is % p = log( DE1/DE2 ) / log( dx1/dx2 ) % as dx gets smaller DE should decrease so that p approaches the order of % accuracy of the discretization scheme

end

% Boundary conditions function Tout=setBC(T,x,y,dx,dy,alpha,egen,k)

Tout=T; nx=size(T,1); ny=size(T,2); % Boundary condition functions % Since boundary conditions are implemented in a 1D sense, it will % simplify the implementation by writing an anonymous function which is % a function of ds, T1, T2, T3, etc. For example, no heat flux boundary % condition is % % dTdx = 0 = (3T1 -4T2 + T3)/dx % % T1=(4T2-T3)/3 % % T_noHeatFlux = @(ds, T2, T3) (4*T2 - T3)/3; % % And for the left boundary you can call % % T(1,:) = T_noHeatFlux(dx, T(2,:), T(3,:)) %Add Boundary conditions

end

% Residual function R=computeResidual(T,dx,dy,alpha,egen,k)

%Add residual calculation end

% Time step function Tout = updateSolution(T,R,dt)

%Add time step

end

%% Helper functions that you don't need to modify % Setup to plot the temperature function plotTemperature(T,x,y,time)

nx=size(T,1); ny=size(T,2); % This extrapolates values to the corner to smooth the plot. The % corners aren't actually used for anything. T(1,1) = 2*T(2,2)-T(3,3); T(nx,1) = 2*T(nx-1,2)-T(nx-2,3); T(1,ny) = 2*T(2,ny-1)-T(3,ny-2); T(nx,ny) = 2*T(nx-1,ny-1)-T(nx-2,ny-2);

surf(x,y,T) shading interp view([0,0,90]);

title(['time=',num2str(time,2)]); colorbar end

% Setup to plot the iterative residual function plotResidual(L2Resid)

semilogy(L2Resid,'k-') xlabel('Iteration') ylabel('Residual') title('Iterative Convergence') end

% Exact solution for code verification function [T,S]=Tmms(x,y,alpha,egen,k) % Exact solution % This function was arbitarily chosen as an exact solution. It's not % technically a solution; however, it can be made one by adding a source % term S to the iterative residual. We can use this exact solution to % compute the discretization error in the numerical solution and make sure % we've programmed everything correctly. % % The current governing equation is % % alpha*(d2T/dx2 + d2T/dy2 + egen/k) = dT/dt % % The source term is computed by plugging the exact solution into the % governing equation % % alpha*(d2T/dx2 + d2T/dy2 + egen/k) - dT/dt = S % % If the chosen solution is actually an exact solution then S=0

Lx=1; Ly=1; cx=35; cy=27; c0=45; T = cx*sin( pi*x/Lx ) + cy*cos(pi*y/Ly) + c0; dTdx= cx*(pi/Lx)*cos(pi*x/Lx); d2Tdx2 = -cx*(pi/Lx)^2*sin(pi*x/Lx); dTdy = -cy*(pi/Ly)*sin(pi*y/Ly); d2Tdy2 = -cy*(pi/Ly)^2*cos(pi*y/Ly); S = alpha*d2Tdx2 + alpha*d2Tdy2 + alpha*egen/k;

end