% COMSOL Multiphysics Model M-file % Generated by COMSOL 3.5a (COMSOL 3.5.0.606, $Date: 2009/04/29 09:11:29 $) % Some geometry objects are stored in a separate file. % The name of this file is given by the variable 'flbinaryfile'. flclear fem % COMSOL version clear vrsn vrsn.name = 'COMSOL 3.5'; vrsn.ext = 'a'; vrsn.major = 0; vrsn.build = 606; vrsn.rcs = '$Name: v35ap $'; vrsn.date = '$Date: 2009/04/29 09:11:29 $'; fem.version = vrsn; % Geometry g1=sphere3('1e-4','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0'); g2=sphere3('2e-4','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0'); g3=geomcomp({g1,g2},'ns',{'g1','g2'},'sf','g1+g2','face','none','edge','all'); g4=block3('1e-3','1e-3','1e-3','base','center','pos',{'0','0','-.5e-4'},'axis',{'0','0','1'},'rot','0'); g4=move(g4,[0,0,-.5e-4]); g5=block3('1e-3','1e-3','1e-3','base','center','pos',{'0','0','-.5e-3'},'axis',{'0','0','1'},'rot','0'); g6=geomcomp({g3,g5},'ns',{'g3','g5'},'sf','g3*g5','face','none','edge','all'); parr={point3(0,0,0)}; g7=geomcoerce('point',parr); % Analyzed geometry clear p s p.objs={g7}; p.name={'PT1'}; p.tags={'g7'}; s.objs={g6}; s.name={'CO2'}; s.tags={'g6'}; fem.draw=struct('p',p,'s',s); fem.geom=geomcsg(fem); % flbinaryfile='nickel_model.mphm'; % Constants fem.const = {'r0','4e-7', ... 't0','1e-9', ... 'r1','1e-4', ... 'r2','2e-4', ... 'I0','1e12'}; % Geometry % fem.geom = flbinary('g1','geom',flbinaryfile); % Application mode 1 clear appl appl.mode.class = 'SmeSolid3'; appl.module = 'SME'; appl.gporder = 4; appl.cporder = 2; appl.assignsuffix = '_smsld'; clear prop prop.analysis='time'; appl.prop = prop; clear equ equ.nu = 0.35; equ.Tflag = 1; equ.rho = 8912; equ.betadK = {'((r-r1)/(r2-r1))^3*1e-6',0}; equ.dampingtype = {'Rayleigh','nodamping'}; equ.D = {{1.6360e11,0.9230e11,0.6792e11,0,0,0;0.9230e11, ... 1.6360e11,0.6792e11,0,0,0;0.6792e11,0.6792e11,1.8520e11, ... 0,0,0;0,0,0,0.3565e11,0,0;0,0,0,0,0.4705e11, ... 0;0,0,0,0,0,0.4705e11}}; equ.D = C_cell; equ.E = 110e9; equ.Tempref = 273.15; equ.materialmodel = 'aniso'; equ.Temp = 'T'; equ.alpha = 17e-6; equ.name = {'','default'}; equ.alphavector = {{1.34e-5;1.34e-5;1.34e-5}}; equ.ind = [1,2]; appl.equ = equ; fem.appl{1} = appl; % Application mode 2 clear appl appl.mode.class = 'HeatTransfer'; appl.assignsuffix = '_ht'; clear prop clear weakconstr weakconstr.value = 'off'; weakconstr.dim = {'lm4'}; prop.weakconstr = weakconstr; appl.prop = prop; clear bnd bnd.q0 = {0,0,'q'}; bnd.type = {'q0','cont','q'}; bnd.ind = [1,1,1,2,3,2,1,2,1,2]; appl.bnd = bnd; clear equ equ.C = '445[J/(kg*K)]'; equ.k = '90.7[W/(m*K)]'; equ.name = 'default'; equ.rho = '8912[kg/m^3]'; equ.ind = [1,1]; appl.equ = equ; fem.appl{2} = appl; fem.frame = {'ref'}; fem.border = 1; fem.outform = 'general'; clear units; units.basesystem = 'SI'; fem.units = units; % Subdomain settings clear equ equ.ind = [1,2]; equ.dim = {'u','v','w','p','T'}; % Subdomain expressions equ.expr = {'r',{'sqrt(x^2+y^2+z^2)',''}}; fem.equ = equ; % Boundary settings clear bnd bnd.ind = [1,1,1,1,2,1,1,1,1,1]; % Boundary expressions bnd.expr = {'rxy',{'','sqrt(x^2+y^2)'}, ... 'q',{'','I0*t*exp(-t/t0)*exp(-rxy^2/r0^2)/t0'}}; fem.bnd = bnd; % Library materials clear lib lib.mat{1}.name='Ni'; lib.mat{1}.varname='mat1'; lib.mat{1}.variables.nu='0.31'; lib.mat{1}.variables.E='219e9[Pa]'; lib.mat{1}.variables.sigma='13.8e6[S/m]'; lib.mat{1}.variables.alpha='13.4e-6[1/K]'; lib.mat{1}.variables.C='445[J/(kg*K)]'; lib.mat{1}.variables.rho='8900[kg/m^3]'; lib.mat{1}.variables.k='90.7[W/(m*K)]'; fem.lib = lib; % ODE Settings clear ode clear units; units.basesystem = 'SI'; ode.units = units; fem.ode=ode; % Multiphysics fem=multiphysics(fem); % Initialise the mesh fem.mesh=meshinit(fem, ... 'hauto',mesh_level, ... 'hmaxvtx',[7,1e-7], ... 'hgradvtx',[7,10], ... 'hmaxfac',[5,4.e-6], ... 'hgradfac',[5,10]); fem.xmesh = meshextend(fem); % Evaluate initial value init = asseminit(fem,'blocksize','auto'); % Update model fem.sol = init; % Save current fem structure for restart purposes fem0=fem; % Solve the model fem.sol=femtime(fem, ... 'solcomp',{'w','T','v','u'}, ... 'outcomp',{'w','v','T','u'}, ... 'blocksize','auto', ... 'odesolver','genalpha', ... 'tlist',times, ... 'tout','tlist', ... 'atol',{'1e-16'}, ... 'incrdelay','off', ... 'nlsolver','manual', ... 'linsolver','spooles', ... 'ntolfact',1, ... 'maxiter',4, ... 'dtech','const', ... 'damp',1.0, ... 'jtech','minimal'); % Plot solution postplot(fem, ... 'tridata',{'w','cont','internal','unit','m'}, ... 'tridlim',[-5e-12 5e-12], ... 'trimap','Rainbow', ... 'title','Time=0 Boundary: z-displacement [m]', ... 'geom','off', ... 'grid','on', ... 'campos',[-0.0021924694954226307,9.008658625091105E-4,0.001798428853438161], ... 'camtarget',[-3.2522338210578045E-5,-8.03394678989124E-6,-7.467791340396796E-5], ... 'camup',[0.6096405526049752,-0.1537799791553065,0.7776182319303567], ... 'camva',14.673518253685748);