-
Notifications
You must be signed in to change notification settings - Fork 0
/
onsasExample_dynamicVonMises.m
124 lines (95 loc) · 3.8 KB
/
onsasExample_dynamicVonMises.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
120
121
122
123
124
%md
%mdBefore defining the structs, the workspace is cleaned, the ONSAS directory is added to the path and scalar auxiliar parameters are defined.
close all, clear all ; addpath( genpath( [ pwd '/../../src'] ) );
% scalar parameters
%%% Parametros Estructura
rho = 7850; % kg/m3 (acero)
Lx = .374/2;
L0 = .205 ;
Lc = .240; %m
Ic = .0254*.0032^3/12; %m4
Ac = .0254*.0032; %m2
E = 200000e6 ; %Pa (acero)
kc = 3*E*Ic/Lc^3; %N/m
mConc = 1.4; %kg Pandeo incipiente en 1.4
c = 0; %kg/s (amortiguamiento por friccion juntas y arrastre pesa)
g = 9.81; %m/s2
tf = 1.0;
dt = .000025; % sec
[u, normalForce, times ] = centralDiffDynVonMises(rho, Lx, L0, Lc, Ic, Ac, E, kc, mConc, c, g, tf, dt);
mb = L0*Ac*rho; %kg
Lz = sqrt( L0^2 - Lx^2 ); %m
rhoBarraMasa = mConc*.5 / (Lz*Ac);
M = [mb 0 ; 0 (mb+mConc)/2] ;
nu = .3 ;
materials(1).modelName = 'elastic-rotEngStr' ;
materials(1).modelParams = [ E nu ] ;
materials(1).density = rho ;
materials(2).modelName = 'elastic-rotEngStr' ;
materials(2).modelParams = [ 1e10*E nu ] ;
materials(2).density = rhoBarraMasa ;
materials(3).modelName = 'elastic-rotEngStr' ;
materials(3).modelParams = [ kc*L0/Ac nu ] ;
materials(3).density = rho ;
elements(1).elemType = 'node' ;
elements(2).elemType = 'truss';
elements(2).elemCrossSecParams{1,1} = 'rectangle' ;
elements(2).elemCrossSecParams{2,1} = [sqrt(Ac) sqrt(Ac) ] ;
elements(2).massMatType = 'lumped' ;
boundaryConds(1).imposDispDofs = [ 3 5 ] ;
boundaryConds(1).imposDispVals = [ 0 0 ] ;
boundaryConds(2).imposDispDofs = [ 1 3 ] ;
boundaryConds(2).imposDispVals = [ 0 0 ] ;
boundaryConds(2).loadsCoordSys = 'global' ;
boundaryConds(2).loadsTimeFact = @(t) 1 ;
boundaryConds(2).loadsBaseVals = [ 0 0 0 0 -(mConc+mb)/2*g 0 ] ;
boundaryConds(3).imposDispDofs = [ 1 3 ] ;
boundaryConds(3).imposDispVals = [ 0 0 ] ;
boundaryConds(4).imposDispDofs = [ 1 3 5 ] ;
boundaryConds(4).imposDispVals = [ 0 0 0 ] ;
mesh.nodesCoords = [ 0 0 0 ; ...
Lx 0 Lz ; ...
Lx 0 Lz-Lz; ...
-L0 0 0 ] ;
mesh.conecCell = { } ;
mesh.conecCell{ 1, 1 } = [ 0 1 1 1 ] ;
mesh.conecCell{ 2, 1 } = [ 0 1 2 2 ] ;
mesh.conecCell{ 3, 1 } = [ 0 1 3 3 ] ;
mesh.conecCell{ 4, 1 } = [ 0 1 4 4 ] ;
mesh.conecCell{ 4+1, 1 } = [ 1 2 0 1 2 ] ;
mesh.conecCell{ 4+2, 1 } = [ 2 2 0 2 3 ] ;
mesh.conecCell{ 4+3, 1 } = [ 3 2 0 1 4 ] ;
initialConds = struct() ;
analysisSettings.methodName = 'newmark' ;
%md and the following parameters correspond to the iterative numerical analysis settings
analysisSettings.deltaT = 0.0005 ;
analysisSettings.finalTime = 1 ;
analysisSettings.stopTolDeltau = 1e-8 ;
analysisSettings.stopTolForces = 1e-8 ;
analysisSettings.stopTolIts = 10 ;
%md
%md### otherParams
otherParams.problemName = 'dynamicVonMisesTruss';
otherParams.plotsFormat = 'vtk' ;
%md
%md### Analysis case 1: NR with Rotated Eng Strain
%md In the first case ONSAS is run and the solution at the dof of interest is stored.
[matUs, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
uONSAS = [ matUs(1,:) ; matUs(6+5,:) ] ;
timesONSAS = 0:analysisSettings.deltaT:analysisSettings.finalTime ;
subplot(3,1,1)
plot(times(1:10:end),1000*u(1,1:10:end))
hold on
plot(timesONSAS,1000*uONSAS(1,:),'r' )
xlabel('t [s]'); ylabel('u_1 [mm]');
axis( [0 2 1e3*min(u(1,:))*1.1 1e3*max(u(1,:))*1.1] );
subplot(3,1,2)
hold on
plot(times(1:10:end),1000*u(2,1:10:end))
plot(timesONSAS,1000*uONSAS(2,:),'r' )
xlabel('t [s]'); ylabel('u_2 [mm]');
axis([0 2 1e3*min(u(2,:))*1.1 1e3*max(u(2,:))*1.1]);
subplot(3,1,3)
plot(times(1:10:end),normalForce(1:10:end))
xlabel('t [s]'); ylabel('Directa [N]')
axis([0 2 min(normalForce)*1.1 max(normalForce)*1.1]);