clear; close all; % Set the values for parameters kair=0.024; % thermal conductivity of air pcpair=1230; % volumetric heat capacity of air thick=0.25; %thickness of the house walls L=5; % length sigma=5.67e-8; % stefan-boltzmann constant %Set some values for different Thermal Conductivities: kaluminum=237; kcopper=401; kdiamond=1000; kmanganese=7.81; kbrick=.6; % Set Values for Thermal Diffusivities daluminum=9.7e-5; dcopper=1.11e-4; ddiamond=2e-4; dmanganese=2.261e-6; dbrick=5.2e-7; % ask the user to choose a material for the walls of the house choice=input( 'What do you want to make your house out of? \n 1) Aluminum \n 2) Copper \n 3) Diamond \n 4) Manganese \n 5) Bricks \n Choice: ') if choice == 1 || choice ==2 || choice == 3 || choice == 4 || choice == 5 if choice == 1 kval=kaluminum; pcpedge=kval/daluminum; end if choice == 2 kval=kcopper; pcpedge=kval/dcopper; end if choice == 3 kval=kdiamond; pcpedge=kval/ddiamond; end if choice == 4 kval=kmanganese; pcpedge=kval/dmanganese; end if choice == 5 kval=kbrick; pcpedge=kval/dbrick; end else fprintf('Fine. Be that way. You just get bricks then. \n \n') kval=kbrick; pcpedge=kval/dbrick; end % build a cell-centered grid with N=50 cells % on the interval x=0 to x=L N=100; xmin=0; xmax=L; h=(xmax-xmin)/(N-1); x=xmin:h:xmax; x=x'; ymin=0; ymax=L; y=ymin:h:ymax; y=y'; [X,Y]=ndgrid(x,y); % Make grids of volumetric heat capacity values and thermal diffusivity values for the house pcp=zeros(length(x),length(y)); allk=zeros(length(x),length(y)); for l=1:N for m=1:N if x(l)> thick && x(l)< L-thick && y(m)> thick && y(m)< L-thick pcp(l,m)=pcpair; allk(l,m)=kair; else pcp(l,m)=pcpedge; allk(l,m)=kval; end end end % Make arrays of average k values plus and minus one half for each dimension kmx=allk*0; kpx=allk*0; kmx(2:N-1,1:N)=.5*(allk(2:N-1,1:N)+allk(1:N-2,1:N)); kmx(1,1:N)=allk(1,1:N); kmx(N,1:N)=allk(N,1:N); kpx(2:N-1,1:N)=.5*(allk(2:N-1,1:N)+allk(3:N,1:N)); kpx(1,1:N)=allk(1,1:N); kpx(N,1:N)=allk(N,1:N); kmy=allk*0; kpy=allk*0; kmy(1:N,2:N-1)=.5*(allk(1:N,2:N-1)+allk(1:N,1:N-2)); kmy(1:N,1)=allk(1:N,1); kmy(1:N,N)=allk(1:N,N); kpy(1:N,2:N-1)=.5*(allk(1:N,2:N-1)+allk(1:N,3:N)); kpy(1:N,1)=allk(1:N,1); kpy(1:N,N)=allk(1:N,N); % make a j variable for later j=2:N-1; k=2:N-1; % define the initial displacement and velocity vs. x z = 300*ones(length(x),length(y)); vz = zeros(length(x),length(y)); taulim=.25*h^2/(max(max(kpx./pcp))); fprintf(' tau limit: %g \n',taulim) % Choose a time step tau=input(' Enter the time step - ') % Choose how long to run and when to plot tfinal=input(' Enter tfinal - ') skip=input(' Enter # of steps to skip between plots (faster) - ') nsteps=tfinal/tau; % here is the loop that steps the solution along figure % open a new frame for the animation for n=1:nsteps time=n*tau; % compute the time % Use leapfrog and the boundary conditions to load znew with z % at the next time step using z znew=300*ones(length(x),length(y)); znew(2:N-1,2:N-1)=z(2:N-1,2:N-1)+(tau./(pcp(j,k).*(h^2))).*((kpx(j,k).*(z(j+1,k)-z(j,k))-kmx(j,k).*(z(j,k)-z(j-1,k)))+... (kpy(j,k).*(z(j,k+1)-z(j,k))-kmy(j,k).*(z(j,k)-z(j,k-1)))); znew(1,:)=znew(2,:)-(sigma*h/kval)*((z(2,:)+z(1,:))/2).^4; znew(N,:)=znew(N-1,:)-(sigma*h/kval)*((z(N,:)+z(N-1,:))/2).^4; znew(:,1)=znew(:,2)-(sigma*h/kval)*((z(:,2)+z(:,1))/2).^4; znew(:,N)=znew(:,N-1)-(sigma*h/kval)*((z(:,N)+z(:,N-1))/2).^4; z=znew; % make plots every skip time steps if mod(n,skip)==0 %contourf(z,10) surf(z) xlabel('x');ylabel('y'); title(['Time=' num2str(time)]) pause(.0001) end end