-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathrotinas_ept.F
312 lines (256 loc) · 11.3 KB
/
rotinas_ept.F
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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
c*************************************************************************
c
c Subrotinas para estado plano de tensões, nomeadas EPT_xxx.
c
c
c
c Eduardo Pagnussat Titello, 2019
c
c*************************************************************************
c-------------------------------------------------------------------------
c
c Tensões principais
c
SUBROUTINE EPT_TensPrinc(stress, TensPrinc, theta)
c Cálcula as tensões principais e o angulo formado por estas em EPT
c
c Parametros:
c stress (dp, 3x1, in) - Estado de tensões
c TensPrinc (dp, 2x1, out) - Tensões principais (1,2)
c theta (dp, sc, out) - Angulo de rotação (rad)
DOUBLE PRECISION stress(3), TensPrinc(2), theta
c Variáveis internas:
c P1 e P2 (dp, sc) - Partes 1 e 2 da expressão para tensões principais
c Pi (dp, sc) - Precisa mesmo?
c cost, sent (dp, sc) - Cosseno e Seno de Theta
c cost2, sent2 (dp, sc) - Cosseno2 e Seno2 de Theta
DOUBLE PRECISION P1, P2, Pi, cost, sent, cost2, sent2
Pi = 4.d0*DATAN(1.d0)
c Partes que compõem as tensões principais
P1 = (stress(1) + stress(2))/2.d0
P2 = (stress(1) - stress(2))/2.d0
c Angulo theta:
c Condições de contorno para arctan
IF(DABS(stress(3)) .LE. 10.d0**-4) THEN
theta = 0.d0
IF(P2 .LT. 0.d0) theta = Pi/2.d0
ELSEIF(DABS(P2) .LE. 10.d0**-4) THEN
theta = Pi/4.d0
IF(stress(3) .LT. 0.d0) theta = -theta
ELSE
theta = stress(3)/P2
theta = 0.5d0*DATAN(theta)
END IF
c Tensões máxima e mínima
P2 = P2**2.d0 + stress(3)**2.d0
P2 = DSQRT(P2)
TensPrinc(1) = P1 + P2
TensPrinc(2) = P1 - P2
c Nova abordagem: usando angulo para determinar tensões principais
cost = DCOS(theta)
sent = DSIN(theta)
cost2 = cost**2.d0
sent2 = sent**2.d0
c Tensões principais iniciais
TensPrinc(1) = stress(1)*cost2 + stress(2)*sent2 +
& 2.d0*stress(3)*sent*cost
TensPrinc(2) = stress(1)*sent2 + stress(2)*cost2 -
& 2.d0*stress(3)*sent*cost
c Teste se angulo não está errado, vai girando até acertar
DO WHILE(TensPrinc(1) .LT. TensPrinc(2))
theta = theta + Pi/2.d0
cost = DCOS(theta)
sent = DSIN(theta)
cost2 = cost**2.d0
sent2 = sent**2.d0
TensPrinc(1) = stress(1)*cost2 + stress(2)*sent2 +
& 2.d0*stress(3)*sent*cost
TensPrinc(2) = stress(1)*sent2 + stress(2)*cost2 -
& 2.d0*stress(3)*sent*cost
END DO
END SUBROUTINE
c-------------------------------------------------------------------------
c-------------------------------------------------------------------------
c
c "Invariantes de tensão" -> I1, J2, J3, theta
c
SUBROUTINE EPT_Invars(stress, Invars)
c Alguns invariantes e semelhantes úteis
c
c Parametros:
c stress (dp, 3x1, in) - Estado de tensões
c Invars (dp, 4x1, out) - Alguns invariantes (uteis) + ang.theta
c Sendo #1=I1, #2=J2, #3=J3, #4=theta
DOUBLE PRECISION stress(3), Invars(4), p
c Primeiro invariante de tensões principais I1
Invars(1) = stress(1) + stress(2)
c Segundo invariante de tensoes desviadoras J2
Invars(2) = (stress(1) - stress(2))**2.d0
Invars(2) = Invars(2) + ( stress(2))**2.d0
Invars(2) = Invars(2) + (-stress(1))**2.d0
Invars(2) = 1.d0/6.d0*Invars(2)
Invars(2) = Invars(2) + (stress(3))**2.d0
c Terceiro Invariante de Tensões desviadoras J3
p = Invars(1)/3.d0
Invars(3) = (stress(1)-p)*(stress(2)-p)*(-p)
Invars(3) = Invars(3) + p*(stress(3))**2.d0
c Não é invariante, mas é útil
c cos(3theta)
IF(Invars(2) .NE. 0.d0) THEN
Invars(4) = 1.5d0*DSQRT(3.d0)*Invars(3)/(Invars(2)**1.5d0)
IF(Invars(4) .LT. -1.d0) Invars(4) = -1.d0
IF(Invars(4) .GT. 1.d0) Invars(4) = 1.d0
c theta = 1/3*acos(cos(3theta))
Invars(4) = 1.d0/3.d0*DACOS(Invars(4))
ELSE
Invars(4) = 0.d0
END IF
END SUBROUTINE
c-------------------------------------------------------------------------
c-------------------------------------------------------------------------
c
c Vetor fluxo plástico de Von-Mises
c
SUBROUTINE EPT_Flux_VM(stress, flux)
c Parametros:
c stress (dp, 3x1, in) - Vetor com Tensões
c flux (dp, 3x1, out) - Vetor de fluxo
DOUBLE PRECISION stress(3), flux(3)
c Variáveis internas
c Sxx, Syy, Sxy (dp, sc) - Tensões em X, Y e XY
c Const (dp, sc) - Constante que pré multiplica a matriz conforme Schmitz
DOUBLE PRECISION Sxx, Syy, Sxy, Const
Sxx = stress(1)
Syy = stress(2)
Sxy = stress(3)
c Cálculo da constante
Const = (Sxx**2.d0 + Syy**2.d0) - (Sxx*Syy)
Const = Const + 3.d0*(Sxy**2.d0)
flux = 0.d0
IF(Const .GT. 0.d0) THEN
Const = 1.d0/(2.d0*DSQRT(Const))
flux(1) = Const*(2.d0*Sxx - Syy)
flux(2) = Const*(2.d0*Syy - Sxx)
flux(3) = Const*(6.d0*Sxy)
END IF
END SUBROUTINE
c-------------------------------------------------------------------------
c-------------------------------------------------------------------------
c
c Matriz constitutiva elástica p/ EPT
c
SUBROUTINE EPT_MatrizD(E, poisson, D)
c Parametros:
c E (dp, sc, in) - Modulo de elasticidade
c poisson (dp, sc, in) - Coeficiente de Poisson
c D (dp, 3x3, out) - Matriz constitutiva D
DOUBLE PRECISION E, poisson, D(3,3)
c
c Variáveis internas
c C1, C2, C3 (dp, sc) - Coeficientes de construção de D
DOUBLE PRECISION C1, C2, C3
c Cálcula coeficientes
C1 = E/(1.d0 - poisson**2.d0)
C2 = C1*poisson
C3 = C1*(1.d0 - poisson)/2.d0
c Monta matriz
D = 0.d0
D(1,1) = C1
D(2,2) = C1
D(1,2) = C2
D(2,1) = C2
D(3,3) = C3
END SUBROUTINE
c-------------------------------------------------------------------------
c-------------------------------------------------------------------------
c
c Rotaciona tensões em EPT
c
SUBROUTINE EPT_Rot_Tens(Stress, theta)
c Parametros:
c Stress (dp, 3x1, io) - Vetor de tensões a ser rotacionado e retorno
c theta (dp, sc, in) - Angulo de rotação
DOUBLE PRECISION Stress(3), theta
c
c Variáveis internas
c
c cost, sent (dp, sc) - Seno e Cosseno de theta
c cos2t, sen2t(dp, sc) - Seno**2 e Cosseno**2 de theta
c StressRot (dp 3x1) - Vetor de tensões rotacionadas
DOUBLE PRECISION cost, sent, StressRot(3), cos2t, sen2t
c Cálcula coeficientes
cost = DCOS(theta)
sent = DSIN(theta)
cos2t = cost**2.d0
sen2t = sent**2.d0
c Monta vetor rotacionado
StressRot(1) = Stress(1)*cos2t + Stress(2)*sen2t +
& 2.d0*Stress(3)*cost*sent
StressRot(2) = Stress(1)*sen2t + Stress(2)*cos2t -
& 2.d0*Stress(3)*cost*sent
StressRot(3) = (Stress(2) - Stress(1))*cost*sent +
& Stress(3)*(cos2t - sen2t)
c Passa vetor rotacionado para entrada/saida
Stress(1) = StressRot(1)
Stress(2) = StressRot(2)
Stress(3) = StressRot(3)
END SUBROUTINE
c-------------------------------------------------------------------------
c-------------------------------------------------------------------------
c
c Rotaciona deformações em EPT
c
SUBROUTINE EPT_Rot_Def(Strain, theta)
c Parametros:
c Strain (dp, 3x1, io) - Vetor de deformações a ser rotacionado e retorno
c theta (dp, sc, in) - Angulo de rotação
DOUBLE PRECISION Strain(3), theta
c Altera Strain(3) e usa rotina de tensões ;D
Strain(3) = Strain(3)/2.d0
CALL EPT_Rot_Tens(Strain, theta)
Strain(3) = Strain(3)*2.d0
END SUBROUTINE
c-------------------------------------------------------------------------
c-------------------------------------------------------------------------
c
c Rotaciona matriz constitutiva EPT
c
SUBROUTINE EPT_Rot_Const(dsdePl, theta)
c Parametros:
c dsdePl (dp, 3x3, io) - Matriz constitutiva EPT
c theta (dp, sc, in) - Angulo de rotação
DOUBLE PRECISION dsdePl(3,3), theta
c
c Variáveis internas
c cost, sent (dp, sc) - Seno e Cosseno de theta
c mrot (dp 3x3) - Matriz de rotação
DOUBLE PRECISION cost, sent, mrot(3,3), dsdePlRot(3,3)
c Cálcula coeficientes e monta matriz de rotação
cost = DCOS(theta)
sent = DSIN(theta)
c Monta matriz rotacionada
dsdePlRot(1,1) = sent*(dsdePl(2,2)*sent+dsdePl(1,2)*cost)
& +cost*(dsdePl(1,2)*sent+dsdePl(1,1)*cost)
dsdePlRot(1,2) = sent*(dsdePl(2,2)*cost-dsdePl(1,2)*sent)
& +cost*(dsdePl(1,2)*cost-dsdePl(1,1)*sent)
dsdePlRot(1,3) = dsdePl(2,3)*sent+dsdePl(1,3)*cost
dsdePlRot(2,1) = cost*(dsdePl(2,2)*sent+dsdePl(1,2)*cost)
& -sent*(dsdePl(1,2)*sent+dsdePl(1,1)*cost)
dsdePlRot(2,2) = cost*(dsdePl(2,2)*cost-dsdePl(1,2)*sent)
& -sent*(dsdePl(1,2)*cost-dsdePl(1,1)*sent)
dsdePlRot(2,3) = dsdePl(2,3)*cost-dsdePl(1,3)*sent
dsdePlRot(3,1) = dsdePl(2,3)*sent+dsdePl(1,3)*cost
dsdePlRot(3,2) = dsdePl(2,3)*cost-dsdePl(1,3)*sent
dsdePlRot(3,3) = dsdePl(3,3)
c Passa para saida
dsdePl(1,1) = dsdePlRot(1,1)
dsdePl(1,2) = dsdePlRot(1,2)
dsdePl(1,3) = dsdePlRot(1,3)
dsdePl(2,1) = dsdePlRot(2,1)
dsdePl(2,2) = dsdePlRot(2,2)
dsdePl(2,3) = dsdePlRot(2,3)
dsdePl(3,1) = dsdePlRot(3,1)
dsdePl(3,2) = dsdePlRot(3,2)
dsdePl(3,3) = dsdePlRot(3,3)
END SUBROUTINE
c-------------------------------------------------------------------------