This tutorial is intended for the full version of the toolbox.
Introduction
This tutorial is continuation of Tutorial 7 on Steady State Heat Conduction in a Two-Component Material and it is adviced that you first analyse in details the previous tutorial, since we assume here that certain concepts and portions of code will be clear for you.
Problem Definition
We will simulate transient heat conduction in a two-component material composed of two distinct phases which differ in heat conductivity. Geometry with ids of Physical Lines and ids of Physical Surfaces are shown in figure below.
The problem is described by the following mathematical model
We assume the following boundary conditions:
- Line 12: Constant heat flux of 10 000 W/m^2
- Line 13: Insulated
- Line 14: Heat convection with htc = 20 W/(m2K) and ambient temperature of 293 K.
The initial condition is a temperature of 293 K everywhere.
Code
In this section we are going to show the Matlab script which solves the problem in QuickerSim CFD Toolbox. We will mainly describe the code in comments, since majority of steps should be clear after completing Tutorial 7.
clc; clear; % Import and display mesh [p,e,t] = importMeshGmsh('simpleRib.msh'); displayMesh2D(p,t); <img class="size-full wp-image-2396 aligncenter" src="http://quickersim.com/cfdtoolbox/wp-content/uploads/2017/07/QuickerSim-Heat-Transfer-Rib-mesh.png" alt="" width="560" height="420" /> % Define timestepping parameters nstep = 120; % number of steps dt = 240; % time step size % Define material parameters rho = 2000; % density cp = 1200; % heat capacity lambda1 = 50; % thermal conductivity of the surrounding material lambda2 = 400; % thermal conductivity of the rib material % Assign thermal conductivities to appropriate subdomains nelements = size(t,2); lambda = lambda2*ones(nelements,1); lambda(t(end,:)==19) = lambda1; % Init solution with constant value [T,convergence] = initScalarSolution(p,293); % Precompute mass matrix and diffusion matrix M = assembleMassMatrix2D(p,t,rho*cp/dt); D = assembleDiffusionMatrix2D(p,t,lambda);
We will shoot a movie of the temperature field evolution. To accompish that, we will use function from our toolbox. Additionally, we exploit Euler method for timestepping. This resuts in following iterative process for timestepping (k in upper index denotes time step index):
The discretized for with mass matrix M and diffusion matrix D reads:
which can be rearranged as
% We realize this algorithm in the timestepping loop below % Init movie movieObj = initMovie('ribHeatCondution.avi',16); % Begin timestepping for step = 1:nstep time = step*dt; % Compute transient part of the right hand side F = M*T; % Assemble problem matrix K = M + D; % Impose boundary conditions [K, F] = imposeScalarBoundaryCondition2D(p,e,K,F,12,'flux',10000); [K, F] = imposeScalarBoundaryCondition2D(p,e,K,F,14,'robin',20*293,20); % Solve for temperature field at the new time level T = K\F; % Report step number (uncomment code line below) % disp(step); % Display current temperature field and shoot movie frame figure(1); displaySolution2D(p,t,T,'Temperature [deg C]'); drawnow; writeVideo(movieObj,getframe(gcf)); end <img class="size-full wp-image-2397 aligncenter" src="http://quickersim.com/cfdtoolbox/wp-content/uploads/2017/07/QuickerSim-Heat-Transfer-Rib-Temperature-Map.png" alt="" width="560" height="420" /> % Close movie file close(movieObj);