function [] = Example_Transients () % revised 4/20/20 % Time stepping example of a transient conducting bar % with a mesh of three L2 linear line elements % PDE: K A d2T/dx2 - Q = rho cp A dT/dt % Matrix: [M]{T_dot} + [S]{T} = c_Q(t) % % Non-dimensional form [S*] + a [M*] = {c*(t)} % a = rho cp Le^2 / dt K A % See text Fig 15.3 for time history graphs %---------------------------------------------------- % EBC: V_0 1 2 3 4 V_L % .***(1)***.***(2)***.***(3)***. Q=0 % ICs: T_0 T_0 T_0 T_0 % Only T_2 and T_3 are free %------------------------MESH------------------------ I_C = 0 ; T_0 = I_C*[1; 1] ; % Initial Conditions V_0 = 10 ; V_L = 20 ; % Constant EBC values x=[0 1 2 3]' ; % x Coordinates n_steps = 30 ; % Number of time steps T = zeros (2, n_steps) ; % Allocate Free temperatures diag = 1 ;% turn on/off use only diagonal mass matrix ave = 0 ;% turn on/off use only averaged mass matrix a = 8.0 ; % *** Bigger for smaller Del_T K = [2 -1; -1 2] ; % Assembled conduction matrix M = [4 1; 1 4]*a/6 ; % Assembled capacity matrix Q = [0 ; 0] ; % Assembled source vector if ( diag == 1 ) ; % then use diagonal mass M = [5 0; 0 5]*a/6 ; % scaled diagonal of {M} end ; % if diagonal mass matrix if ( ave == 1 ) ; % then use averaged mass M = ([4 1; 1 4]*a/6 + [5 0; 0 5]*a/6)/2;% ave end ; % if averaged mass matrix S = K + M ; % New linear system at each time S_inv = inv(S) ; % Invert only once (constant Del_T) BC_1 = V_0*[-1;0]; BC_4 = V_L*[0;-1]; % c_EBC vectors c = M*T_0 - BC_1 - BC_4 ; % first c(t=Del_T) T(:,1) = S_inv * c ; % solve for first history result for n = 2:n_steps ; % time step loops ==> ==> ==> ==> c = M*T(:, n-1)-BC_1-BC_4 ; % update the RHS source T(:,n) = S_inv * c ; % solve for T at this new time end ; % for each step after the first <== <== <== <== % ----------------end time marching ----------------- % Restore EBC to graph mesh Show values through time Show = zeros (n_steps+1, 4) ; % allocate graph data Show(1,:) = I_C ; % place ICs in the first row, t=0 Show(:,1) = V_0 ; % put left BC in col 1 for all t Show(:,4) = V_L; ; % put right BC in col 4 for all t Show(2:n_steps+1, 2:3) = T' ; % insert computed T(t) % Scale the axes, label axes, set title xmin = 0; xmax = 3; % set length axis bounds ymin = min([T_0; V_0; V_L]) ; % set min temperature ymax = max([T_0; V_0; V_L]) ; % set max temperature clf; hold on ; grid on ; axis([xmin xmax ymin ymax]); xlabel('x'); ylabel('Temperature') ; % label axes title('FEA Transient Conduction Solution') ; % title % Add stready state temperatures for reference end_x=[0 3];end_t=[V_0 V_L];% steady state end values plot(end_x, end_t,'k--') ; % show dashed reference % Graph each t state approach steady state value for n = 1:n_steps+1; % loop over time --> --> --> --> plot (x, Show(n,:),'k-') ; % show element values plot (x, Show(n,:),'k*') ; % show node values end; % for inital and following steps <-- <-- <-- <-- plot (x, Show(1,:),'ko'); % repeat initial conditions ave=(V_0+V_L)/2 ; % text location text(0.5, ave,[' IC = ', num2str(I_C)]) ; % list IC print -dpng bar_transients; % save to a png plot file fprintf ('Created file bar_transients.png \n') hold off ; % all done % end Example_Transients ============================