-
-
Notifications
You must be signed in to change notification settings - Fork 3
/
Benchmark_init.m
119 lines (103 loc) · 2.6 KB
/
Benchmark_init.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.3 (COMSOL 3.3.0.405, $Date: 2006/08/31 18:03:47 $)
% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.3';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 405;
vrsn.rcs = '$Name: $';
vrsn.date = '$Date: 2006/08/31 18:03:47 $';
fem.version = vrsn;
% Geometry
%g5=rect2(0.405,0.015,'base','corner','pos',[0.0,-0.0075]);
g5=rect2(0.4,0.02,'base','corner','pos',[0.0,-0.01]);
g8=circ2(0.05,'pos',[0.0,0.0]);
% Geometry
g1=geomcomp({g5,g8},'ns',{'g5','g8'},'sf','g5+g8','edge','none');
% Geometry
g2=geomdel(g1);
% Analyzed geometry
clear s
s.objs={g2};
s.name={'CO2'};
s.tags={'g2'};
fem.draw=struct('s',s);
fem.geom=geomcsg(fem);
% Initialize mesh
fem.mesh=meshinit(fem, ...
'hauto',5);
% ======= HERE CAN WE DEFINE THE RESOLUTION OF THE MESH !!!! =========
% Refine mesh
%fem.mesh=meshrefine(fem, ...
% 'mcase',0, ...
% 'rmethod','regular');
% Refine mesh
%fem.mesh=meshrefine(fem, ...
% 'mcase',0, ...
% 'rmethod','regular');
% (Default values are not included)
% Application mode 1
clear appl
appl.mode.class = 'SmePlaneStrain';
appl.module = 'SME';
appl.gporder = 4;
appl.cporder = 2;
appl.assignsuffix = '_smpn';
clear prop
prop.analysis='time';
prop.largedef='on';
appl.prop = prop;
clear bnd
bnd.Fx = 'force_x(x,y)';
bnd.constrcond = {'free','free','fixed','fixed'};
%bnd.constrcond = {'fixed','fixed','fixed','fixed'};
bnd.Hy = {0,0,1,1};
bnd.Hx = {0,0,1,1};
bnd.loadtype = {'length','length','length','length'};
bnd.Fy = 'force_y(x,y)';
bnd.ind = [2,1,1,3,3,4,4];
appl.bnd = bnd;
%clear equ
%equ.loadtype = 'area';
%equ.constrtype = 'general';
%equ.dampingtype = 'nodamping';
%equ.ind = [1];
clear equ
%equ.rho = {'1000[kg/m^3]'}; !!!!
equ.rho = {'10000[kg/m^3]'};
equ.loadtype = {'area'};
equ.constrtype = {'general'};
equ.dampingtype = 'nodamping'; %if you want to have damping turn on this
equ.constrcond = {'free'};
%equ.batadK={0.05};
%equ.alphadM={20};
equ.nu = {0.4};
equ.Fy = {0,0};
%equ.E = {'1.4E6[Pa]'};
equ.E = {'1.4E6[Pa]'};
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;
% Multiphysics
fem=multiphysics(fem);
% Extend mesh
fem.xmesh=meshextend(fem);
% store and initialize variables
init = asseminit(fem);
fem.sol = init;
fem0 = fem;
% ------ Solve problem -----------
fem.sol=femtime(fem, ...
'solcomp',{'u','v'}, ...
'outcomp',{'u','v'}, ...
'tlist',[0:0:0], ...
'atol',{'0.000010'}, ...
'rtol',0.0001, ...
'tout','tlist', ...
'tsteps','strict');