Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Matous_Const_Model.f90
Go to the documentation of this file.
1 !*********************************************************************
2 !* Illinois Open Source License *
3 !* *
4 !* University of Illinois/NCSA *
5 !* Open Source License *
6 !* *
7 !* Copyright@2008, University of Illinois. All rights reserved. *
8 !* *
9 !* Developed by: *
10 !* *
11 !* Center for Simulation of Advanced Rockets *
12 !* *
13 !* University of Illinois *
14 !* *
15 !* www.csar.uiuc.edu *
16 !* *
17 !* Permission is hereby granted, free of charge, to any person *
18 !* obtaining a copy of this software and associated documentation *
19 !* files (the "Software"), to deal with the Software without *
20 !* restriction, including without limitation the rights to use, *
21 !* copy, modify, merge, publish, distribute, sublicense, and/or *
22 !* sell copies of the Software, and to permit persons to whom the *
23 !* Software is furnished to do so, subject to the following *
24 !* conditions: *
25 !* *
26 !* *
27 !* @ Redistributions of source code must retain the above copyright *
28 !* notice, this list of conditions and the following disclaimers. *
29 !* *
30 !* @ Redistributions in binary form must reproduce the above *
31 !* copyright notice, this list of conditions and the following *
32 !* disclaimers in the documentation and/or other materials *
33 !* provided with the distribution. *
34 !* *
35 !* @ Neither the names of the Center for Simulation of Advanced *
36 !* Rockets, the University of Illinois, nor the names of its *
37 !* contributors may be used to endorse or promote products derived *
38 !* from this Software without specific prior written permission. *
39 !* *
40 !* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
41 !* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
42 !* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
43 !* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
44 !* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
45 !* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
46 !* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
47 !* USE OR OTHER DEALINGS WITH THE SOFTWARE. *
48 !*********************************************************************
49 !* Please acknowledge The University of Illinois Center for *
50 !* Simulation of Advanced Rockets in works and publications *
51 !* resulting from this software or its derivatives. *
52 !*********************************************************************
53 SUBROUTINE matous_const_model(cd, cb, &
54  p1, p2, yin, softparam, km, k2, gm, g2, vm, nu2, a_eta, a_zeta,&
55  cm, c2, lo, lm, l2, mm, m2, l_bar, strain, delstrain,&
56  s11, s22, s33, s12, s13, s23,ielem)
57 
58 
59  USE precision
60 
61  IMPLICIT NONE
62 
63 ! input
64 
65  INTEGER :: ielem
66  REAL(KIND=wp) :: cd, cb, softparam
67 
68  REAL(KIND=wp), DIMENSION(1:6,1:6) :: lo, l_bar
69  REAL(KIND=wp) :: km, k2, gm, g2, vm, nu2, a_eta, a_zeta, cm
70 
71 ! added
72 
73  REAL(KIND=wp) :: p1, p2, yin, c2, y
74 
75  REAL(KIND=wp) :: gy, tot
76 
77  REAL(KIND=wp), DIMENSION(1:6,1:6) :: scriptj, tmpm
78  INTEGER :: kk,ii,jj
79  REAL(KIND=wp), DIMENSION(1:6) :: strain, scriptjstrain, totv
80 
81  REAL(KIND=wp), DIMENSION(1:6,1:6) :: lm, l2, mm, m2
82 
83  REAL(KIND=wp), DIMENSION(1:6,1:6) :: dbm,dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd
84 
85  REAL(KIND=wp), DIMENSION(1:6,1:6) :: a2, am
86 
87  REAL(KIND=wp), DIMENSION(1:6,1:6) :: llbar
88 
89  REAL(KIND=wp) :: skk
90 
91  REAL(KIND=wp), DIMENSION(1:6,1:6) :: scriptam, scriptab, scriptad
92 
93  REAL(KIND=wp), DIMENSION(1:6,1:6) :: scriptk
94 
95  REAL(KIND=wp), DIMENSION(1:6,1:6) :: rinv, ninv, w, qm, q2, r, n
96 
97  REAL(KIND=wp), DIMENSION(1:6) :: pt
98 
99  REAL(KIND=wp) :: h, etke, dcd
100 
101  REAL(KIND=wp), PARAMETER :: tol = 1.e-6
102 
103  REAL(KIND=wp), DIMENSION(1:6,1:6) :: checkstress1, checkstress2, test
104 
105  REAL(KIND=wp), DIMENSION(1:6) :: mu_d, strain_m, strain_d, strain_b
106  REAL(KIND=wp), DIMENSION(1:6) :: stress_d, stress_b, stress_m
107 
108  REAL(KIND=wp), DIMENSION(1:6,1:6) :: dident = reshape( &
109  (/1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
110  0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
111  0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
112  0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, &
113  0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, &
114  0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp /),(/6,6/) )
115 
116 
117  REAL(KIND=wp) :: delscriptk
118 
119  REAL(KIND=wp), DIMENSION(1:6) :: delstrain
120 
121  REAL(KIND=wp) :: sum1, sum2, scalar1, scalar2, cd_gauss_old, nor, nor1
122 
123  integer :: icnt
124 ! local
125  REAL(KIND=wp) :: dcd_old
126 
127  REAL(KIND=wp) :: s11, s22, s33, s12, s13, s23
128 
129  REAL(KIND=wp) :: cd_n, gytmp
130 
131  INTEGER,PARAMETER :: monnd = 175
132 
133 
134  dcd_old = 0.0_wp
135 
136  cb = c2 - cd
137 
138  cd_n = cd
139 
140 ! PRINT*,'&&&', ielem,cd, Dcd, cd_n
141  CALL a_d_tensors(lo, lm, l2, mm, m2, cm, cb, cd, &
142  dbm,dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd, a2, am, l_bar)
143 ! PRINT*,'ljlk'
144 
145 ! present context
146 ! Lb = L2 ; Ad = A2
147 ! Ld = L2 ; Ab = A2
148  !Ld, Lb
149 
150 
151  CALL operatorj(lm, l2, l2, l2, l_bar, mm, m2, cm, cb, cd, &
152  dbm, dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd, &
153  am, a2, a2, a2, km, k2, gm, g2, vm, nu2, a_eta, a_zeta, &
154  scriptj, n, ninv, r, rinv, q2, qm, w)
155 
156 
157 
158 ! PRINT*,'Am,A2',Am(1,1),A2(1,1)
159 ! PRINT*,'Lbar', L_bar(1,1)
160 
161 ! ScriptK = 2*ScriptJ*Rinv*(L2*Q2-W*Lm*Qm)*TRANSPOSE(A2) * L2 * Ninv*R
162 
163  scriptk = matmul(l2,q2) - matmul(w,matmul(lm,qm))
164 
165  scriptk = matmul(rinv,scriptk)
166 
167  scriptk = 2.d0*matmul(scriptj,scriptk)
168 
169  scriptk = matmul(scriptk,transpose(a2) )
170 
171  scriptk = matmul(scriptk,l2 )
172 
173  scriptk = matmul(scriptk,ninv )
174 
175  scriptk = matmul(scriptk,r )
176 
177  scriptam = am + matmul(dmd,matmul(ninv,r))
178  scriptab = a2 + matmul(dbd,matmul(ninv,r))
179  scriptad = a2 + matmul(ddd-dident, matmul(ninv,r))
180 
181  llbar = cm*matmul(lm,scriptam) + cb*matmul(l2,scriptab) + cd*matmul(l2,scriptad)
182 
183  s11 = llbar(1,1)*delstrain(1) + llbar(1,2)*delstrain(2) + llbar(1,3)*delstrain(3)
184  s22 = llbar(2,1)*delstrain(1) + llbar(2,2)*delstrain(2) + llbar(2,3)*delstrain(3)
185  s33 = llbar(3,1)*delstrain(1) + llbar(3,2)*delstrain(2) + llbar(3,3)*delstrain(3)
186 
187  skk = ( s11 + s22 + s33 )/3.d0
188 
189 ! _ T _
190 ! Y = - 0.5 * E * J * E ... (24)
191 !
192 
193  CALL vecmatvecmul( strain, scriptj, strain, 6, y)
194 
195  CALL vecmatvecmul( strain, scriptj, delstrain, 6, scalar1)
196 
197  y = -0.5d0 * y
198 
199 ! The function G(y) characterizes the damage process given by the Weibull distribution
200  IF(y.LT.yin)THEN
201  gy = 0.d0
202  ELSE
203  gy = c2 - c2 * exp( - ( (y - yin)/(p1*yin) )**p2)
204  ENDIF
205 
206 ! IF(ielem.EQ.10) PRINT*,Y, cd, Gy, SoftParam, scalar1, Skk
207  IF((gy - softparam).GT.0.d0.AND.scalar1<0.d0.AND.skk>0.d0.AND.cd.LT.c2)THEN ! damage
208 ! IF((Gy - SoftParam).GT.0.d0)THEN ! damage
209 
210 
211 ! IF(i.EQ.357) PRINT*,'DAMAGE************',
212 
213  icnt = 0
214  converge: DO
215 
216 !!$ IF(Y.LT.Yin) THEN
217 !!$ H = 0.d0
218 !!$ ELSE
219 !!$ H = c2*EXP( - ( (Y - Yin)/(p1*Yin) )**p2)*p2*(Y-Yin)**(p2-1.d0)/(p1*Yin)**p2
220 !!$ ENDIF
221 !!$
222 !!$ CALL VecMatVecMul( strain, ScriptK, strain, 6, scalar1)
223 !!$ scalar1 = H/(1.d0 + 0.5d0*H*scalar1)
224 !!$
225 !!$ CALL VecMatVecMul( strain, ScriptJ, Delstrain, 6, scalar2)
226 !!$
227 !!$ Dcd = - scalar1 * scalar2
228 
229  dcd = gy - softparam
230 
231 ! IF(i.EQ.MonNd) PRINT*,'cd(1,i) , Dcd', cd(1,i) , Dcd,cd(1,i) + Dcd
232 
233 ! IF(i.EQ.MonNd) PRINT*,'sc1, sc2',scalar1, scalar2
234 
235  cd = cd_n + dcd
236 
237 
238  cb = c2 - cd
239 
240  CALL a_d_tensors(lo, lm, l2, mm, m2, cm, cb, cd, &
241  dbm,dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd, a2, am, l_bar)
242 
243 
244  CALL operatorj(lm, l2, l2, l2, l_bar, mm, m2, cm, cb, cd, &
245  dbm, dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd, &
246  am, a2, a2, a2, km, k2, gm, g2, vm, nu2, a_eta, a_zeta, &
247  scriptj, n, ninv, r, rinv, q2, qm, w)
248 
249 ! ScriptK = 2*ScriptJ*Rinv*(L2*Q2-W*Lm*Qm)*TRANSPOSE(A2) * L2 * Ninv*R
250 
251  scriptk = matmul(l2,q2) - matmul(w,matmul(lm,qm))
252  scriptk = matmul(rinv,scriptk)
253  scriptk = 2.d0*matmul(scriptj,scriptk)
254  scriptk = matmul(scriptk,transpose(a2) )
255  scriptk = matmul(scriptk,l2 )
256  scriptk = matmul(scriptk,ninv )
257  scriptk = matmul(scriptk,r )
258 
259  nor = sqrt( (dcd - dcd_old)**2)
260 
261  IF(icnt.EQ.0)THEN
262  nor1 = nor
263  nor = nor/nor1
264  ELSE
265  nor = nor/nor1
266  ENDIF
267 
268 ! _ T _
269 ! Y = - 0.5 * E * J * E ... (24)
270 !
271 
272  CALL vecmatvecmul( strain, scriptj, strain, 6, y)
273 
274  y = -0.5d0 * y
275 
276 
277  IF(y.LT.yin)THEN
278  gy = 0.d0
279  ELSE
280  gy = c2 - c2 * exp( - ( (y - yin)/(p1*yin) )**p2)
281  ENDIF
282 
283 !!$ IF(ielem.EQ.MonNd)THEN
284 !!$ PRINT*, 'nor,cd,gy =', nor, cd, Gy
285 !!$ PRINT*, 'Y,H,Dcd,kappaH',Y, H, Dcd,DelScriptK
286 !!$ ENDIF
287 
288  IF(y.LT.yin)THEN
289  cd = cd_n
290  EXIT
291  ENDIF
292 
293  icnt = icnt + 1
294 
295  IF(icnt.GT.100)THEN
296  print*,'Did not converge'
297  print*,'Element',ielem
298  stop
299  ENDIF
300  IF( nor .LE. tol) exit! converged
301 
302  dcd_old = dcd
303 
304  ENDDO converge
305 
306  IF(cd.GT.0.99999d0*c2) cd = c2
307 
308 ! Store cd at gauss point
309 
310  cb = c2 - cd
311 
312  IF(dcd.LT.0.d0.OR.cd.GT.c2) THEN
313  print*,'Dcd < 0 or cd > c2',ielem
314  print*,' Dcd, cd', dcd, cd
315  stop
316  ENDIF
317 
318  CALL a_d_tensors(lo, lm, l2, mm, m2, cm, cb, cd, &
319  dbm,dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd, a2, am, l_bar)
320 
321 
322  CALL operatorj(lm, l2, l2, l2, l_bar, mm, m2, cm, cb, cd, &
323  dbm, dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd, &
324  am, a2, a2, a2, km, k2, gm, g2, vm, nu2, a_eta, a_zeta, &
325  scriptj, n, ninv, r, rinv, q2, qm, w)
326 
327 
328 ! ScriptK = 2*ScriptJ*Rinv*(L2*Q2-W*Lm*Qm)*TRANSPOSE(A2) * L2 * Ninv*R
329 
330  scriptk = matmul(l2,q2) - matmul(w,matmul(lm,qm))
331  scriptk = matmul(rinv,scriptk)
332  scriptk = 2.d0*matmul(scriptj,scriptk)
333  scriptk = matmul(scriptk,transpose(a2) )
334  scriptk = matmul(scriptk,l2 )
335  scriptk = matmul(scriptk,ninv )
336  scriptk = matmul(scriptk,r )
337 
338 
339 !!$! COMPUTE:
340 !!$
341 !!$! _ T _ _ T -
342 !!$! DK = - e * J * De - .5 * e * K * e * Dcd
343 !!$!
344 !!$
345 !!$! a)
346 !!$! _ T _
347 !!$! e * J * De
348 !!$
349 !!$
350 !!$ CALL VecMatVecMul( strain, ScriptJ, Delstrain,6, scalar1)
351 !!$
352 !!$! b)
353 !!$! _ T -
354 !!$! .5 * e * K * e * Dcd
355 !!$!
356 !!$
357 !!$
358 !!$ CALL VecMatVecMul( strain, ScriptK, strain,6, scalar2)
359 !!$
360 !!$
361 !!$! DK
362 !!$
363 !!$ DelScriptK = - scalar1 - 0.5d0* scalar2 * Dcd
364 !!$
365 
366  CALL vecmatvecmul( strain, scriptj, strain, 6, y)
367 
368  y = -0.5d0 * y
369 
370  IF(y.LT.yin) THEN
371  gy = 0.d0
372  ELSE
373  gy = c2 - c2 * exp( - ( (y - yin)/(p1*yin) )**p2)
374  ENDIF
375 
376 
377  ! Xi = DK * H
378 
379  ! Update Softening Parameter
380  ! SoftParam = SoftParam + DelScriptK*H
381 
382  IF(gy.gt.softparam) softparam = gy
383 
384 ! IF(i.EQ.MonNd) PRINT*,'cd after',cd,SoftParam, Gy
385 
386  scriptam = am + matmul(dmd,matmul(ninv,r))
387  scriptab = a2 + matmul(dbd,matmul(ninv,r))
388  scriptad = a2 + matmul((ddd-dident),matmul(ninv,r))
389 
390  ENDIF ! end of damage loop
391 
392 ! Check micro and macro stresses
393 
394  tmpm = matmul(ninv,r)
395 
396  DO ii = 1, 6
397  tot = 0.d0
398  DO jj = 1, 6
399  tot = tot + tmpm(ii,jj)*strain(jj)
400  ENDDO
401  mu_d(ii) = tot
402  ENDDO
403 
404  DO ii = 1, 6
405  tot = 0.d0
406  DO jj = 1, 6
407  tot = tot + am(ii,jj)*strain(jj)
408  ENDDO
409  strain_m(ii) = tot
410  ENDDO
411 
412  DO ii = 1, 6
413  tot = 0.d0
414  DO jj = 1, 6
415  tot = tot + dmd(ii,jj)*mu_d(jj)
416  ENDDO
417  strain_m(ii) = strain_m(ii) + tot
418  ENDDO
419 
420 
421  DO ii = 1, 6
422  tot = 0.d0
423  DO jj = 1, 6
424  tot = tot + a2(ii,jj)*strain(jj)
425  ENDDO
426  strain_b(ii) = tot
427  ENDDO
428 
429  DO ii = 1, 6
430  tot = 0.d0
431  DO jj = 1, 6
432  tot = tot + dbd(ii,jj)*mu_d(jj)
433  ENDDO
434  strain_b(ii) = strain_b(ii) + tot
435  ENDDO
436 
437  DO ii = 1, 6
438  tot = 0.d0
439  DO jj = 1, 6
440  tot = tot + a2(ii,jj)*strain(jj)
441  ENDDO
442  strain_d(ii) = tot
443  ENDDO
444 
445  DO ii = 1, 6
446  tot = 0.d0
447  DO jj = 1, 6
448  tot = tot + ddd(ii,jj)*mu_d(jj)
449  ENDDO
450  strain_d(ii) = strain_d(ii) + tot
451  ENDDO
452 
453  DO ii = 1, 6
454  tot = 0.d0
455  DO jj = 1, 6
456  tot = tot + lm(ii,jj)*strain_m(jj)
457  ENDDO
458  stress_m(ii) = tot
459  ENDDO
460 
461  DO ii = 1, 6
462  tot = 0.d0
463  DO jj = 1, 6
464  tot = tot + l2(ii,jj)*strain_b(jj)
465  ENDDO
466  stress_b(ii) = tot
467  ENDDO
468 
469  DO ii = 1, 6
470  tot = 0.d0
471  DO jj = 1, 6
472  tot = tot + l2(ii,jj)*(strain_d(jj)-mu_d(jj))
473  ENDDO
474  stress_d(ii) = tot
475  ENDDO
476 
477 
478 ! calculate the stress -1
479 ! [S] = [C] {E}
480 !
481 
482 !!$ IF(ielem.EQ.MonNd)THEN
483 !!$
484 !!$ DO ii = 1, 6
485 !!$ PRINT*,llbar(ii,:)
486 !!$ ENDDO
487 !!$ PRINT*,'**',cb,cd,cm
488 !!$
489 !!$ ENDIF
490 
491  llbar = cm*matmul(lm,scriptam) + cb*matmul(l2,scriptab) + cd*matmul(l2,scriptad)
492 
493  s11 = llbar(1,1)*strain(1) + llbar(1,2)*strain(2) + llbar(1,3)*strain(3)
494  s22 = llbar(2,1)*strain(1) + llbar(2,2)*strain(2) + llbar(2,3)*strain(3)
495  s33 = llbar(3,1)*strain(1) + llbar(3,2)*strain(2) + llbar(3,3)*strain(3)
496  s12 = strain(6)*llbar(4,4)
497  s23 = strain(4)*llbar(5,5)
498  s13 = strain(5)*llbar(6,6)
499 
500  IF(abs(s33 - (cm*stress_m(3) + cb*stress_b(3) + cd*stress_d(3))).GT.0.0001d0)THEN
501  print*,cd, s33, cm*stress_m(3) + cb*stress_b(3) + cd*stress_d(3)
502  ENDIF
503 
504 
505 END SUBROUTINE matous_const_model
506 
Color::r
unsigned char r() const
Definition: Color.h:68
y
void int int REAL REAL * y
Definition: read.cpp:74
vecmatvecmul
subroutine vecmatvecmul(a, B, c, ndim, scaler)
Definition: A_D_tensors.f90:556
p1
NT p1
Definition: rational_rotation.h:108
sqrt
double sqrt(double d)
Definition: double.h:73
operatorj
subroutine operatorj(Lm, Ld, Lb, L2, Lbar, Mm, M2, cm, cb, cd, Dbm, Dbb, Dbd, Dmm, Dmb, Dmd, Ddm, Ddb, Ddd, Am, Ab, Ad, A2, Km, K2, Gm, G2, vm, v2, a_eta, a_zeta, OpJ, N, Ninv, R, Rinv, Q2, Qm, W)
Definition: OperatorJ.f90:53
a_d_tensors
subroutine a_d_tensors(Lo, Lm, L2, Mm, M2, cm, cb, cd, Dbm, Dbb, Dbd, Dmm, Dmb, Dmd, Ddm, Ddb, Ddd, A2, Am, L_bar)
Definition: A_D_tensors.f90:53
matous_const_model
subroutine matous_const_model(cd, cb, p1, p2, Yin, SoftParam, Km, K2, Gm, G2, vm, nu2, a_eta, a_zeta, cm, c2, Lo, Lm, L2, Mm, M2, L_bar, strain, DelStrain, s11, s22, s33, s12, s13, s23, ielem)
Definition: Matous_Const_Model.f90:53
test
void test(void)
Definition: flotsam.C:99
INode::pt
Parent_type pt(const Point_2 &p)
Definition: HDS_overlay.h:365
n
const NT & n
Definition: rational_rotation.h:72
Aff_transformation_repS2::print
virtual std::ostream & print(std::ostream &os) const
Definition: Aff_transformationS2.h:277