1! CalculiX - A 3-dimensional finite element program 2! Copyright (C) 1998-2021 Guido Dhondt 3! 4! This program is free software; you can redistribute it and/or 5! modify it under the terms of the GNU General Public License as 6! published by the Free Software Foundation(version 2); 7! 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12! GNU General Public License for more details. 13! 14! You should have received a copy of the GNU General Public License 15! along with this program; if not, write to the Free Software 16! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 17! 18 subroutine materialdata_cp_sec(imat,ntmat_,t1l,shcon,nshcon,cp, 19 & physcon) 20! 21 implicit none 22! 23! determines the secant specific heat at constant pressure cp 24! 25! the difference with materialdata_cp is that the specific heat at 26! constant pressure cp as returned from the present routine 27! is the secant value and not the differential value. 28! For the differential value we have: 29! dh=cp*dT 30! and consequently 31! h=int_from_0_to_T cp*dT 32! For the secant value one has: 33! h=cp_secant*T 34! 35 integer imat,ntmat_,id,nshcon(*),four,i 36! 37 real*8 t1l,shcon(0:3,ntmat_,*),cp,physcon(*) 38! 39 four=4 40! 41! calculating the tangent specific heat 42! 43 call ident2(shcon(0,1,imat),t1l,nshcon(imat),four,id) 44 if(nshcon(imat).eq.0) then 45 continue 46 elseif(nshcon(imat).eq.1) then 47 cp=shcon(1,1,imat) 48 elseif(id.eq.0) then 49 cp=shcon(1,1,imat) 50 elseif(id.eq.nshcon(imat)) then 51 cp=(shcon(0,1,imat)-physcon(1))*shcon(1,1,imat) 52 do i=2,nshcon(imat) 53 cp=cp+(shcon(0,i,imat)-shcon(0,i-1,imat))* 54 & (shcon(1,i,imat)+shcon(1,i-1,imat))/2.d0 55 enddo 56 cp=cp+(t1l-shcon(0,nshcon(imat),imat))* 57 & (shcon(1,nshcon(imat),imat)) 58 cp=cp/(t1l-physcon(1)) 59 else 60 cp=shcon(1,id,imat)+ 61 & (shcon(1,id+1,imat)-shcon(1,id,imat))* 62 & (t1l-shcon(0,id,imat))/ 63 & (shcon(0,id+1,imat)-shcon(0,id,imat)) 64 cp=(t1l-shcon(0,id,imat))*(cp+shcon(1,id,imat))/2.d0 65 do i=2,id 66 cp=cp+(shcon(0,i,imat)-shcon(0,i-1,imat))* 67 & (shcon(1,i,imat)+shcon(1,i-1,imat))/2.d0 68 enddo 69 cp=cp+(shcon(0,1,imat)-physcon(1))*shcon(1,1,imat) 70 cp=cp/(t1l-physcon(1)) 71 endif 72! 73 return 74 end 75 76 77 78 79 80 81 82