FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
8contains
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
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
320end 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_get_coefficient(tpoi, imat, coef, ntab, temp, funca, funcb)
subroutine heat_capacity_shell_741(etype, nn, ecoord, temperature, imat, thick, 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_c1(etype, nn, ecoord, temperature, imat, surf, lumped, mass, ntab1, temp1, funca1, funcb1, ntab2, temp2, funca2, funcb2)
subroutine heat_capacity_c3(etype, nn, ecoord, temperature, imat, 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