%FFFF_MET_CONV the code %Copyright 2019 Jie ZHANG (ORCID: 0000-0002-8023-8144) % estimate the storage and time needed before calculation with fine meshes % play any PPT can prevent your PC from going sleep clear; tic, global nl nw nh %------------ settings --------------------------------------------------- filename2use = 'test'; % alwayse change a filename, otherwise you will over-write previous results %filename2use = 'conv_v10h60Tn203Tp571ta_correction_FULL'; for oo=0 % this for-loop is used to claspe the code nl = 81; % element number in the L (length) direction. Do not set too high, recommendation: nl<200 lengthh = 18; width = 0.8; height = 12; % bulk dimension of the geometry, unit mm dx = lengthh/nl; % element length, unit mm dy = .4; % nozzle diameter, unit mm dz = .3; % layer thickness, unit mm vnozzle = 10; % speed of nozzle movement, unit mm/s dt = dx/vnozzle; % time step nw = round(width/dy); % element number in the W (width) direction nh = round(height/dz); % element number in the H (height) direction % please make sure min{nl, nw, nh}>=2 (for now). nw can be set to 1, e.g. to model printing single wall structure, but a different algorithm will be needed regarding the BC treatment. ne = nl*nw*nh; % total amounts of elements totaltime = dt*ne; % total time for THE FFF process (not the time for the calculation) tn = 203; % actual nozzle temperature (initial condition) tp = 57.1; % actual plate temperature ta = 21.2; % air temperature (far environment temperature) TaCal = @(Vd) 30.25*exp(-0.1059*Vd) + 21.2; % air temperature above the plate as a function of vertical distance (Vd) (near environment temperature) datacp = load ('placp.txt'); tem_Cp = datacp(:,1)-273.15; Cp = datacp(:,2)/72.06; % specific heat capacity as a function od temperature, unit: J/g/K tem_Rho = [38.79 50.13 62.41 75.25 88.43 102.4 116.7 132 147.7 163.3 179.3 195.3 211.2 230.4]; Rho = 0.001./[.8052 .8108 .818 .8259 .8369 .8526 .8638 .8753 .8879 .9005 .9142 .9279 .9435 .9601]; % temperature-dependent density, unit: g/mm3 lambda = 1.95e-4; % temperature-dependent value can be formulized, unit: W/mm/oC [the unit of lambda in THE paper, is wrong, it should be W/mm/K. But is is purely a typo ] h_conv = 60e-6; % convective heat transfer coefficient, unit: W/mm^2/K (h = 60e-6 W/mm^2/K corresponds to forced convection; for natural convection h = 10e-6 W/mm^2/K can be taken. But more dedicated mesaurement on these coefficients are needed for more accuracy description) h_conv_end = 10e-6; % convective heat transfer coefficient when FFF is finished, unit: W/mm^2/K (h = 10e-6 W/mm^2/K corresponds to natural convection) tem_Nu = linspace(20,215); Nu = lambda./(interp1(tem_Rho,Rho,tem_Nu,'pchip') .* interp1(tem_Cp,Cp,tem_Nu,'pchip')); % Nu denotes the thermal diffusivity, unit: mm^2/s. It is taken as a function of temperature. But in fact, its dependence on temperature is not strong for PLA. age = [1:-1:-ne+2]; % ages of element at t= dt. It is designed to judge if en element has been deposited or not. See the reference T =[tn; nan + zeros(ne-1,1)] ; % temperature of ALL elements at t= dt. NaN means the element has not yet been deposited, thus its temperature is 'Not A Nmber' end %% -------------the actual calculations-------------------------------------- for n = 2:2*ne tem = T(:,end); % take temperature fields everywhere at the previous time spot if n <= ne % indicating FFF is still in progress for m = 1:n-1 % from the first element to the (n-1)th elements deposited so far t_itself = tem(m); % find the temperature of element m. [i,j,k] = n2ijk(m); % find the spatial indexes of element m. (find where m is) ta = TaCal(k*dz); % update the air temperature around m % prepare the temperature of NEIGHBOORs of elements m for oo=1 % down and right % down ---------------- if k-1<1 t_down = tp; else id = ijk2n(i,j,k-1); t_down = tem(id); end %----------------------- % right ---------------- if j-1<1 % j=1 on the right boundary % first type BC: t_right = ta; if j+1<=nw % neighbor within the region id = ijk2n(i,j+1,k); % ID of its left neighbor if age(id) >0 % use O(h^2) t_right = tem(id) - 2*dy*h_conv/lambda * (t_itself -ta) ; else % use O(h) t_right = t_itself - dy*h_conv/lambda *( t_itself -ta); end end else id = ijk2n(i,j-1,k); t_right = tem(id); end %----------------------- end for oo=2 % up % up ------------------------ if k+1> nh % k=nk, the last layer, make sure the calculate lasts for at least 2 layers, use O(h^2) id = ijk2n(i,j,k-1); t_up = tem(id)-2*dz*h_conv/lambda*(t_itself-ta); else id = ijk2n(i,j,k+1); if age(id) >0 t_up = tem(id); else if k==1 % use O(h) t_up = t_itself - dz*h_conv/lambda*(t_itself-ta); else %use O(h^2) id = ijk2n(i,j,k-1); t_up = tem(id) - 2*dz*h_conv/lambda*(t_itself-ta); end end end end for oo=3 % left % left ------------------------ if j+1> nw % there is always element on the right, use O(h^2) id = ijk2n(i,j-1,k); t_left = tem(id)-2*dy*h_conv/lambda*(t_itself-ta); else id = ijk2n(i,j+1,k); if age(id) >0 t_left = tem(id); else if j==1 % use O(h) t_left = t_itself - dy*h_conv/lambda*(t_itself-ta); else %use O(h^2) id = ijk2n(i,j-1,k); t_up = tem(id) - 2*dy*h_conv/lambda*(t_itself-ta); end end end end for oo=4 % front % front ------------------------ if i+1> nl id = ijk2n(i-1,j,k); if age(id) >0 % use O(h^2) t_front = tem(id)-2*dx*h_conv/lambda*(t_itself-ta); else % use O(h) t_front = t_itself - dx*h_conv/lambda*(t_itself-ta); end else id = ijk2n(i+1,j,k); if age(id) >0 t_front = tem(id); else if i==1 % use O(h) t_front = t_itself - dx*h_conv/lambda*(t_itself-ta); else %use O(h^2) id = ijk2n(i-1,j,k); t_front = tem(id) - 2*dx*h_conv/lambda*(t_itself-ta); end end end end for oo=5 % back % back ------------------------ if i-1<1 id = ijk2n(i+1,j,k); if age(id) >0 % use O(h^2) t_back = tem(id)-2*dx*h_conv/lambda*(t_itself-ta); else % use O(h) t_back = t_itself - dx*h_conv/lambda*(t_itself-ta); end else id = ijk2n(i-1,j,k); if age(id) >0 t_back = tem(id); else if i==nl % use O(h) t_back = t_itself - dx*h_conv/lambda*(t_itself-ta); else %use O(h^2) id = ijk2n(i+1,j,k); t_back = tem(id) - 2*dx*h_conv/lambda*(t_itself-ta); end end end end tem_previous = tem(m); a = interp1(tem_Nu,Nu,tem_previous,'pchip'); % update the diffusivity % this interp1 function takes a lot of time!!! fx = a * dt/(dx*dx); fy = a * dt/(dy*dy); fz = a * dt/(dz*dz); % intermediates. See THE paper tem(m) = fx*(t_front + t_back) + fy*(t_left + t_right) + ... fz*(t_up + t_down) + (1-2*(fx + fy +fz)) * tem(m); % the temperature of element at this moment based on the temperature of all its neighbors and itself at previous moment end tem(n) = tn; % there is a new element deposition, whose temperature is determined by the initial condition. age = age +1; % update the age of ALL elments T = [T,tem]; % store the temperature of ALL elementst at the current moments else % when FFF is finished, and cooling continues, with a LOWer convective coefficient h_conv = h_conv_end; % W/mm^2/K natural convection for m = 1:ne % for ALL elements t_itself = tem(m); [i,j,k] = n2ijk(m); ta = TaCal(k*dz); % prepare the temperature of NEIGHBOORs for oo=1 % down and right % down ---------------- if k-1<1 t_down = tp; else id = ijk2n(i,j,k-1); t_down = tem(id); end %----------------------- % right ---------------- if j-1<1 % j=1 on the right boundary % first type BC: t_right = ta; if j+1<=nw % neighbor within the region id = ijk2n(i,j+1,k); % ID of its left neighbor if age(id) >0 % use O(h^2) t_right = tem(id) - 2*dy*h_conv/lambda * (t_itself -ta) ; else % use O(h) t_right = t_itself - dy*h_conv/lambda *( t_itself -ta); end end else id = ijk2n(i,j-1,k); t_right = tem(id); end %----------------------- end for oo=2 % up % up ------------------------ if k+1> nh % k=nk, the last layer, make sure the calculate lasts for at least 2 layers, use O(h^2) id = ijk2n(i,j,k-1); t_up = tem(id)-2*dz*h_conv/lambda*(t_itself-ta); else id = ijk2n(i,j,k+1); if age(id) >0 t_up = tem(id); else if k==1 % use O(h) t_up = t_itself - dz*h_conv/lambda*(t_itself-ta); else %use O(h^2) id = ijk2n(i,j,k-1); t_up = tem(id) - 2*dz*h_conv/lambda*(t_itself-ta); end end end end for oo=3 % left % left ------------------------ if j+1> nw % there is always element on the right, use O(h^2) id = ijk2n(i,j-1,k); t_left = tem(id)-2*dy*h_conv/lambda*(t_itself-ta); else id = ijk2n(i,j+1,k); if age(id) >0 t_left = tem(id); else if j==1 % use O(h) t_left = t_itself - dy*h_conv/lambda*(t_itself-ta); else %use O(h^2) id = ijk2n(i,j-1,k); t_up = tem(id) - 2*dy*h_conv/lambda*(t_itself-ta); end end end end for oo=4 % front % front ------------------------ if i+1> nl id = ijk2n(i-1,j,k); if age(id) >0 % use O(h^2) t_front = tem(id)-2*dx*h_conv/lambda*(t_itself-ta); else % use O(h) t_front = t_itself - dx*h_conv/lambda*(t_itself-ta); end else id = ijk2n(i+1,j,k); if age(id) >0 t_front = tem(id); else if i==1 % use O(h) t_front = t_itself - dx*h_conv/lambda*(t_itself-ta); else %use O(h^2) id = ijk2n(i-1,j,k); t_front = tem(id) - 2*dx*h_conv/lambda*(t_itself-ta); end end end end for oo=5 % back % back ------------------------ if i-1<1 id = ijk2n(i+1,j,k); if age(id) >0 % use O(h^2) t_back = tem(id)-2*dx*h_conv/lambda*(t_itself-ta); else % use O(h) t_back = t_itself - dx*h_conv/lambda*(t_itself-ta); end else id = ijk2n(i-1,j,k); if age(id) >0 t_back = tem(id); else if i==nl % use O(h) t_back = t_itself - dx*h_conv/lambda*(t_itself-ta); else %use O(h^2) id = ijk2n(i+1,j,k); t_back = tem(id) - 2*dx*h_conv/lambda*(t_itself-ta); end end end end tem_previous = tem(m); a = interp1(tem_Nu,Nu,tem_previous,'pchip'); % diffusivity fx = a * dt/(dx*dx); fy = a * dt/(dy*dy); fz = a * dt/(dz*dz); tem(m) = fx*(t_front + t_back) + fy*(t_left + t_right) + ... fz*(t_up + t_down) + (1-2*(fx + fy +fz)) * tem(m); end T = [T,tem]; end end save(filename2use) % save matrix T to a differnt data type to reduce the storage if necessary endtime = toc % notification ms2 = ['Running time ', num2str(endtime),'s']; msgh = msgbox({'Done!',ms2},'Over'); %%- voice notification, use only when long running time is required % load handel % sound(y,Fs)