FrontISTR  6.0.0
Large-scale structural analysis program with finit element method
ElasticLinear.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
6 module m_elasticlinear
7  use hecmw_util
8  use mmaterial
9 
10  implicit none
11 
12 contains
13 
16  subroutine calelasticmatrix( matl, sectType, D, temp, hdflag )
17  type( tmaterial ), intent(in) :: matl
18  integer, intent(in) :: sectType
19  real(kind=kreal), intent(out) :: d(:,:)
20  real(kind=kreal), intent(in) :: temp
21  real(kind=kreal) :: ee, pp, coef1, coef2, ina(1), outa(2)
22  integer(kind=kint), intent(in), optional :: hdflag
23  logical :: ierr
24  real(kind=kreal) :: k, g
25 
26  d(:,:)=0.d0
27 
28  ina(1) = temp
29  call fetch_tabledata( mc_isoelastic, matl%dict, outa, ierr, ina )
30  if( ierr ) then
31  ee = matl%variables(m_youngs)
32  pp = matl%variables(m_poisson)
33  else
34  ee = outa(1)
35  pp = outa(2)
36  endif
37 
38  select case (secttype)
39  case (d3)
40  k = ee/(1.d0-2.d0*pp)/3.d0
41  g = ee/(1.d0+pp)*0.5d0
42  if( present(hdflag) ) then
43  if( hdflag == 1 ) k = 0.d0
44  if( hdflag == 2 ) g = 0.d0
45  end if
46  d(1,1)=k+(4.d0/3.d0)*g
47  d(1,2)=k-(2.d0/3.d0)*g
48  d(4,4)=g
49  d(1,3)=d(1,2)
50  d(2,1)=d(1,2)
51  d(2,2)=d(1,1)
52  d(2,3)=d(1,2)
53  d(3,1)=d(1,3)
54  d(3,2)=d(2,3)
55  d(3,3)=d(1,1)
56  d(5,5)=d(4,4)
57  d(6,6)=d(4,4)
58  case (planestress)
59  coef1=ee/(1.d0-pp*pp)
60  coef2=0.5d0*(1.d0-pp)
61  d(1,1)=coef1
62  d(1,2)=coef1*pp
63  d(1,3)=0.d0
64  d(2,1)=d(1,2)
65  d(2,2)=d(1,1)
66  d(2,3)=0.d0
67  d(3,1)=0.d0
68  d(3,2)=0.d0
69  d(3,3)=coef1*coef2
70  case (planestrain)
71  coef1=ee/((1.d0+pp)*(1.d0-2.d0*pp))
72  coef2=ee/(2.d0*(1.d0+pp))
73  d(1,1)=coef1*(1.d0-pp)
74  d(1,2)=coef1*pp
75  d(1,3)=0.d0
76  d(2,1)=d(1,2)
77  d(2,2)=d(1,1)
78  d(2,3)=0.d0
79  d(3,1)=0.d0
80  d(3,2)=0.d0
81  d(3,3)=coef2
82  case (axissymetric)
83  coef1=ee*(1.d0-pp)/((1.d0+pp)*(1.d0-2.d0*pp))
84  coef2=(1.d0-2.d0*pp)/(2.d0*(1.d0-pp))
85  d(1,1)=coef1
86  d(1,2)=coef1*pp/(1.d0-pp)
87  d(1,3)=0.d0
88  d(1,4)=d(1,2)
89  d(2,1)=d(1,2)
90  d(2,2)=d(1,1)
91  d(2,3)=0.d0
92  d(2,4)=d(1,2)
93  d(3,1)=0.d0
94  d(3,2)=0.d0
95  d(3,3)=coef1*coef2
96  d(3,4)=0.d0
97  d(4,1)=d(1,4)
98  d(4,2)=d(2,4)
99  d(4,3)=0.d0
100  d(4,4)=d(1,1)
101  case default
102  stop "Section type not defined"
103  end select
104 
105  end subroutine
106 
107 
112  subroutine calelasticmatrix_ortho( matl, sectType, bij, DMAT, temp, hdflag )
113  use m_utilities
114  type( tmaterial ), intent(in) :: matl
115  integer, intent(in) :: secttype
116  real(kind=kreal), intent(in) :: bij(3,3)
117  real(kind=kreal), intent(out) :: dmat(:,:)
118  real(kind=kreal), intent(in) :: temp
119  integer(kind=kint), intent(in), optional :: hdflag
120  real(kind=kreal) :: e1, e2, e3, g12, g23, g13, nyu12, nyu23,nyu13
121  real(kind=kreal) :: nyu21,nyu32,nyu31, delta1, ina(1), outa(9)
122  real(kind=kreal) :: tm(6,6)
123  real(kind=kreal) :: cvol, dvol(6,6), amat(6,6), uvec(6)
124  integer :: ii, jj
125  logical :: ierr
126 
127  ina(1) = temp
128  call fetch_tabledata( mc_orthoelastic, matl%dict, outa, ierr, ina )
129  if( ierr ) then
130  stop "Fails in fetching orthotropic elastic constants!"
131  endif
132 
133  e1 = outa(1)
134  e2 = outa(2)
135  e3 = outa(3)
136  nyu12 = outa(4)
137  nyu13 = outa(5)
138  nyu23 = outa(6)
139  g12 = outa(7)
140  g13 = outa(8)
141  g23 = outa(9)
142  nyu21 = e2/e1*nyu12
143  nyu32 = e3/e2*nyu23
144  nyu31 = e3/e1*nyu13
145  delta1 = 1.d0/(1.d0 -nyu12*nyu21 -nyu23*nyu32 -nyu31*nyu13 -2.d0*nyu21*nyu32*nyu13)
146 
147  dmat(:,:)=0.d0
148  dmat(1,1) = e1*(1.d0-nyu23*nyu32)*delta1
149  dmat(2,2) = e2*(1.d0-nyu13*nyu31)*delta1
150  dmat(3,3) = e3*(1.d0-nyu12*nyu21)*delta1
151  dmat(1,2) = e1*(nyu21+nyu31*nyu23)*delta1
152  dmat(1,3) = e1*(nyu31+nyu21*nyu32)*delta1
153  dmat(2,3) = e2*(nyu32+nyu12*nyu31)*delta1
154  dmat(4,4) = g12
155  dmat(5,5) = g23
156  dmat(6,6) = g13
157 
158  dmat(2,1) = dmat(1,2)
159  dmat(3,2) = dmat(2,3)
160  dmat(3,1) = dmat(1,3)
161 
162  call transformation(bij, tm)
163 
164  dmat = matmul( transpose(tm), dmat)
165  dmat = matmul( dmat, (tm) )
166 
167  ! Projector split for selective ES/NS: c = 1/(m^T D^{-1} m), D_vol = c m m^T
168  if( present(hdflag) ) then
169  if( hdflag == 1 .or. hdflag == 2 ) then
170  amat(:,:) = dmat(:,:)
171  uvec(:) = 0.d0; uvec(1) = 1.d0; uvec(2) = 1.d0; uvec(3) = 1.d0
172  call calsolve(6, amat, uvec)
173  cvol = 1.d0 / (uvec(1) + uvec(2) + uvec(3))
174  dvol(:,:) = 0.d0
175  do ii = 1, 3
176  do jj = 1, 3
177  dvol(ii, jj) = cvol
178  end do
179  end do
180  if( hdflag == 1 ) then
181  dmat(:,:) = dmat(:,:) - dvol(:,:)
182  else if( hdflag == 2 ) then
183  dmat(:,:) = dvol(:,:)
184  end if
185  end if
186  end if
187 
188  end subroutine
189 
190 
191  !####################################################################
192  subroutine linearelastic_shell &
193  (matl, secttype, c, &
194  e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
195  alpha, n_layer, cdsys)
196  !####################################################################
197 
198  type( tmaterial ), intent(in) :: matl
199  integer, intent(in) :: sectType, n_layer
200  real(kind = kreal), intent(out) :: c(:, :, :, :)
201  real(kind = kreal), intent(in) :: e1_hat(3), e2_hat(3), e3_hat(3)
202  real(kind = kreal), intent(in) :: cg1(3), cg2(3), cg3(3)
203  real(kind = kreal), intent(out) :: alpha
204  real(kind = kreal), intent(in), optional :: cdsys(3, 3)
205 
206  !--------------------------------------------------------------------
207 
208  real(kind = kreal) :: ee, pp, ee2, g12, g23, g31, theta, pp2
209  real(kind = kreal) :: outa(2)
210  real(kind = kreal) :: lambda1, lambda2, mu, k_correction, r44, r55
211  real(kind = kreal) :: c_hat(3, 3, 3, 3), d(5,5), d_hat(5,5), d_temp(5,5), t(5,5)
212  real(kind = kreal) :: e_hat_dot_cg(3, 3)
213  real(kind = kreal) :: alpha_over_mu
214  real(kind = kreal) :: v_ref(3), v_proj(3), v_proj_abs
215  real(kind = kreal) :: v_ref_y(3), v_sense
216  real(kind = kreal) :: theta_ref, phi, c_ref, s_ref, c_phi, s_phi
217  real(kind = kreal) :: d44_ref, d45_ref, d55_ref
218 
219  integer :: index_i(5), index_j(5), index_k(5), index_l(5)
220  integer :: is, js, ii, ij, ik, il
221  integer :: i, j, k, l, total_layer, matl_type
222 
223  logical :: ierr
224 
225  !--------------------------------------------------------------------
226 
227  call fetch_tabledata(mc_isoelastic, matl%dict, outa, ierr)
228 
229  !--------------------------------------------------------------------
230 
231  matl_type = matl%shell_var(n_layer)%ortho
232  if(matl_type == 0)then
233  if( ierr ) then
234 
235  ee = matl%shell_var(n_layer)%ee
236  pp = matl%shell_var(n_layer)%pp
237  alpha_over_mu = matl%variables(m_alpha_over_mu)
238 
239  else
240 
241  ee = outa(1)
242  pp = outa(2)
243  alpha_over_mu = matl%variables(m_alpha_over_mu)
244 
245  end if
246 
247  !--------------------------------------------------------------------
248 
249  ! Elastic constant
250  lambda1 = ee/( 1.0d0-pp*pp )
251  lambda2 = pp*lambda1
252  mu = 0.5d0*ee/( 1.0d0+pp )
253 
254  !--------------------------------------------------------------------
255 
256  ! Shear correction factor
257  k_correction = 5.0d0/6.0d0
258  !k_correction = 1.0D0
259 
260  ! Laminate transverse shear ratio (layers share ee and pp when table data is given)
261  r44 = 1.0d0
262  r55 = 1.0d0
263  if( ierr ) call shelltransshearratio( matl, r44, r55 )
264 
265  !--------------------------------------------------------------------
266 
267  alpha = alpha_over_mu*mu
268 
269  !--------------------------------------------------------------------
270 
271  ! Constitutive tensor
272  c_hat(:, :, :, :) = 0.0d0
273 
274  !--------------------------------------------------------
275 
276  c_hat(1, 1, 1, 1) = lambda1
277  c_hat(1, 1, 2, 2) = lambda2
278  c_hat(2, 2, 1, 1) = lambda2
279  c_hat(2, 2, 2, 2) = lambda1
280  c_hat(1, 2, 1, 2) = mu
281  c_hat(1, 2, 2, 1) = mu
282  c_hat(2, 1, 1, 2) = mu
283  c_hat(2, 1, 2, 1) = mu
284  c_hat(1, 3, 1, 3) = k_correction*mu*r55
285  c_hat(1, 3, 3, 1) = k_correction*mu*r55
286  c_hat(2, 3, 2, 3) = k_correction*mu*r44
287  c_hat(2, 3, 3, 2) = k_correction*mu*r44
288  c_hat(3, 1, 3, 1) = k_correction*mu*r55
289  c_hat(3, 1, 1, 3) = k_correction*mu*r55
290  c_hat(3, 2, 3, 2) = k_correction*mu*r44
291  c_hat(3, 2, 2, 3) = k_correction*mu*r44
292 
293  !--------------------------------------------------------
294 
295  elseif(matl_type == 1)then
296  total_layer = matl%totallyr
297 
298  ee = matl%shell_var(n_layer)%ee
299  pp = matl%shell_var(n_layer)%pp
300  ee2 = matl%shell_var(n_layer)%ee2
301  g12 = matl%shell_var(n_layer)%g12
302  g23 = matl%shell_var(n_layer)%g23
303  g31 = matl%shell_var(n_layer)%g31
304  theta = matl%shell_var(n_layer)%angle
305 
306  v_ref = 0.0d0
307  v_ref_y = 0.0d0
308  if( present(cdsys) ) v_ref(1:3) = cdsys(1, 1:3)
309  if( present(cdsys) ) v_ref_y(1:3) = cdsys(2, 1:3)
310  v_proj(1:3) = v_ref(1:3) - dot_product(v_ref, e3_hat) * e3_hat(1:3)
311  v_proj_abs = dsqrt( dot_product(v_proj, v_proj) )
312  v_sense = dot_product(v_proj, e1_hat) * dot_product(v_ref_y, e2_hat) &
313  - dot_product(v_proj, e2_hat) * dot_product(v_ref_y, e1_hat)
314  ! theta_ref: layer angle from the projected orientation axis
315  ! phi: angle from e1_hat to the projected orientation axis
316  theta_ref = sign(1.0d0, v_sense) * theta
317  phi = datan2( dot_product(v_proj, e2_hat), dot_product(v_proj, e1_hat) )
318  theta = theta_ref + phi
319 
320  alpha_over_mu = matl%variables(m_alpha_over_mu)
321 
322  !--------------------------------------------------------------------
323 
324  ! Shear correction factor
325  k_correction = 5.0d0/6.0d0
326 
327  ! Laminate transverse shear ratio, defined in the orientation frame
328  call shelltransshearratio( matl, r44, r55 )
329 
330  !--------------------------------------------------------------------
331 
332  alpha = alpha_over_mu * 0.5d0 * ee / ( 1.0d0+pp )
333 
334  !--------------------------------------------------------------------
335 
336  ! write(*,*) ee,pp,ee2,g12,g23,g31,theta,alpha_over_mu,n_layer
337 
338  d(:,:) = 0.0d0
339  d_hat(:,:) = 0.0d0
340  d_temp(:,:) = 0.0d0
341  t(:,:) = 0.0d0
342 
343  pp2 = pp * ee2 / ee
344 
345  d(1,1) = ee / (1.0d0 - pp * pp2)
346  d(1,2) = pp2 * ee / (1.0d0- pp * pp2)
347  d(2,1) = pp2 * ee / (1.0d0- pp * pp2)
348  d(2,2) = ee2 / (1.0d0 - pp * pp2)
349  d(3,3) = g12
350  d(4,4) = k_correction * g23
351  d(5,5) = k_correction * g31
352 
353  t(1,1) = cos(theta) * cos(theta)
354  t(1,2) = sin(theta) * sin(theta)
355  t(2,1) = sin(theta) * sin(theta)
356  t(2,2) = cos(theta) * cos(theta)
357  t(3,3) = cos(theta) * cos(theta) - sin(theta) * sin(theta)
358  t(1,3) = sin(theta) * cos(theta)
359  t(2,3) = -sin(theta) * cos(theta)
360  t(3,1) = -2.0d0 * sin(theta) * cos(theta)
361  t(3,2) = 2.0d0 * sin(theta) * cos(theta)
362  t(4,4) = cos(theta)
363  t(4,5) = -sin(theta)
364  t(5,4) = sin(theta)
365  t(5,5) = cos(theta)
366 
367  !--------------------- D_temp = [D]*[T]
368 
369  d_temp(1,1) = d(1,1)*t(1,1)+d(1,2)*t(2,1)
370  d_temp(1,2) = d(1,1)*t(1,2)+d(1,2)*t(2,2)
371  d_temp(2,1) = d(2,1)*t(1,1)+d(2,2)*t(2,1)
372  d_temp(2,2) = d(2,1)*t(1,2)+d(2,2)*t(2,2)
373  d_temp(3,1) = d(3,3)*t(3,1)
374  d_temp(3,2) = d(3,3)*t(3,2)
375  d_temp(1,3) = d(1,1)*t(1,3)+d(1,2)*t(2,3)
376  d_temp(2,3) = d(2,1)*t(1,3)+d(2,2)*t(2,3)
377  d_temp(3,3) = d(3,3)*t(3,3)
378  d_temp(4,4) = d(4,4)*t(4,4)
379  d_temp(4,5) = d(4,4)*t(4,5)
380  d_temp(5,4) = d(5,5)*t(5,4)
381  d_temp(5,5) = d(5,5)*t(5,5)
382 
383  !--------------------- D_hat = [trans_T]*[D_temp]
384 
385  d_hat(1,1) = t(1,1)*d_temp(1,1)+t(1,2)*d_temp(2,1)+t(3,1)*d_temp(3,1)
386  d_hat(1,2) = t(1,1)*d_temp(1,2)+t(1,2)*d_temp(2,2)+t(3,1)*d_temp(3,2)
387  d_hat(2,1) = t(2,1)*d_temp(1,1)+t(2,2)*d_temp(2,1)+t(3,2)*d_temp(3,1)
388  d_hat(2,2) = t(2,1)*d_temp(1,2)+t(2,2)*d_temp(2,2)+t(3,2)*d_temp(3,2)
389  d_hat(3,1) = t(1,3)*d_temp(1,1)+t(2,3)*d_temp(2,1)+t(3,3)*d_temp(3,1)
390  d_hat(3,2) = t(1,3)*d_temp(1,2)+t(2,3)*d_temp(2,2)+t(3,3)*d_temp(3,2)
391  d_hat(1,3) = t(1,1)*d_temp(1,3)+t(1,2)*d_temp(2,3)+t(3,1)*d_temp(3,3)
392  d_hat(2,3) = t(2,1)*d_temp(1,3)+t(2,2)*d_temp(2,3)+t(3,2)*d_temp(3,3)
393  d_hat(3,3) = t(1,3)*d_temp(1,3)+t(2,3)*d_temp(2,3)+t(3,3)*d_temp(3,3)
394  d_hat(4,4) = t(4,4)*d_temp(4,4)+t(5,4)*d_temp(5,4)
395  d_hat(4,5) = t(4,4)*d_temp(4,5)+t(5,4)*d_temp(5,5)
396  d_hat(5,4) = t(4,5)*d_temp(4,4)+t(5,5)*d_temp(5,4)
397  d_hat(5,5) = t(4,5)*d_temp(4,5)+t(5,5)*d_temp(5,5)
398 
399  ! Apply the correction in the orientation frame, preserving positive definiteness,
400  ! then rotate the transverse shear block by phi into the shell frame
401  c_ref = cos(theta_ref)
402  s_ref = sin(theta_ref)
403  d44_ref = r44 * ( c_ref * c_ref * d(4,4) + s_ref * s_ref * d(5,5) )
404  d45_ref = sqrt(r44 * r55) * c_ref * s_ref * ( d(5,5) - d(4,4) )
405  d55_ref = r55 * ( s_ref * s_ref * d(4,4) + c_ref * c_ref * d(5,5) )
406 
407  c_phi = cos(phi)
408  s_phi = sin(phi)
409  d_hat(4,4) = c_phi * c_phi * d44_ref + 2.0d0 * c_phi * s_phi * d45_ref + s_phi * s_phi * d55_ref
410  d_hat(4,5) = c_phi * s_phi * ( d55_ref - d44_ref ) + ( c_phi * c_phi - s_phi * s_phi ) * d45_ref
411  d_hat(5,4) = d_hat(4,5)
412  d_hat(5,5) = s_phi * s_phi * d44_ref - 2.0d0 * c_phi * s_phi * d45_ref + c_phi * c_phi * d55_ref
413 
414  !--------------------------------------------------------------------
415 
416  ! Constitutive tensor
417  c_hat(:, :, :, :) = 0.0d0
418  !write(*,*) 'Elastic linear.f90 make c_hat'
419 
420  !-------D_hat to c_hat
421 
422  index_i(1) = 1
423  index_i(2) = 2
424  index_i(3) = 1
425  index_i(4) = 2
426  index_i(5) = 3
427 
428  index_j(1) = 1
429  index_j(2) = 2
430  index_j(3) = 2
431  index_j(4) = 3
432  index_j(5) = 1
433 
434  index_k(1) = 1
435  index_k(2) = 2
436  index_k(3) = 1
437  index_k(4) = 2
438  index_k(5) = 3
439 
440  index_l(1) = 1
441  index_l(2) = 2
442  index_l(3) = 2
443  index_l(4) = 3
444  index_l(5) = 1
445 
446  !--------------------------------------------------------------------
447 
448  do js = 1, 5
449  do is = 1, 5
450 
451  ii = index_i(is)
452  ij = index_j(is)
453  ik = index_k(js)
454  il = index_l(js)
455 
456  c_hat(ii, ij, ik, il) = d_hat(is, js)
457  c_hat(ij, ii, ik, il) = d_hat(is, js) ! minor symmetry i<->j (Bug #764/#765)
458  c_hat(ii, ij, il, ik) = d_hat(is, js) ! minor symmetry k<->l
459  c_hat(ij, ii, il, ik) = d_hat(is, js) ! minor symmetry i<->j & k<->l
460 
461  end do
462  end do
463 
464  !--------------------------------------------------------
465 
466  else
467  !$omp critical
468  write(*,*) 'shell matl type isnot collect'
469  !$omp end critical
470  stop
471  endif
472 
473  select case( secttype )
474  case( shell )
475 
476  e_hat_dot_cg(1, 1) &
477  = e1_hat(1)*cg1(1)+e1_hat(2)*cg1(2) &
478  +e1_hat(3)*cg1(3)
479  e_hat_dot_cg(2, 1) &
480  = e2_hat(1)*cg1(1)+e2_hat(2)*cg1(2) &
481  +e2_hat(3)*cg1(3)
482  e_hat_dot_cg(3, 1) = 0.0d0
483  e_hat_dot_cg(1, 2) &
484  = e1_hat(1)*cg2(1)+e1_hat(2)*cg2(2) &
485  +e1_hat(3)*cg2(3)
486  e_hat_dot_cg(2, 2) &
487  = e2_hat(1)*cg2(1)+e2_hat(2)*cg2(2) &
488  +e2_hat(3)*cg2(3)
489  e_hat_dot_cg(3, 2) = 0.0d0
490  e_hat_dot_cg(1, 3) &
491  = e1_hat(1)*cg3(1)+e1_hat(2)*cg3(2) &
492  +e1_hat(3)*cg3(3)
493  e_hat_dot_cg(2, 3) &
494  = e2_hat(1)*cg3(1)+e2_hat(2)*cg3(2) &
495  +e2_hat(3)*cg3(3)
496  e_hat_dot_cg(3, 3) &
497  = e3_hat(1)*cg3(1)+e3_hat(2)*cg3(2) &
498  +e3_hat(3)*cg3(3)
499 
500  !--------------------------------------------------------
501 
502  ! Constitutive tensor
503 
504  c(1, 1, 1, 1) = 0.0d0
505  c(2, 2, 1, 1) = 0.0d0
506  c(1, 2, 1, 1) = 0.0d0
507  c(2, 2, 2, 2) = 0.0d0
508  c(1, 2, 2, 2) = 0.0d0
509  c(1, 2, 1, 2) = 0.0d0
510  c(3, 1, 1, 1) = 0.0d0
511  c(3, 1, 2, 2) = 0.0d0
512  c(3, 1, 1, 2) = 0.0d0
513  c(2, 3, 1, 1) = 0.0d0
514  c(2, 3, 2, 2) = 0.0d0
515  c(2, 3, 1, 2) = 0.0d0
516  c(3, 1, 3, 1) = 0.0d0
517  c(3, 1, 2, 3) = 0.0d0
518  c(2, 3, 2, 3) = 0.0d0
519 
520  do l = 1, 2
521 
522  do k = 1, 2
523 
524  do j = 1, 2
525 
526  do i = 1, 2
527 
528  c(1, 1, 1, 1) &
529  = c(1, 1, 1, 1) &
530  +c_hat(i, j, k, l) &
531  *e_hat_dot_cg(i, 1)*e_hat_dot_cg(j ,1) &
532  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
533  c(2, 2, 1, 1) &
534  = c(2, 2, 1, 1) &
535  +c_hat(i, j, k, l) &
536  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 2) &
537  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
538  c(1, 2, 1, 1) &
539  = c(1, 2, 1, 1) &
540  +c_hat(i, j, k, l) &
541  *e_hat_dot_cg(i, 1)*e_hat_dot_cg(j, 2) &
542  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
543  c(2, 2, 2, 2) &
544  = c(2, 2, 2, 2) &
545  +c_hat(i, j, k, l) &
546  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 2) &
547  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 2)
548 
549  c(1, 2, 2, 2) &
550  = c(1, 2, 2, 2) &
551  +c_hat(i, j, k, l) &
552  *e_hat_dot_cg(i, 1)*e_hat_dot_cg(j, 2) &
553  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 2)
554  c(1, 2, 1, 2) &
555  = c(1, 2, 1, 2) &
556  +c_hat(i, j, k, l) &
557  *e_hat_dot_cg(i, 1)*e_hat_dot_cg(j, 2) &
558  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 2)
559 
560  end do
561 
562  do i = 1, 3
563 
564  c(3, 1, 1, 1) &
565  = c(3, 1, 1, 1) &
566  +c_hat(i, j, k, l) &
567  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
568  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
569  c(3, 1, 2, 2) &
570  = c(3, 1, 2, 2) &
571  +c_hat(i, j, k, l) &
572  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
573  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 2)
574  c(3, 1, 1, 2) &
575  = c(3, 1, 1, 2) &
576  +c_hat(i, j, k, l) &
577  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
578  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 2)
579 
580  end do
581 
582  end do
583 
584  do j = 1, 3
585 
586  do i = 1, 2
587 
588  c(2, 3, 1, 1) &
589  = c(2, 3, 1, 1) &
590  +c_hat(i, j, k, l) &
591  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 3) &
592  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
593  c(2, 3, 2, 2) &
594  = c(2, 3, 2, 2) &
595  +c_hat(i, j, k, l) &
596  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 3) &
597  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 2)
598  c(2, 3, 1, 2) &
599  = c(2, 3, 1, 2) &
600  +c_hat(i, j, k, l) &
601  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 3) &
602  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 2)
603 
604  end do
605 
606  end do
607 
608  end do
609 
610  do k = 1, 3
611 
612  do j = 1, 2
613 
614  do i = 1, 3
615 
616  c(3, 1, 3, 1) &
617  = c(3, 1, 3, 1) &
618  +c_hat(i, j, k, l) &
619  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
620  *e_hat_dot_cg(k, 3)*e_hat_dot_cg(l, 1)
621 
622  end do
623 
624  end do
625 
626  end do
627 
628  end do
629 
630  do l = 1, 3
631 
632  do k = 1, 2
633 
634  do j = 1, 2
635 
636  do i = 1, 3
637 
638  c(3, 1, 2, 3) &
639  = c(3, 1, 2, 3) &
640  +c_hat(i, j, k, l) &
641  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
642  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 3)
643 
644  end do
645 
646  end do
647 
648  do j = 1, 3
649 
650  do i = 1, 2
651 
652  c(2, 3, 2, 3) &
653  = c(2, 3, 2, 3) &
654  +c_hat(i, j, k, l) &
655  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 3) &
656  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 3)
657 
658  end do
659 
660  end do
661 
662  end do
663 
664  end do
665 
666  c(1, 1, 2, 2) = c(2, 2, 1, 1)
667  c(1, 1, 1, 2) = c(1, 2, 1, 1)
668  c(1, 1, 2, 3) = c(2, 3, 1, 1)
669  c(1, 1, 3, 1) = c(3, 1, 1, 1)
670  c(2, 2, 1, 2) = c(1, 2, 2, 2)
671  c(2, 2, 2, 3) = c(2, 3, 2, 2)
672  c(2, 2, 3, 1) = c(3, 1, 2, 2)
673  c(1, 2, 2, 3) = c(2, 3, 1, 2)
674  c(1, 2, 3, 1) = c(3, 1, 1, 2)
675  c(2, 3, 3, 1) = c(3, 1, 2, 3)
676 
677  !--------------------------------------------------------
678 
679  ! DO l = 1, 3
680  !
681  ! DO k = 1, 3
682  !
683  ! DO j = 1, 3
684  !
685  ! DO i = 1, 3
686  !
687  ! c(i, j, k, l) = 0.0D0
688  !
689  ! DO ll = 1, 3
690  !
691  ! DO kk = 1, 3
692  !
693  ! DO jj = 1, 3
694  !
695  ! DO ii = 1, 3
696  !
697  ! c(i, j, k, l) &
698  ! = c(i, j, k, l) &
699  ! +c_hat(ii, jj, kk, ll) &
700  ! *e_hat_dot_cg(ii, i)*e_hat_dot_cg(jj, j) &
701  ! *e_hat_dot_cg(kk, k)*e_hat_dot_cg(ll, l)
702  !
703  ! END DO
704  !
705  ! END DO
706  !
707  ! END DO
708  !
709  ! END DO
710  !
711  ! END DO
712  !
713  ! END DO
714  !
715  ! END DO
716  !
717  ! END DO
718 
719  !--------------------------------------------------------
720 
721  end select
722 
723  !--------------------------------------------------------------------
724 
725  return
726 
727  !####################################################################
728  end subroutine linearelastic_shell
729  !####################################################################
730  ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
731 
735  subroutine shelltransshearratio( matl, r44, r55 )
736  type( tmaterial ), intent(in) :: matl
737  real(kind = kreal), intent(out) :: r44, r55
738 
739  real(kind = kreal) :: e_layer(2, matl%totallyr), g_layer(2, matl%totallyr)
740  real(kind = kreal) :: ee, pp, pp2, k_correction, ratio(2)
741  real(kind = kreal) :: zb, zt, zn, z, ea, es, d_bend, g_bot, g_pt, gg_int, gt_sum
742  real(kind = kreal) :: gauss_xi(3), gauss_w(3)
743  real(kind = kreal) :: q11, q22, q12h, cs, sn
744  integer :: nlyr, idir, i, ip
745 
746  r44 = 1.0d0
747  r55 = 1.0d0
748  nlyr = matl%totallyr
749  if( nlyr <= 1 ) return
750 
751  k_correction = 5.0d0/6.0d0
752  gauss_xi = (/ -sqrt(0.6d0), 0.0d0, sqrt(0.6d0) /)
753  gauss_w = (/ 5.0d0/9.0d0, 8.0d0/9.0d0, 5.0d0/9.0d0 /)
754 
755  do i = 1, nlyr
756  ee = matl%shell_var(i)%ee
757  pp = matl%shell_var(i)%pp
758  if( matl%shell_var(i)%ortho == 1 ) then
759  pp2 = pp * matl%shell_var(i)%ee2 / ee
760  q11 = ee / (1.0d0 - pp * pp2)
761  q22 = matl%shell_var(i)%ee2 / (1.0d0 - pp * pp2)
762  q12h = 2.0d0 * (pp2 * q11 + 2.0d0 * matl%shell_var(i)%g12)
763  cs = cos( matl%shell_var(i)%angle )
764  sn = sin( matl%shell_var(i)%angle )
765  e_layer(1, i) = sn**4 * q11 + cs**2 * sn**2 * q12h + cs**4 * q22
766  e_layer(2, i) = cs**4 * q11 + cs**2 * sn**2 * q12h + sn**4 * q22
767  g_layer(1, i) = cs**2 * matl%shell_var(i)%g23 + sn**2 * matl%shell_var(i)%g31
768  g_layer(2, i) = sn**2 * matl%shell_var(i)%g23 + cs**2 * matl%shell_var(i)%g31
769  else
770  e_layer(:, i) = ee / (1.0d0 - pp * pp)
771  g_layer(:, i) = 0.5d0 * ee / (1.0d0 + pp)
772  endif
773  end do
774 
775  ! Layer 1 is at the bottom surface; weights are thickness fractions summing to 1
776  ratio(:) = 1.0d0
777  do idir = 1, 2
778  if( any( g_layer(idir, :) <= 0.0d0 ) ) cycle
779 
780  ea = 0.0d0
781  es = 0.0d0
782  zb = -0.5d0
783  do i = 1, nlyr
784  zt = zb + matl%shell_var(i)%weight
785  ea = ea + e_layer(idir, i) * (zt - zb)
786  es = es + e_layer(idir, i) * 0.5d0 * (zt * zt - zb * zb)
787  zb = zt
788  end do
789  if( ea <= 0.0d0 ) cycle
790  zn = es / ea
791 
792  d_bend = 0.0d0
793  gg_int = 0.0d0
794  gt_sum = 0.0d0
795  g_bot = 0.0d0
796  zb = -0.5d0
797  do i = 1, nlyr
798  zt = zb + matl%shell_var(i)%weight
799  d_bend = d_bend + e_layer(idir, i) * ((zt - zn)**3 - (zb - zn)**3) / 3.0d0
800  do ip = 1, 3
801  z = 0.5d0 * (zb + zt) + 0.5d0 * (zt - zb) * gauss_xi(ip)
802  g_pt = g_bot + 0.5d0 * e_layer(idir, i) * ((z - zn)**2 - (zb - zn)**2)
803  gg_int = gg_int + 0.5d0 * (zt - zb) * gauss_w(ip) * g_pt * g_pt / g_layer(idir, i)
804  end do
805  g_bot = g_bot + 0.5d0 * e_layer(idir, i) * ((zt - zn)**2 - (zb - zn)**2)
806  gt_sum = gt_sum + g_layer(idir, i) * (zt - zb)
807  zb = zt
808  end do
809  if( gg_int > 0.0d0 ) ratio(idir) = d_bend * d_bend / gg_int / (k_correction * gt_sum)
810  end do
811 
812  r44 = ratio(1)
813  r55 = ratio(2)
814  end subroutine shelltransshearratio
815 
816 
817 
818 end module m_elasticlinear
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
m_elasticlinear
This module provides functions for elastic material.
Definition: ElasticLinear.f90:6
m_elasticlinear::linearelastic_shell
subroutine linearelastic_shell(matl, sectType, c, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, alpha, n_layer, cdsys)
Definition: ElasticLinear.f90:196
m_elasticlinear::calelasticmatrix
subroutine calelasticmatrix(matl, sectType, D, temp, hdflag)
Calculate isotropic elastic matrix hdflag: 0=full, 1=deviatoric only (edge-smoothed,...
Definition: ElasticLinear.f90:17
m_elasticlinear::shelltransshearratio
subroutine shelltransshearratio(matl, r44, r55)
Ratio of the energy-equivalent transverse shear stiffness D**2/int(g**2/G dz) of a laminate,...
Definition: ElasticLinear.f90:736
m_elasticlinear::calelasticmatrix_ortho
subroutine calelasticmatrix_ortho(matl, sectType, bij, DMAT, temp, hdflag)
Calculate orthotropic elastic matrix bij: rotation matrix from material axes to global axes,...
Definition: ElasticLinear.f90:113
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
m_utilities::calsolve
subroutine calsolve(NN, A, b)
Solve A x = b for dense NN x NN system (Gaussian elimination with partial pivoting)
Definition: utilities.f90:345
m_utilities::transformation
subroutine transformation(jacob, tm)
Definition: utilities.f90:392
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
mmaterial::m_youngs
integer(kind=kint), parameter m_youngs
Definition: material.f90:95
mmaterial::planestress
integer(kind=kint), parameter planestress
Definition: material.f90:88
mmaterial::d3
integer(kind=kint), parameter d3
Definition: material.f90:87
mmaterial::planestrain
integer(kind=kint), parameter planestrain
Definition: material.f90:89
mmaterial::shell
integer(kind=kint), parameter shell
Definition: material.f90:91
mmaterial::m_poisson
integer(kind=kint), parameter m_poisson
Definition: material.f90:96
mmaterial::axissymetric
integer(kind=kint), parameter axissymetric
Definition: material.f90:90
mmaterial::mc_isoelastic
character(len=dict_key_length) mc_isoelastic
Definition: material.f90:142
mmaterial::mc_orthoelastic
character(len=dict_key_length) mc_orthoelastic
Definition: material.f90:143
mmaterial::m_alpha_over_mu
integer(kind=kint), parameter m_alpha_over_mu
Definition: material.f90:109
mmaterial::tmaterial
Structure to manage all material related data.
Definition: material.f90:169