Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cst_coh.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 cst_coh(coor,lmc,matc,R_co,d,deltan,deltat, &
54  sigmax,taumax,sinit,sthresh,nstep,numnp,numco,numat_coh, &
55  delta, numploadelem, idpressload, pressload, &
56  numboundmesh, idmesh, rmesh, nboundtype,num_fail_coh, &
57  index_fail_coh, npress, nmesh, ielem_coh, iface_coh)
58 
59  IMPLICIT NONE
60 
61 !-----global variables
62  INTEGER :: nstep ! time step
63  INTEGER :: idstrt ! time step to start output
64  INTEGER :: idint ! time interval for output
65  INTEGER :: numnp ! number of nodal points
66  INTEGER :: numco ! number of CST cohesive elements
67  INTEGER :: numat_coh ! number of cohesive materials
68  INTEGER :: start_ccst ! number of cohesive elements
69  ! on border
70  INTEGER :: iflag
71 
72 !-----S/F interaction variables (ALE explicit code)
73  INTEGER :: npress
74  INTEGER :: nmesh
75  INTEGER :: numploadelem
76  INTEGER :: numboundmesh
77  INTEGER :: num_fail_coh
78  INTEGER, DIMENSION(1:5,1:npress) :: idpressload
79  REAL*8 , DIMENSION(1:3,1:npress) :: pressload
80  INTEGER, DIMENSION(1:4,1:nmesh) :: idmesh
81  REAL*8 , DIMENSION(1:3,1:nmesh) :: rmesh
82  INTEGER, DIMENSION(1:numnp) :: nboundtype
83  INTEGER, DIMENSION(1:numco) :: index_fail_coh
84  INTEGER, DIMENSION(1:2,1:numco) :: ielem_coh
85  INTEGER, DIMENSION(1:2,1:numco) :: iface_coh
86 
87 ! -- global coordinates
88  REAL*8, DIMENSION(1:3,1:numnp) :: coor
89 ! -- connectivity table for cohesive el
90  INTEGER, DIMENSION(1:6,1:numco) :: lmc
91 ! -- mat number for cohesive element
92  INTEGER, DIMENSION(1:numco) :: matc
93 ! -- cohesive traction loads
94  real*8,DIMENSION(1:3*numnp) :: r_co
95 ! -- nodal displacements
96  REAL*8, DIMENSION(1:3*numnp) :: d
97 ! -- the threshold, then damage of edge
98  REAL*8, DIMENSION(1:3,1:numco) :: sthresh
99 ! -- normal characteristic length
100  REAL*8, DIMENSION(1:numat_coh) :: deltan
101 ! -- tangential characteristic length
102  REAL*8, DIMENSION(1:numat_coh) :: deltat
103 ! -- maximum normal stress
104  REAL*8, DIMENSION(1:numat_coh) :: sigmax
105 ! -- maximum normal stress
106  REAL*8, DIMENSION(1:numat_coh) :: taumax
107 ! -- initial Sthresh
108  REAL*8, DIMENSION(1:numat_coh) :: sinit
109  REAL*8 :: delta ! time step
110 !-----local variables
111  REAL*8 :: deltn, deltt ! the char length of current element
112  REAL*8 :: sig ! max normal stress of element
113  REAL*8 :: tau ! max shearing stress of element
114  REAL*8 :: area ! area of face
115  REAL*8 :: dx1,dy1,dz1 ! nd 1 openings in the x,y,z direction
116  REAL*8 :: dx2,dy2,dz2 ! nd 2 openings in the x,y,z direction
117  REAL*8 :: dx3,dy3,dz3 ! nd 3 openings in the x,y,z direction
118  REAL*8 :: dn1,dn2,dn3 ! openings in the normal direction
119  REAL*8 :: dt1,dt2,dt3 ! openings in the tan direction
120  REAL*8 :: dn,dt ! opening at the sample point
121  REAL*8 :: tn,tt ! cohesive traction at sample point
122  REAL*8 :: tx1,tx2,tx3 ! cohesive traction at sample point
123  REAL*8 :: ty1,ty2,ty3 ! cohesive traction at sample point
124  REAL*8 :: tz1,tz2,tz3 ! cohesive traction at sample point
125  REAL*8 :: rx1,rx2,rx3 ! cohesive forces at nodes 1,2,3
126  REAL*8 :: ry1,ry2,ry3 ! cohesive forces at nodes 1,2,3
127  REAL*8 :: rz1,rz2,rz3 ! cohesive forces at nodes 1,2,3
128 ! -- openings in tagent directions
129  REAL*8 :: dt_x, dt_y, dt_z
130 ! -- Components of the tangent vector
131  REAL*8 :: xtangent_1,ytangent_1, ztangent_1
132  REAL*8 :: xtangent_2,ytangent_2, ztangent_2
133  REAL*8 :: xtangent_3,ytangent_3, ztangent_3
134 ! -- Components of two vectors parallel to plane
135  REAL*8 :: vec1_x,vec1_y,vec1_z,vec2_x,vec2_y,vec2_z
136 ! -- Componets of the normal vector of the face
137  REAL*8 :: xnormal,ynormal,znormal
138 ! -- Norm of the normal
139  REAL*8 :: xnorm
140 ! --
141 
142  REAL*8 :: g1,g2 ! gauss quadrature points
143  REAL*8 :: w1 ! weights
144  INTEGER :: m ! material number of current element
145  INTEGER :: n1,n2,n3,n4,n5,n6 ! 3 times the node number of element
146  INTEGER :: i ! loop counter
147  INTEGER :: j ! loop counter
148  INTEGER :: iaux
149  REAL*8 :: x ! dummy
150  parameter(g1 = 0.666666666666667, &
151  g2 = 0.166666666666667, &
152  w1 = 0.333333333333333)
153 !-----The cohesive law
154  DO i = 1,numco
155  IF (index_fail_coh(i) .eq. 0) then
156 
157  m = matc(i)
158  deltn = 1.d0/deltan(m)
159  deltt = 1.d0/deltat(m)
160  sig = sigmax(m)
161  tau = taumax(m)
162  n1 = lmc(1,i)*3 ! n1 = node number 1 (*3)
163  n2 = lmc(2,i)*3 ! n2 = node number 2 (*3)
164  n3 = lmc(3,i)*3 ! n3 = node number 3 (*3)
165  n4 = lmc(4,i)*3 ! n4 = node number 4 (*3)
166  n5 = lmc(5,i)*3 ! n5 = node number 5 (*3)
167  n6 = lmc(6,i)*3 ! n6 = node number 6 (*3)
168 
169 ! --- Calculate the normal
170 !
171 ! Components of two vectors parallel to plane
172  vec1_x = coor(1,lmc(3,i)) - coor(1,lmc(1,i))
173  vec1_y = coor(2,lmc(3,i)) - coor(2,lmc(1,i))
174  vec1_z = coor(3,lmc(3,i)) - coor(3,lmc(1,i))
175  vec2_x = coor(1,lmc(2,i)) - coor(1,lmc(1,i))
176  vec2_y = coor(2,lmc(2,i)) - coor(2,lmc(1,i))
177  vec2_z = coor(3,lmc(2,i)) - coor(3,lmc(1,i))
178 ! Take the cross product
179  xnormal = vec1_y*vec2_z - vec2_y*vec1_z
180  ynormal = - (vec1_x*vec2_z - vec2_x*vec1_z)
181  znormal = vec1_x*vec2_y - vec2_x*vec1_y
182 ! Length of the normal
183  xnorm = sqrt(xnormal**2 + ynormal**2 + znormal**2)
184 
185 ! Normalizing the normal
186 ! _ _
187 ! _ Vec1 x Vec2
188 ! n = ------------------
189 ! || Vec1 x Vec2 ||
190 !
191  xnormal = xnormal/xnorm
192  ynormal = ynormal/xnorm
193  znormal = znormal/xnorm
194 
195  dx1 = d(n5-2) - d(n1-2)
196  dy1 = d(n5-1) - d(n1-1)
197  dz1 = d(n5) - d(n1)
198 !
199 ! _ _
200 ! COD = COD dot n
201 ! n
202 !
203  dn1 = dx1*xnormal + dy1*ynormal + dz1*znormal
204 !
205 ! _ _ _
206 ! COD = COD - COD n
207 ! t n
208 !
209  dt_x = dx1 - dn1*xnormal
210  dt_y = dy1 - dn1*ynormal
211  dt_z = dz1 - dn1*znormal
212 !
213 ! _
214 ! COD = || COD ||
215 ! t t
216 !
217  dt1 = sqrt( dt_x**2 + dt_y**2 + dt_z**2)
218 ! _
219 ! COD
220 ! t
221 ! Tangent = t = ----
222 ! COD
223 ! t
224 !
225  IF(dt1.EQ.0.d0)THEN
226  xtangent_1 = 0.d0
227  ytangent_1 = 0.d0
228  ztangent_1 = 0.d0
229  ELSE
230  xtangent_1 = dt_x / dt1
231  ytangent_1 = dt_y / dt1
232  ztangent_1 = dt_z / dt1
233  ENDIF
234 
235  dt1 = dt1*deltt !
236  dn1 = dn1*deltn ! normalized openings
237 
238  dx2 = d(n6-2) - d(n2-2)
239  dy2 = d(n6-1) - d(n2-1)
240  dz2 = d(n6) - d(n2)
241 
242  dn2 = dx2*xnormal + dy2*ynormal + dz2*znormal
243 
244  dt_x = dx2 - dn2*xnormal
245  dt_y = dy2 - dn2*ynormal
246  dt_z = dz2 - dn2*znormal
247  dt2 = sqrt( dt_x**2 + dt_y**2 + dt_z**2)
248 
249  IF(dt2.EQ.0.d0)THEN
250  xtangent_2 = 0.d0
251  ytangent_2 = 0.d0
252  ztangent_2 = 0.d0
253  ELSE
254  xtangent_2 = dt_x / dt2
255  ytangent_2 = dt_y / dt2
256  ztangent_2 = dt_z / dt2
257  ENDIF
258 
259  dt2 = dt2*deltt
260  dn2 = dn2*deltn
261 
262  dx3 = d(n4-2) - d(n3-2)
263  dy3 = d(n4-1) - d(n3-1)
264  dz3 = d(n4) - d(n3)
265  dn3 = dx3*xnormal + dy3*ynormal + dz3*znormal
266 
267  dt_x = dx3 - dn3*xnormal
268  dt_y = dy3 - dn3*ynormal
269  dt_z = dz3 - dn3*znormal
270  dt3 = sqrt( dt_x**2 + dt_y**2 + dt_z**2)
271 
272  IF(dt3.EQ.0.d0)THEN
273  xtangent_3 = 0.d0
274  ytangent_3 = 0.d0
275  ztangent_3 = 0.d0
276  ELSE
277  xtangent_3 = dt_x / dt3
278  ytangent_3 = dt_y / dt3
279  ztangent_3 = dt_z / dt3
280  ENDIF
281 
282  dt3 = dt3*deltt
283  dn3 = dn3*deltn
284 
285 !----gauss point 1
286  dt = g1*dt1 + g2*dt2 + g2*dt3
287  dn = g1*dn1 + g2*dn2 + g2*dn3
288  x = 1.d0 - sqrt( dn*dn + dt*dt )
289  IF (x.LE.0.d0) THEN
290  sthresh(1,i) = 0.d0
291  ELSEIF (x.LE.sthresh(1,i)) THEN
292  sthresh(1,i) = x
293  ENDIF
294  IF (dn.GT.0.d0) THEN
295  tn = sthresh(1,i)/(1.d0-sthresh(1,i))*sig*dn
296  ELSE
297  tn = sinit(m)/(1.d0-sinit(m))*sig*dn
298  ENDIF
299  tt = sthresh(1,i)/(1.d0-sthresh(1,i))*tau*dt
300 ! _ _ _
301 ! T = T n + T t
302 ! n t
303 !
304  tx1 = tn*xnormal + tt*xtangent_1
305  ty1 = tn*ynormal + tt*ytangent_1
306  tz1 = tn*znormal + tt*ztangent_1
307 !----gauss point 2
308  dt = g2*dt1 + g1*dt2 + g2*dt3
309  dn = g2*dn1 + g1*dn2 + g2*dn3
310  x = 1.d0 - sqrt( dn*dn + dt*dt )
311  IF (x.LE.0.d0) THEN
312  sthresh(2,i) = 0.d0
313  ELSEIF (x.LE.sthresh(2,i)) THEN
314  sthresh(2,i) = x
315  ENDIF
316  IF (dn.GT.0.d0) THEN
317  tn = sthresh(2,i)/(1.d0-sthresh(2,i))*sig*dn
318  ELSE
319  tn = sinit(m)/(1.d0-sinit(m))*sig*dn
320  ENDIF
321  tt = sthresh(2,i)/(1.d0-sthresh(2,i))*tau*dt
322 ! _ _ _
323 ! T = T n + T t
324 ! n t
325 !
326  tx2 = tn*xnormal + tt*xtangent_2
327  ty2 = tn*ynormal + tt*ytangent_2
328  tz2 = tn*znormal + tt*ztangent_2
329 !----gauss point 3
330  dt = g2*dt1 + g2*dt2 + g1*dt3
331  dn = g2*dn1 + g2*dn2 + g1*dn3
332  x = 1.d0 - sqrt( dn*dn + dt*dt )
333  IF (x.LE.0.d0) THEN
334  sthresh(3,i) = 0.d0
335  ELSEIF (x.LE.sthresh(3,i)) THEN
336  sthresh(3,i) = x
337  ENDIF
338  IF (dn.GT.0.d0) THEN
339  tn = sthresh(3,i)/(1.d0-sthresh(3,i))*sig*dn
340  ELSE
341  tn = sinit(m)/(1.d0-sinit(m))*sig*dn
342  ENDIF
343  tt = sthresh(3,i)/(1.d0-sthresh(3,i))*tau*dt
344 
345 ! _ _ _
346 ! T = T n + T t
347 ! n t
348 !
349  tx3 = tn*xnormal + tt*xtangent_3
350  ty3 = tn*ynormal + tt*ytangent_3
351  tz3 = tn*znormal + tt*ztangent_3
352 !
353 ! Area of the triangle is half the area of the parallelogram
354 ! determined by the cross product of the vectors P12 and P13
355 
356  area = 0.5d0*xnorm
357 
358  rx1 = area*w1*(g1*tx1+g2*tx2+g2*tx3)
359  ry1 = area*w1*(g1*ty1+g2*ty2+g2*ty3)
360  rz1 = area*w1*(g1*tz1+g2*tz2+g2*tz3)
361  rx2 = area*w1*(g2*tx1+g1*tx2+g2*tx3)
362  ry2 = area*w1*(g2*ty1+g1*ty2+g2*ty3)
363  rz2 = area*w1*(g2*tz1+g1*tz2+g2*tz3)
364  rx3 = area*w1*(g2*tx1+g2*tx2+g1*tx3)
365  ry3 = area*w1*(g2*ty1+g2*ty2+g1*ty3)
366  rz3 = area*w1*(g2*tz1+g2*tz2+g1*tz3)
367 
368  r_co(n1-2) = r_co(n1-2) + rx1
369  r_co(n1-1) = r_co(n1-1) + ry1
370  r_co(n1) = r_co(n1) + rz1
371  r_co(n5-2) = r_co(n5-2) - rx1
372  r_co(n5-1) = r_co(n5-1) - ry1
373  r_co(n5) = r_co(n5) - rz1
374  r_co(n2-2) = r_co(n2-2) + rx2
375  r_co(n2-1) = r_co(n2-1) + ry2
376  r_co(n2) = r_co(n2) + rz2
377  r_co(n6-2) = r_co(n6-2) - rx2
378  r_co(n6-1) = r_co(n6-1) - ry2
379  r_co(n6) = r_co(n6) - rz2
380  r_co(n3-2) = r_co(n3-2) + rx3
381  r_co(n3-1) = r_co(n3-1) + ry3
382  r_co(n3) = r_co(n3) + rz3
383  r_co(n4-2) = r_co(n4-2) - rx3
384  r_co(n4-1) = r_co(n4-1) - ry3
385  r_co(n4) = r_co(n4) - rz3
386 !
387 ! find out free surface as a result of crack propagation
388 !
389 !------ Update mesh BC conditions
390 
391  IF ((sthresh(1,i) .le. 1.0e-3) .AND.&
392  (sthresh(2,i) .le. 1.0e-3) .AND.&
393  (sthresh(3,i) .le. 1.0e-3) ) Then
394 
395  IF (index_fail_coh(i) .eq. 0) then
396  num_fail_coh = num_fail_coh + 1
397  index_fail_coh(i) = 1
398 
399  print *,' cohesive element ',i,' failed'
400  print *,' total number of failed coh elems ',&
401  num_fail_coh
402 
403 !------------ Solid Propellant part
404  numploadelem = numploadelem + 1
405  idpressload(1,numploadelem) = ielem_coh(1,i)
406  idpressload(2,numploadelem) = iface_coh(1,i)
407  idpressload(3,numploadelem) = lmc(1,i)
408  idpressload(4,numploadelem) = lmc(2,i)
409  idpressload(5,numploadelem) = lmc(3,i)
410  pressload(1,numploadelem) = 1.0d0
411  pressload(2,numploadelem) = 1.0d0
412  pressload(3,numploadelem) = 1.0d0
413 
414 !------------ Case part
415  numploadelem = numploadelem + 1
416 !-- be careful of node numbering !
417  idpressload(1,numploadelem) = ielem_coh(2,i)
418  idpressload(2,numploadelem) = iface_coh(2,i)
419  idpressload(3,numploadelem) = lmc(5,i)
420  idpressload(4,numploadelem) = lmc(4,i)
421  idpressload(5,numploadelem) = lmc(6,i)
422  pressload(1,numploadelem) = 1.0d0
423  pressload(2,numploadelem) = 1.0d0
424  pressload(3,numploadelem) = 1.0d0
425 
426 !------------ store newly-generated crack face subjected to burning
427 
428  do j = 1, 3 ! only for SP
429  numboundmesh = numboundmesh + 1
430  idmesh(1,numboundmesh) = lmc(j,i)
431  if(nboundtype(lmc(j,i)) .eq. 8) then ! corner node
432  idmesh(2,numboundmesh) = 0
433  idmesh(3,numboundmesh) = 0
434  idmesh(4,numboundmesh) = 1
435  rmesh( 1,numboundmesh) = 1.0d0
436  rmesh( 2,numboundmesh) = -1.0d0
437  rmesh( 3,numboundmesh) = 0.0d0
438  else
439  idmesh(2,numboundmesh) = 1
440  idmesh(3,numboundmesh) = 0
441  idmesh(4,numboundmesh) = 1
442  rmesh( 1,numboundmesh) = 0.0d0
443  rmesh( 2,numboundmesh) = -1.0d0
444  rmesh( 3,numboundmesh) = 0.0d0
445  endif
446  enddo
447 
448  ENDIF
449  ENDIF
450 
451  ENDIF
452  ENDDO
453 
454  print *,' total number of failed coh elems ',num_fail_coh
455 
456  RETURN
457  END
458 
Aff_transformation_2::m
FT m(int i, int j) const
Definition: Aff_transformation_2.h:228
d
const NT & d
Definition: rational_rotation.h:73
sqrt
double sqrt(double d)
Definition: double.h:73
cst_coh
subroutine cst_coh(coor, lmc, matc, R_co, d, deltan, deltat, sigmax, taumax, Sinit, Sthresh, nstep, numnp, numco, numat_coh, delta, numploadelem, idpressload, pressload, numboundmesh, idmesh, rmesh, nboundtype, num_fail_coh, index_fail_coh, npress, nmesh, ielem_coh, iface_coh)
Definition: cst_coh.f90:53
i
blockLoc i
Definition: read.cpp:79
x
void int int REAL * x
Definition: read.cpp:74
Direction_2::delta
RT delta(int i) const
Definition: Direction_2.h:159
Aff_transformation_repS2::print
virtual std::ostream & print(std::ostream &os) const
Definition: Aff_transformationS2.h:277
j
j indices j
Definition: Indexing.h:6