FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
heat_LIB_CAPACITY.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 !-------------------------------------------------------------------------------
7 
8 contains
9 
10  subroutine heat_get_coefficient(Tpoi, imat, coef, ntab, temp, funcA, funcB)
11  use hecmw
12  implicit none
13  real(kind=kreal) :: coef(3)
14  integer(kind=kint) :: i, in, imat, itab, ntab
15  real(kind=kreal) :: tpoi, temp(ntab), funca(ntab+1), funcb(ntab+1)
16 
17  itab = 0
18  if(tpoi < temp(1)) then
19  itab = 1
20  elseif(temp(ntab) <= tpoi)then
21  itab = ntab + 1
22  else
23  do in = 1, ntab - 1
24  if(temp(in) <= tpoi .and. tpoi < temp(in+1))then
25  itab = in + 1
26  exit
27  endif
28  enddo
29  endif
30 
31  coef(1) = funca(itab)*tpoi + funcb(itab)
32  coef(2) = coef(1)
33  coef(3) = coef(1)
34  end subroutine heat_get_coefficient
35 
36  subroutine heat_capacity_c1(etype, nn, ecoord, temperature, IMAT, surf, lumped, mass, &
37  ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
38  use hecmw
39  implicit none
40  integer(kind=kint), intent(in) :: etype
41  integer(kind=kint), intent(in) :: nn
42  real(kind=kreal), intent(in) :: temperature(nn)
43  real(kind=kreal), intent(out) :: mass(:,:)
44  real(kind=kreal), intent(out) :: lumped(:)
45  integer(kind=kint), parameter :: ndof = 1
46  integer(kind=kint) :: i, j, IMAT
47  integer(kind=kint) :: ntab1, ntab2
48  real(kind=kreal) :: ecoord(3, nn)
49  real(kind=kreal) :: dx, dy, dz, surf, val, temp, length
50  real(kind=kreal) :: spe(3), den(3)
51  real(kind=kreal) :: temp1(ntab1), funca1(ntab1+1), funcb1(ntab1+1)
52  real(kind=kreal) :: temp2(ntab2), funca2(ntab2+1), funcb2(ntab2+1)
53 
54  if(etype == 611)then
55  ecoord(1,2) = ecoord(1,3); ecoord(2,2) = ecoord(2,3); ecoord(3,2) = ecoord(3,3)
56  endif
57 
58  dx = ecoord(1,2) - ecoord(1,1)
59  dy = ecoord(2,2) - ecoord(2,1)
60  dz = ecoord(3,2) - ecoord(3,1)
61  length = dsqrt(dx*dx + dy*dy + dz*dz)
62  val = surf*length
63 
64  temp = (temperature(1) + temperature(2))*0.5d0
65  call heat_get_coefficient(temp, imat, spe, ntab1, temp1, funca1, funcb1)
66  call heat_get_coefficient(temp, imat, den, ntab2, temp2, funca2, funcb2)
67 
68  lumped(1) = val*spe(1)*den(1)*0.5d0
69  lumped(2) = val*spe(1)*den(1)*0.5d0
70  end subroutine heat_capacity_c1
71 
72  subroutine heat_capacity_c2(etype, nn, ecoord, temperature, IMAT, thick, lumped, mass, &
73  ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
74  use mmechgauss
75  use m_matmatrix
76  use elementinfo
77  implicit none
78  !type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points
79  integer(kind=kint), intent(in) :: etype
80  integer(kind=kint), intent(in) :: nn
81  real(kind=kreal), intent(in) :: ecoord(2,nn)
82  real(kind=kreal), intent(out) :: mass(:,:)
83  real(kind=kreal), intent(out) :: lumped(:)
84  real(kind=kreal), intent(in) :: temperature(nn)
85  type(tmaterial), pointer :: matl
86  integer(kind=kint) :: i, j, lx, imat, ntab1, ntab2
87  real(kind=kreal) :: naturalcoord(2)
88  real(kind=kreal) :: func(nn), thick, temp
89  real(kind=kreal) :: det, wg, rho, diag_mass, total_mass
90  real(kind=kreal) :: d(1,1), n(1,nn), dn(1,nn)
91  real(kind=kreal) :: gderiv(nn,2)
92  real(kind=kreal) :: spe(3), den(3)
93  real(kind=kreal) :: temp1(ntab1), funca1(ntab1+1), funcb1(ntab1+1)
94  real(kind=kreal) :: temp2(ntab2), funca2(ntab2+1), funcb2(ntab2+1)
95  logical :: is_lumped
96 
97  mass = 0.0d0
98  lumped = 0.0d0
99  !matl => gausses(1)%pMaterial
100 
101  do lx = 1, numofquadpoints(etype)
102  call getquadpoint(etype, lx, naturalcoord)
103  call getshapefunc(etype, naturalcoord, func)
104  call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
105 
106  temp = dot_product(func, temperature)
107  call heat_get_coefficient(temp, imat, spe, ntab1, temp1, funca1, funcb1)
108  call heat_get_coefficient(temp, imat, den, ntab2, temp2, funca2, funcb2)
109 
110  d(1,1) = spe(1)*den(1)*thick
111  wg = getweight(etype, lx)*det
112 
113  n = 0.0d0
114  do i = 1, nn
115  n(1,i) = func(i)
116  enddo
117 
118  dn = matmul(d, n)
119  forall(i = 1:nn, j = 1:nn)
120  mass(i,j) = mass(i,j) + dot_product(n(:,i), dn(:,j))*wg
121  end forall
122  enddo
123 
124  is_lumped = .true.
125  if(is_lumped) call get_lumped_mass(nn, 1, mass, lumped)
126  end subroutine heat_capacity_c2
127 
128  subroutine heat_capacity_c3(etype, nn, ecoord, temperature, IMAT, lumped, mass, &
129  ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
130  use mmechgauss
131  use m_matmatrix
132  use elementinfo
133  implicit none
134  !type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points
135  integer(kind=kint), intent(in) :: etype
136  integer(kind=kint), intent(in) :: nn
137  real(kind=kreal), intent(in) :: ecoord(3,nn)
138  real(kind=kreal), intent(out) :: mass(:,:)
139  real(kind=kreal), intent(out) :: lumped(:)
140  real(kind=kreal), intent(in), optional :: temperature(nn)
141  type(tmaterial), pointer :: matl
142  integer(kind=kint) :: i, j, lx, imat, ntab1, ntab2
143  real(kind=kreal) :: naturalcoord(3)
144  real(kind=kreal) :: func(nn), temp
145  real(kind=kreal) :: det, wg, rho, diag_mass, total_mass
146  real(kind=kreal) :: d(1, 1), n(1, nn), dn(1, nn)
147  real(kind=kreal) :: gderiv(nn, 3)
148  real(kind=kreal) :: spe(3), den(3)
149  real(kind=kreal) :: temp1(ntab1), funca1(ntab1+1), funcb1(ntab1+1)
150  real(kind=kreal) :: temp2(ntab2), funca2(ntab2+1), funcb2(ntab2+1)
151  logical :: is_lumped
152 
153  mass = 0.0d0
154  lumped = 0.0d0
155  !matl => gausses(1)%pMaterial
156 
157  do lx = 1, numofquadpoints(etype)
158  call getquadpoint(etype, lx, naturalcoord)
159  call getshapefunc(etype, naturalcoord, func)
160  call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
161 
162  temp = dot_product(func, temperature)
163  call heat_get_coefficient(temp, imat, spe, ntab1, temp1, funca1, funcb1)
164  call heat_get_coefficient(temp, imat, den, ntab2, temp2, funca2, funcb2)
165 
166  d = 0.0d0
167  d(1,1) = spe(1)*den(1)
168  wg = getweight(etype, lx)*det
169 
170  n = 0.0d0
171  do i = 1, nn
172  n(1,i) = func(i)
173  enddo
174 
175  dn = matmul(d, n)
176  forall(i = 1:nn, j = 1:nn)
177  mass(i,j) = mass(i,j) + dot_product(n(:,i), dn(:,j))*wg
178  end forall
179  enddo
180 
181  is_lumped = .true.
182  if(is_lumped) call get_lumped_mass(nn, 1, mass, lumped)
183  end subroutine heat_capacity_c3
184 
185  subroutine heat_capacity_shell_731(etype, nn, ecoord, temperature, IMAT, thick, lumped, mass, &
186  ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
187  use mmechgauss
188  use m_matmatrix
189  use elementinfo
190  use m_dynamic_mass
191  implicit none
192  !type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points
193  integer(kind=kint), intent(in) :: etype
194  integer(kind=kint), intent(in) :: nn
195  real(kind=kreal), intent(in) :: ecoord(3,nn)
196  real(kind=kreal), intent(out) :: mass(:,:)
197  real(kind=kreal), intent(out) :: lumped(:)
198  real(kind=kreal), intent(in), optional :: temperature(nn)
199  type(tmaterial), pointer :: matl
200  integer(kind=kint) :: i, j, lx, imat, ntab1, ntab2
201  real(kind=kreal) :: surf, thick, temp
202  real(kind=kreal) :: det, wg, rho, diag_mass, total_mass
203  real(kind=kreal) :: spe(3), den(3)
204  real(kind=kreal) :: temp1(ntab1), funca1(ntab1+1), funcb1(ntab1+1)
205  real(kind=kreal) :: temp2(ntab2), funca2(ntab2+1), funcb2(ntab2+1)
206  logical :: is_lumped
207 
208  mass = 0.0d0
209  lumped = 0.0d0
210  !matl => gausses(1)%pMaterial
211 
212  surf = get_face3(ecoord)
213 
214  do i = 1, nn
215  temp = temperature(i)
216  call heat_get_coefficient(temp, imat, spe, ntab1, temp1, funca1, funcb1)
217  call heat_get_coefficient(temp, imat, den, ntab2, temp2, funca2, funcb2)
218 
219  mass(i,i) = mass(i,i) + surf*thick*spe(1)*den(1)/3.0d0
220  enddo
221 
222  is_lumped = .true.
223  if(is_lumped) call get_lumped_mass(nn, 1, mass, lumped)
224  end subroutine heat_capacity_shell_731
225 
226  subroutine heat_capacity_shell_741(etype, nn, ecoord, temperature, IMAT, thick, lumped, mass, &
227  ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
228  use mmechgauss
229  use m_matmatrix
230  use elementinfo
231  implicit none
232  !type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points
233  integer(kind=kint), intent(in) :: etype
234  integer(kind=kint), intent(in) :: nn
235  real(kind=kreal), intent(in) :: ecoord(3,nn)
236  real(kind=kreal), intent(out) :: mass(:,:)
237  real(kind=kreal), intent(out) :: lumped(:)
238  real(kind=kreal), intent(in), optional :: temperature(nn)
239  type(tmaterial), pointer :: matl
240  integer(kind=kint) :: i, j, lx, ly, imat, ntab1, ntab2
241  real(kind=kreal) :: naturalcoord(2)
242  real(kind=kreal) :: func(nn), thick, temp
243  real(kind=kreal) :: det, wg, rho, diag_mass, total_mass
244  real(kind=kreal) :: d(1,1), n(1, nn), dn(1, nn)
245  real(kind=kreal) :: gderiv(nn,2)
246  real(kind=kreal) :: spe(3), den(3)
247  real(kind=kreal) :: temp1(ntab1), funca1(ntab1+1), funcb1(ntab1+1)
248  real(kind=kreal) :: temp2(ntab2), funca2(ntab2+1), funcb2(ntab2+1)
249  real(kind=kreal) :: xg(2), ri, si, rp, sp, rm, sm, hr(4), hs(4)
250  real(kind=kreal) :: xr, xs, yr, ys, zr, zs
251  real(kind=kreal) :: h(4), x(4), y(4), z(4)
252  logical :: is_lumped
253 
254  mass = 0.0d0
255  lumped = 0.0d0
256  !matl => gausses(1)%pMaterial
257 
258  x(1) = ecoord(1,1); y(1) = ecoord(2,1); z(1) = ecoord(3,1)
259  x(2) = ecoord(1,2); y(2) = ecoord(2,2); z(2) = ecoord(3,2)
260  x(3) = ecoord(1,3); y(3) = ecoord(2,3); z(3) = ecoord(3,3)
261  x(4) = ecoord(1,4); y(4) = ecoord(2,4); z(4) = ecoord(3,4)
262 
263  xg(1) = -0.5773502691896258d0
264  xg(2) = -xg(1)
265 
266  do lx = 1, 2
267  ri = xg(lx)
268  do ly = 1, 2
269  si = xg(ly)
270  rp = 1.0d0 + ri
271  sp = 1.0d0 + si
272  rm = 1.0d0 - ri
273  sm = 1.0d0 - si
274 
275  !C*INTERPOLATION FUNCTION
276  h(1) = 0.25d0*rp*sp
277  h(2) = 0.25d0*rm*sp
278  h(3) = 0.25d0*rm*sm
279  h(4) = 0.25d0*rp*sm
280 
281  !C* FOR R-COORDINATE
282  hr(1) = 0.25d0*sp
283  hr(2) = -0.25d0*sp
284  hr(3) = -0.25d0*sm
285  hr(4) = 0.25d0*sm
286 
287  !C* FOR S-COORDINATE
288  hs(1) = 0.25d0*rp
289  hs(2) = 0.25d0*rm
290  hs(3) = -0.25d0*rm
291  hs(4) = -0.25d0*rp
292 
293  !C*JACOBI MATRIX
294  xr = hr(1)*x(1) + hr(2)*x(2) + hr(3)*x(3) + hr(4)*x(4)
295  xs = hs(1)*x(1) + hs(2)*x(2) + hs(3)*x(3) + hs(4)*x(4)
296  yr = hr(1)*y(1) + hr(2)*y(2) + hr(3)*y(3) + hr(4)*y(4)
297  ys = hs(1)*y(1) + hs(2)*y(2) + hs(3)*y(3) + hs(4)*y(4)
298  zr = hr(1)*z(1) + hr(2)*z(2) + hr(3)*z(3) + hr(4)*z(4)
299  zs = hs(1)*z(1) + hs(2)*z(2) + hs(3)*z(3) + hs(4)*z(4)
300 
301  det = (yr*zs - zr*ys)**2 + (zr*xs - xr*zs)**2 + (xr*ys - yr*xs)**2
302  det = dsqrt(det)
303 
304  temp = 0.0d0
305  do i = 1, 4
306  temp = temp + temperature(i)*h(i)
307  enddo
308  call heat_get_coefficient(temp, imat, spe, ntab1, temp1, funca1, funcb1)
309  call heat_get_coefficient(temp, imat, den, ntab2, temp2, funca2, funcb2)
310 
311  do i = 1, 4
312  mass(i,i) = mass(i,i) + spe(1)*den(1)*h(i)*det*thick
313  enddo
314  enddo
315  enddo
316 
317  is_lumped = .true.
318  if(is_lumped) call get_lumped_mass(nn, 1, mass, lumped)
319  end subroutine heat_capacity_shell_741
320 end module m_heat_lib_capacity
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
Definition: hecmw.f90:6
This module contains subroutines used in 3d eigen analysis for.
Definition: dynamic_mass.f90:6
real(kind=kreal) function get_face3(ecoord)
subroutine get_lumped_mass(nn, ndof, mass, lumped)
subroutine heat_capacity_c2(etype, nn, ecoord, temperature, IMAT, thick, lumped, mass, ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
subroutine heat_capacity_c1(etype, nn, ecoord, temperature, IMAT, surf, lumped, mass, ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
subroutine heat_get_coefficient(Tpoi, imat, coef, ntab, temp, funcA, funcB)
subroutine heat_capacity_c3(etype, nn, ecoord, temperature, IMAT, lumped, mass, ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
subroutine heat_capacity_shell_731(etype, nn, ecoord, temperature, IMAT, thick, lumped, mass, ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
subroutine heat_capacity_shell_741(etype, nn, ecoord, temperature, IMAT, thick, lumped, mass, ntab1, temp1, funcA1, funcB1, ntab2, temp2, funcA2, funcB2)
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6