% COMSOL Multiphysics Model M-file % Generated by COMSOL 3.5 snapshot (COMSOL 3.5.0.489, $Date: 2008/09/12 16:03:11 $) flclear fem % COMSOL version clear vrsn vrsn.name = 'COMSOL 3.5'; vrsn.ext = ' snapshot'; vrsn.major = 0; vrsn.build = 489; vrsn.rcs = '$Name: $'; vrsn.date = '$Date: 2008/09/12 16:03:11 $'; fem.version = vrsn; % Geometry g1=rect2('0.065','0.1','base','corner','pos',{'0','0'},'rot','0'); g2=rect2('0.0625','0.025','base','corner','pos',{'0','0.1'},'rot','0'); g3=rect2('0.11575','0.04','base','corner','pos',{'0','0.165'},'rot','0'); g4=rect2('0.11575','0.3675','base','corner','pos',{'0','0.205'},'rot','0'); g5=rect2('0.11575','0.4','base','corner','pos',{'0','0.5725'},'rot','0'); g6=rect2('0.11575','0.6','base','corner','pos',{'0','0.9725'},'rot','0'); g7=curve2([0,0],[0.125,0.165]); g8=curve2([0.11575,0.0625],[0.165,0.125]); g9=geomcoerce('solid',{g1,g2,g3,g4,g5,g6,g7,g8}); g10=mirror(g9,[0,0],[1,0]); [g11]=geomcopy({g10}); % Constants fem.const = {'T0','300[K]', ... 'T_in','1473[K]', ... 'v_cast','1.6[mm/s]', ... 'T_m','1356[K]', ... 'dT','20[K]', ... 'dH','205[kJ/kg]'}; % Geometry % Analyzed geometry clear s s.objs={g9}; s.name={'CO1'}; s.tags={'g9'}; fem.draw=struct('s',s); fem.geom=geomcsg(fem); % Initialize mesh fem.mesh=meshinit(fem, ... 'hauto',5); % Refine mesh fem.mesh=meshrefine(fem, ... 'mcase',0, ... 'rmethod','regular'); % (Default values are not included) % Application mode 1 clear appl appl.mode.class = 'HtNavierStokes'; appl.mode.type = 'axi'; appl.module = 'HT'; appl.shape = {'shlag(1,''ul_chns'')','shlag(1,''vl_chns'')','shlag(1,''p'')','shbub(2,''ub_chns'')','shbub(2,''vb_chns'')','shlag(2,''u'')','shlag(2,''v'')'}; appl.gporder = {4,2}; appl.cporder = {1,2}; appl.assignsuffix = '_chns'; clear prop prop.analysis='static'; prop.weakcompflow='On'; appl.prop = prop; clear bnd bnd.symtype = {'sym','sym','ax','sym','sym'}; bnd.v0 = {0,0,0,0,'v_cast'}; bnd.type = {'walltype','int','sym','inlet','outlet'}; bnd.outtype = {'p','p','p','p','uv'}; bnd.velType = {'U0in','U0in','U0in','U0in','u0'}; bnd.intype = {'uv','uv','uv','p','uv'}; bnd.ind = [3,4,3,2,3,2,3,2,3,2,3,2,3,2,5,1,1,1,1,1,5,5,5]; appl.bnd = bnd; clear equ equ.eta = {0.0434,1}; equ.rhofcnT = 1; equ.gporder = {{1;1;2}}; equ.name = {'default','Solid domain'}; equ.F_z = {'-Sz',0}; equ.Tflowtype = 'nonisoT'; equ.rho = 'rho_htgh'; equ.cporder = {1,{2;2;1}}; equ.F_r = {'-Sr',0}; equ.init = {{0;'v_cast';0;0;0;0;0;0;0},0}; equ.rhofcnp = 1; equ.shape = {[1;2;3;4;5],[6;7;3]}; equ.cdon = {0,1}; equ.sdon = {0,1}; equ.rhofcnTname = 'T'; equ.usage = {1,0}; equ.ind = [1,1,1,1,1,1,1]; appl.equ = equ; fem.appl{1} = appl; % Application mode 2 clear appl appl.mode.class = 'GeneralHeat'; appl.mode.type = 'axi'; appl.module = 'HT'; appl.shape = {'shlag(1,''J'')','shlag(1,''T'')','shlag(2,''T'')'}; appl.gporder = {2,4}; appl.cporder = {1,2}; appl.assignsuffix = '_htgh'; clear prop prop.analysis='static'; clear weakconstr weakconstr.value = 'off'; weakconstr.dim = {'lm8'}; prop.weakconstr = weakconstr; appl.prop = prop; clear bnd bnd.type = {'ax','T','cont','qc','q0','q','q','q'}; bnd.radType = {'none','none','none','none','none','none','none','surf2amb'}; bnd.epsilon = {0,0,0,0,0,0,0,0.8}; bnd.Tamb = {0,0,0,0,0,0,0,'T0'}; bnd.shape = 1; bnd.T0 = {273.15,'T_in',273.15,273.15,273.15,273.15,273.15, ... 273.15}; bnd.Tinf = {273.15,273.15,273.15,273.15,273.15,'T0','T0','T0'}; bnd.h = {0,0,0,0,0,25,800,10}; bnd.ind = [1,2,1,3,1,3,1,3,1,3,1,3,1,3,4,5,5,5,5,6,7,8,8]; appl.bnd = bnd; clear equ equ.eta = {'eta_chns',0}; equ.cdtypesystem = 'system'; equ.gporder = {1,2}; equ.pA = {'101325[Pa]+p',0}; equ.rho = {8500,8700}; equ.cporder = {1,2}; equ.init = {{'T_in';0},{273.15;0}}; equ.shape = {2,3}; equ.cdon = 1; equ.sdon = {1,0}; equ.C = {'Cp1+D*dH',385}; equ.matterstate = {'userdefined','gasliquid'}; equ.name = {'default','Solid domain'}; equ.usedElement = {'Lag1','Lag2'}; equ.convOn = {1,0}; equ.k = {200,400}; equ.v = {'v',0}; equ.u = {'u',0}; equ.ind = [1,1,1,1,1,1,1]; appl.equ = equ; fem.appl{2} = appl; fem.sdim = {'r','z'}; fem.frame = {'ref'}; fem.border = 1; clear units; units.basesystem = 'SI'; fem.units = units; % Scalar expressions fem.expr = {'D','exp(-(T-T_m)^2/(dT^2))/sqrt(pi*dT^2)', ... 'B','(T-T_m+dT)/(2*dT)*((T<=(T_m+dT))*(T>=(T_m-dT)))+(T>(T_m+dT))', ... 'Sr','(1-B)^2/(B^3+1e-3)*1e5[kg/(m^3*s)]*u', ... 'Sz','(1-B)^2/(B^3+1e-3)*1e5[kg/(m^3*s)]*(v-v_cast)', ... 'H','flc2hs((T-T_m)[1/K],dT[1/K])', ... 'Cp1','380[J/(kg*K)]+dH/T_m*H'}; % Descriptions clear descr descr.const= {'v_cast','Casting rate','T_in','Melt inlet temperature','dH','Latent heat','T_m','Melting temperature','T0','Ambient temperature','dT','Temperature transition zone half width'}; fem.descr = descr; % ODE Settings clear ode clear units; units.basesystem = 'SI'; ode.units = units; fem.ode=ode; % Multiphysics fem=multiphysics(fem); % Extend mesh fem.xmesh=meshextend(fem); % Solve problem fem.sol=femstatic(fem, ... 'symmetric','off', ... 'solcomp',{'vl_chns','vb_chns','T','ul_chns','p','ub_chns'}, ... 'outcomp',{'vl_chns','T','vb_chns','ul_chns','p','ub_chns'}, ... 'blocksize','auto', ... 'pname','dT', ... 'plist',[300 100 50 20], ... 'oldcomp',{}, ... 'maxiter',50, ... 'hnlin','on', ... 'linsolver','pardiso', ... 'uscale','none'); % Save current fem structure for restart purposes fem0=fem; % Plot solution postplot(fem, ... 'tridata',{'U_chns','cont','internal','unit','m/s'}, ... 'trimap','jet(1024)', ... 'solnum','end', ... 'title','dT(4)=20 Surface: Velocity field [m/s]', ... 'axis',[-0.8879239709318034,1.0036739708125941,-0.07862499952316285,1.6511249899864198]); % (Default values are not included) % Application mode 1 clear appl appl.mode.class = 'HtNavierStokes'; appl.mode.type = 'axi'; appl.module = 'HT'; appl.shape = {'shlag(1,''ul_chns'')','shlag(1,''vl_chns'')','shlag(1,''p'')','shbub(2,''ub_chns'')','shbub(2,''vb_chns'')','shlag(2,''u'')','shlag(2,''v'')'}; appl.gporder = {4,2}; appl.cporder = {1,2}; appl.assignsuffix = '_chns'; clear prop prop.analysis='static'; prop.weakcompflow='On'; appl.prop = prop; clear bnd bnd.symtype = {'sym','sym','ax','sym','sym'}; bnd.v0 = {0,0,0,0,'v_cast'}; bnd.type = {'walltype','int','sym','inlet','outlet'}; bnd.outtype = {'p','p','p','p','uv'}; bnd.velType = {'U0in','U0in','U0in','U0in','u0'}; bnd.intype = {'uv','uv','uv','p','uv'}; bnd.ind = [3,4,3,2,3,2,3,2,3,2,3,2,3,2,5,1,1,1,1,1,5,5,5]; appl.bnd = bnd; clear equ equ.eta = {0.0434,1}; equ.rhofcnT = 1; equ.gporder = {{1;1;2}}; equ.name = {'default','Solid domain'}; equ.F_z = {'-Sz',0}; equ.Tflowtype = 'nonisoT'; equ.rho = 'rho_htgh'; equ.cporder = {1,{2;2;1}}; equ.F_r = {'-Sr',0}; equ.init = {{0;'v_cast';0;0;0;0;0;0;0},0}; equ.rhofcnp = 1; equ.shape = {[1;2;3;4;5],[6;7;3]}; equ.cdon = {0,1}; equ.sdon = {0,1}; equ.rhofcnTname = 'T'; equ.usage = {1,0}; equ.ind = [1,1,1,1,1,1,1]; appl.equ = equ; fem.appl{1} = appl; % Application mode 2 clear appl appl.mode.class = 'GeneralHeat'; appl.mode.type = 'axi'; appl.module = 'HT'; appl.shape = {'shlag(1,''J'')','shlag(1,''T'')','shlag(2,''T'')'}; appl.gporder = {2,4}; appl.cporder = {1,2}; appl.assignsuffix = '_htgh'; clear prop prop.analysis='static'; clear weakconstr weakconstr.value = 'off'; weakconstr.dim = {'lm8'}; prop.weakconstr = weakconstr; appl.prop = prop; clear bnd bnd.type = {'ax','T','cont','qc','q0','q','q','q'}; bnd.radType = {'none','none','none','none','none','none','none','surf2amb'}; bnd.epsilon = {0,0,0,0,0,0,0,0.8}; bnd.Tamb = {0,0,0,0,0,0,0,'T0'}; bnd.shape = 1; bnd.T0 = {273.15,'T_in',273.15,273.15,273.15,273.15,273.15, ... 273.15}; bnd.Tinf = {273.15,273.15,273.15,273.15,273.15,'T0','T0','T0'}; bnd.h = {0,0,0,0,0,25,800,10}; bnd.ind = [1,2,1,3,1,3,1,3,1,3,1,3,1,3,4,5,5,5,5,6,7,8,8]; appl.bnd = bnd; clear equ equ.eta = {'eta_chns',0}; equ.cdtypesystem = 'system'; equ.gporder = {1,2}; equ.pA = {'101325[Pa]+p',0}; equ.rho = {8500,8700}; equ.cporder = {1,2}; equ.init = {{'T_in';0},{273.15;0}}; equ.shape = {2,3}; equ.cdon = 1; equ.sdon = {1,0}; equ.C = {'Cp1+D*dH',385}; equ.matterstate = {'userdefined','gasliquid'}; equ.name = {'default','Solid domain'}; equ.usedElement = {'Lag1','Lag2'}; equ.convOn = {1,0}; equ.k = {200,400}; equ.v = {'v',0}; equ.u = {'u',0}; equ.ind = [1,1,1,1,1,1,1]; appl.equ = equ; fem.appl{2} = appl; fem.sdim = {'r','z'}; fem.frame = {'ref'}; fem.border = 1; clear units; units.basesystem = 'SI'; fem.units = units; % Scalar expressions fem.expr = {'D','exp(-(T-T_m)^2/(dT^2))/sqrt(pi*dT^2)', ... 'B','(T-T_m+dT)/(2*dT)*((T<=(T_m+dT))*(T>=(T_m-dT)))+(T>(T_m+dT))', ... 'Sr','(1-B)^2/(B^3+1e-3)*1e5[kg/(m^3*s)]*u', ... 'Sz','(1-B)^2/(B^3+1e-3)*1e5[kg/(m^3*s)]*(v-v_cast)', ... 'H','flc2hs((T-T_m)[1/K],dT[1/K])', ... 'Cp1','380[J/(kg*K)]+dH/T_m*H'}; % Descriptions clear descr descr.const= {'v_cast','Casting rate','T_in','Melt inlet temperature','dH','Latent heat','T_m','Melting temperature','T0','Ambient temperature','dT','Temperature transition zone half width'}; fem.descr = descr; % ODE Settings clear ode clear units; units.basesystem = 'SI'; ode.units = units; fem.ode=ode; % Multiphysics fem=multiphysics(fem); % Extend mesh fem.xmesh=meshextend(fem); % Solve problem fem=adaption(fem, ... 'init',fem0.sol, ... 'symmetric','off', ... 'solcomp',{'vl_chns','T','vb_chns','ul_chns','p','ub_chns'}, ... 'outcomp',{'vl_chns','vb_chns','T','ul_chns','p','ub_chns'}, ... 'blocksize','auto', ... 'maxiter',50, ... 'hnlin','on', ... 'solver','stationary', ... 'eigselect',[1], ... 'maxt',10000000, ... 'ngen',2, ... 'rmethod','longest', ... 'resorder','auto', ... 'eefun','l2', ... 'l2scale',[1], ... 'l2staborder',[2], ... 'tppar',1.7, ... 'linsolver','pardiso', ... 'uscale','none', ... 'geomnum',1); % Save current fem structure for restart purposes fem0=fem; % Plot solution postplot(fem, ... 'tridata',{'U_chns','cont','internal','unit','m/s'}, ... 'trimap','jet(1024)', ... 'title','Surface: Velocity field [m/s]', ... 'axis',[-0.7953995850486834,0.9111495849294741,-0.07862499952316285,1.6511249899864198]); % (Default values are not included) % Application mode 1 clear appl appl.mode.class = 'HtNavierStokes'; appl.mode.type = 'axi'; appl.module = 'HT'; appl.shape = {'shlag(1,''ul_chns'')','shlag(1,''vl_chns'')','shlag(1,''p'')','shbub(2,''ub_chns'')','shbub(2,''vb_chns'')','shlag(2,''u'')','shlag(2,''v'')'}; appl.gporder = {4,2}; appl.cporder = {1,2}; appl.assignsuffix = '_chns'; clear prop prop.analysis='static'; prop.weakcompflow='On'; appl.prop = prop; clear bnd bnd.symtype = {'sym','sym','ax','sym','sym'}; bnd.v0 = {0,0,0,0,'v_cast'}; bnd.type = {'walltype','int','sym','inlet','outlet'}; bnd.outtype = {'p','p','p','p','uv'}; bnd.velType = {'U0in','U0in','U0in','U0in','u0'}; bnd.intype = {'uv','uv','uv','p','uv'}; bnd.ind = [3,4,3,2,3,2,3,2,3,2,3,2,3,2,5,1,1,1,1,1,5,5,5]; appl.bnd = bnd; clear equ equ.eta = {0.0434,1}; equ.rhofcnT = 1; equ.gporder = {{1;1;2}}; equ.name = {'default','Solid domain'}; equ.F_z = {'-Sz',0}; equ.Tflowtype = 'nonisoT'; equ.rho = 'rho_htgh'; equ.cporder = {1,{2;2;1}}; equ.F_r = {'-Sr',0}; equ.init = {{0;'v_cast';0;0;0;0;0;0;0},0}; equ.rhofcnp = 1; equ.shape = {[1;2;3;4;5],[6;7;3]}; equ.cdon = {0,1}; equ.sdon = {0,1}; equ.rhofcnTname = 'T'; equ.usage = {1,0}; equ.ind = [1,1,1,1,1,1,1]; appl.equ = equ; fem.appl{1} = appl; % Application mode 2 clear appl appl.mode.class = 'GeneralHeat'; appl.mode.type = 'axi'; appl.module = 'HT'; appl.shape = {'shlag(1,''J'')','shlag(1,''T'')','shlag(2,''T'')'}; appl.gporder = {2,4}; appl.cporder = {1,2}; appl.assignsuffix = '_htgh'; clear prop prop.analysis='static'; clear weakconstr weakconstr.value = 'off'; weakconstr.dim = {'lm8'}; prop.weakconstr = weakconstr; appl.prop = prop; clear bnd bnd.type = {'ax','T','cont','qc','q0','q','q','q'}; bnd.radType = {'none','none','none','none','none','none','none','surf2amb'}; bnd.epsilon = {0,0,0,0,0,0,0,0.8}; bnd.Tamb = {0,0,0,0,0,0,0,'T0'}; bnd.shape = 1; bnd.T0 = {273.15,'T_in',273.15,273.15,273.15,273.15,273.15, ... 273.15}; bnd.Tinf = {273.15,273.15,273.15,273.15,273.15,'T0','T0','T0'}; bnd.h = {0,0,0,0,0,25,800,10}; bnd.ind = [1,2,1,3,1,3,1,3,1,3,1,3,1,3,4,5,5,5,5,6,7,8,8]; appl.bnd = bnd; clear equ equ.eta = {'eta_chns',0}; equ.cdtypesystem = 'system'; equ.gporder = {1,2}; equ.pA = {'101325[Pa]+p',0}; equ.rho = {8500,8700}; equ.cporder = {1,2}; equ.init = {{'T_in';0},{273.15;0}}; equ.shape = {2,3}; equ.cdon = 1; equ.sdon = {1,0}; equ.C = {'Cp1+D*dH',385}; equ.matterstate = {'userdefined','gasliquid'}; equ.name = {'default','Solid domain'}; equ.usedElement = {'Lag1','Lag2'}; equ.convOn = {1,0}; equ.k = {200,400}; equ.v = {'v',0}; equ.u = {'u',0}; equ.ind = [1,1,1,1,1,1,1]; appl.equ = equ; fem.appl{2} = appl; fem.sdim = {'r','z'}; fem.frame = {'ref'}; fem.border = 1; clear units; units.basesystem = 'SI'; fem.units = units; % Scalar expressions fem.expr = {'D','exp(-(T-T_m)^2/(dT^2))/sqrt(pi*dT^2)', ... 'B','(T-T_m+dT)/(2*dT)*((T<=(T_m+dT))*(T>=(T_m-dT)))+(T>(T_m+dT))', ... 'Sr','(1-B)^2/(B^3+1e-3)*1e5[kg/(m^3*s)]*u', ... 'Sz','(1-B)^2/(B^3+1e-3)*1e5[kg/(m^3*s)]*(v-v_cast)', ... 'H','flc2hs((T-T_m)[1/K],dT[1/K])', ... 'Cp1','380[J/(kg*K)]+dH/T_m*H'}; % Descriptions clear descr descr.const= {'v_cast','Casting rate','T_in','Melt inlet temperature','dH','Latent heat','T_m','Melting temperature','T0','Ambient temperature','dT','Temperature transition zone half width'}; fem.descr = descr; % ODE Settings clear ode clear units; units.basesystem = 'SI'; ode.units = units; fem.ode=ode; % Multiphysics fem=multiphysics(fem); % Extend mesh fem.xmesh=meshextend(fem); % Solve problem fem.sol=femstatic(fem, ... 'init',fem0.sol, ... 'symmetric','off', ... 'solcomp',{'vl_chns','vb_chns','T','ul_chns','p','ub_chns'}, ... 'outcomp',{'vl_chns','T','vb_chns','ul_chns','p','ub_chns'}, ... 'blocksize','auto', ... 'pname','dT', ... 'plist',[20 10 5], ... 'oldcomp',{}, ... 'pinitstep',5, ... 'pminstep',2.5, ... 'pmaxstep',5, ... 'maxiter',50, ... 'linsolver','pardiso', ... 'uscale','none'); % Save current fem structure for restart purposes fem0=fem; % Plot solution postplot(fem, ... 'tridata',{'U_chns','cont','internal','unit','m/s'}, ... 'trimap','jet(1024)', ... 'solnum','end', ... 'title','dT(3)=5 Surface: Velocity field [m/s]', ... 'axis',[-0.8879239709318034,1.0036739708125941,-0.07862499952316285,1.6511249899864198]); % Plot solution postplot(fem, ... 'tridata',{'T','cont','internal','unit','K'}, ... 'trimap','jet(1024)', ... 'arrowdata',{'u','v'}, ... 'arrowxspacing',15, ... 'arrowyspacing',15, ... 'arrowtype','arrow', ... 'arrowstyle','proportional', ... 'arrowcolor',[1.0,0.0,0.0], ... 'solnum','end', ... 'title','dT(3)=5 Surface: Temperature [K] Arrow: Velocity field', ... 'axis',[-0.7953995850486834,0.9111495849294741,-0.07862499952316285,1.6511249899864198]); % Plot solution postplot(fem, ... 'tridata',{'B','cont','internal'}, ... 'trimap','jet(1024)', ... 'solnum','end', ... 'title','dT(3)=5 Surface: B', ... 'axis',[-0.2677221280109281,0.4066729474231441,-0.05026844231808769,0.6251331929300652]); % Plot solution postplot(fem, ... 'tridata',{'U_chns','cont','internal','unit','m/s'}, ... 'trimap','jet(1024)', ... 'flowdata',{'u','v'}, ... 'flowcolor',[1.0,0.0,0.0], ... 'flowdens','velocity', ... 'flowlines',33, ... 'solnum','end', ... 'title','dT(3)=5 Surface: Velocity field [m/s] Streamline: Velocity field', ... 'axis',[-0.274768046709493,0.413718866121709,-0.05026844231808769,0.6251331929300652]); % Plot solution postplot(fem, ... 'tridata',{'dflux_htgh','cont','internal','unit','W/m^2'}, ... 'trimap','jet(1024)', ... 'arrowdata',{'tflux_r_htgh','tflux_z_htgh'}, ... 'arrowxspacing',15, ... 'arrowyspacing',15, ... 'arrowtype','arrow', ... 'arrowstyle','proportional', ... 'arrowcolor',[1.0,0.0,0.0], ... 'solnum','end', ... 'title','dT(3)=5 Surface: Conductive heat flux [W/m^2] Arrow: Total heat flux', ... 'axis',[-0.26369588875460526,0.40264670816682124,-0.05026844231808769,0.6251331929300652]); % Plot in cross-section or along domain postcrossplot(fem,1,[20,21,22,23], ... 'lindata','ndflux_htgh', ... 'cont','internal', ... 'linxdata',{'z','unit','m'}, ... 'title','Normal conductive heat flux [W/m^2]', ... 'axislabel',{'z-coordinate [m]','Normal conductive heat flux [W/m^2]'}, ... 'refine','auto');