OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
gruneisen.F File Reference
#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine gruneisen (iflag, nel, pm, off, eint, mu, mu2, espe, dvol, df, vnew, mat, rho0, pnew, dpdm, dpde, psh, nummat, npropm)

Function/Subroutine Documentation

◆ gruneisen()

subroutine gruneisen ( integer, intent(in) iflag,
integer, intent(in) nel,
dimension(npropm,nummat), intent(in) pm,
dimension(nel), intent(in) off,
dimension(nel), intent(inout) eint,
dimension(nel), intent(in) mu,
dimension(nel), intent(in) mu2,
dimension(nel), intent(in) espe,
dimension(nel), intent(in) dvol,
dimension(nel), intent(in) df,
dimension(nel), intent(in) vnew,
integer, dimension(nel), intent(in) mat,
dimension(nel), intent(in) rho0,
dimension(nel), intent(inout) pnew,
dimension(nel), intent(inout) dpdm,
dimension(nel), intent(inout) dpde,
dimension(nel), intent(inout) psh,
integer, intent(in) nummat,
integer, intent(in) npropm )
Parameters
[in]nummatnumber of material laws (size for PM array)
[in]npropmsize for PM array
[in]iflagflag for the operation (staggered scheme IFLAG=0,1 ; collocated scheme IFLAG=2)
[in]nelnumber of elements in the current group
[in]matmaterial identifier of each element

Definition at line 28 of file gruneisen.F.

32C-----------------------------------------------
33C D e s c r i p t i o n
34C-----------------------------------------------
35!----------------------------------------------------------------------------
36!! \details STAGGERED SCHEME IS EXECUTED IN TWO PASSES IN EOSMAIN : IFLG=0 THEN IFLG=1
37!! \details COLLOCATED SCHEME IS DOING A SINGLE PASS : IFLG=2
38!! \details
39!! \details STAGGERED SCHEME
40!! \details EOSMAIN / IFLG = 0 : DERIVATIVE CALCULATION FOR SOUND SPEED ESTIMATION c[n+1] REQUIRED FOR PSEUDO-VISCOSITY (DPDE:partial derivative, DPDM:total derivative)
41!! \details MQVISCB : PSEUDO-VISCOSITY Q[n+1]
42!! \details MEINT : INTERNAL ENERGY INTEGRATION FOR E[n+1] : FIRST PART USING P[n], Q[n], and Q[n+1] CONTRIBUTIONS
43!! \details EOSMAIN / IFLG = 1 : UPDATE P[n+1], T[N+1]
44!! \details INTERNAL ENERGY INTEGRATION FOR E[n+1] : LAST PART USING P[n+1] CONTRIBUTION
45!! \details (second order integration dE = -P.dV where P = 0.5(P[n+1] + P[n]) )
46!! \details COLLOCATED SCHEME
47!! \details EOSMAIN / IFLG = 2 : SINGLE PASS FOR P[n+1] AND DERIVATIVES
48!----------------------------------------------------------------------------
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56C
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER,INTENT(IN) :: NUMMAT !< number of material laws (size for PM array)
61 INTEGER,INTENT(IN) :: NPROPM !< size for PM array
62 INTEGER,INTENT(IN) :: IFLAG !< flag for the operation (staggered scheme IFLAG=0,1 ; collocated scheme IFLAG=2)
63 INTEGER,INTENT(IN) :: NEL !< number of elements in the current group
64 INTEGER,INTENT(IN) :: MAT(NEL) !< material identifier of each element
65 my_real,INTENT(IN) :: pm(npropm,nummat) !< material buffer
66 my_real,INTENT(IN) :: off(nel) !< OFF=1 if element is active, OFF=0 if element is inactive
67 my_real,INTENT(IN) :: mu(nel) !< volumetric strain : rho/rho0-1
68 my_real,INTENT(IN) :: mu2(nel) !< MU**2 IF MU > 0 , 0 otherwise
69 my_real,INTENT(IN) :: espe(nel) !< specific internal energy
70 my_real,INTENT(IN) :: dvol(nel) !< volume change
71 my_real,INTENT(IN) :: vnew(nel) !< current volume
72 my_real,INTENT(IN) :: rho0(nel) !< reference density
73 my_real,INTENT(IN) :: df(nel) !< relative volume (v/v0 = rho0/rho)
74 my_real,INTENT(INOUT) :: eint(nel) !< internal energy
75 my_real,INTENT(INOUT) :: psh(nel) !< pressure shift (for relative pressure modeling)
76 my_real,INTENT(INOUT) :: pnew(nel) !< current pressure
77 my_real,INTENT(INOUT) :: dpde(nel) !< partial derivative dP/dE where E=Eint/V0
78 my_real,INTENT(INOUT) :: dpdm(nel) !< total derivative : DP/Dmu = dPdmu + dPdE* P/(1+mu)**2
79C-----------------------------------------------
80C L o c a l V a r i a b l e s
81C-----------------------------------------------
82 INTEGER I, MX
83 my_real aa, bb, dvv, ff, fg, fac, xx, dff, dfg, fac1, pp
84 my_real cc(nel),s1(nel),s2(nel),s3(nel),g0(nel),ga(nel),pc(nel)
85C-----------------------------------------------
86C B o d y
87C-----------------------------------------------
88 IF(iflag == 0) THEN !--- SOUND SPEED DERIVATIVE (DPDE is partial derivative, DPDM is total derivative)
89 ! user parameters
90 DO i=1,nel
91 mx = mat(i)
92 cc(i) = pm(33,mx)
93 s1(i) = pm(34,mx)
94 s2(i) = pm(160,mx)
95 s3(i) = pm(161,mx)
96 g0(i) = pm(35,mx)
97 ga(i) = pm(36,mx)
98 pc(i) = pm(37,mx)
99 psh(i) = pm(88,mx)
100 ENDDO
101
102 ! Derivative
103 DO i=1,nel
104 fac=one
105 fac1=one
106 IF(mu(i) > zero) THEN
107 xx= mu(i)/(one+mu(i))
108 ff=one+(one-half*g0(i))*mu(i)-half*ga(i)*mu2(i)
109 fg=one-(s1(i)-one+s2(i)*xx+s3(i)*xx*xx)*mu(i)
110 fac=ff/(fg*fg)
111 dff=one-half*g0(i)-ga(i)*mu(i)
112 dfg=one-s1(i)+xx*(-two*s2(i)+xx*(s2(i)-three*s3(i))+two*s3(i)*xx*xx)
113 fac1=fac*(one+mu(i)*(dff/ff-two*dfg/fg))
114 ENDIF
115 aa=fac*rho0(i)*cc(i)*cc(i)*mu(i)
116 bb=g0(i)+ga(i)*mu(i)
117 pp=max(aa+bb*espe(i),pc(i))*off(i) ! total pressure need for sound speed calculation (=> no pressure shift)
118 dpdm(i)=fac1*rho0(i)*cc(i)*cc(i)+pp*df(i)*df(i)*bb+ga(i)*espe(i)
119 dpde(i)=bb
120 pnew(i) = max(pp,pc(i))*off(i)! P(mu[n+1],E[n])
121 pnew(i) = pnew(i) - psh(i)
122 ENDDO
123
124 ELSEIF(iflag == 1) THEN !--- P[n+1] and EINT[n+1]
125 ! user parameters
126 DO i=1,nel
127 mx = mat(i)
128 cc(i) = pm(33,mx)
129 s1(i) = pm(34,mx)
130 s2(i) = pm(160,mx)
131 s3(i) = pm(161,mx)
132 g0(i) = pm(35,mx)
133 ga(i) = pm(36,mx)
134 pc(i) = pm(37,mx)
135 psh(i) = pm(88,mx)
136 ENDDO
137
138 ! Pressure P[n+1]
139 ! Internal Energy Eint[n+1] (2nd order integration requires P[n+1])
140 DO i=1,nel
141 fac=one
142 IF(mu(i) > zero) THEN
143 xx= mu(i)/(one+mu(i))
144 ff=one+(one-half*g0(i))*mu(i)-half*ga(i)*mu2(i)
145 fg=one-(s1(i)-one+s2(i)*xx+s3(i)*xx*xx)*mu(i)
146 fac=ff/(fg*fg)
147 ENDIF
148 aa=fac*rho0(i)*cc(i)*cc(i)*mu(i)
149 bb=g0(i)+ga(i)*mu(i)
150 dpde(i) = bb
151 dvv=half*dvol(i)*df(i) / max(em15,vnew(i))
152 pnew(i)=(aa+bb*espe(i))/(one+bb*dvv)
153 pnew(i)= max(pnew(i),pc(i))*off(i) ! P(mu[n+1],E[n+1])
154 eint(i)=eint(i) - half*dvol(i)*pnew(i) !total pressure for energy integration
155 pnew(i) = pnew(i) - psh(i)
156 ENDDO
157
158 ELSEIF(iflag == 2) THEN !--- collocated scheme (law151)
159 ! user parameters
160 DO i=1, nel
161 mx = mat(i)
162 cc(i) = pm(33,mx)
163 s1(i) = pm(34,mx)
164 s2(i) = pm(160,mx)
165 s3(i) = pm(161,mx)
166 g0(i) = pm(35,mx)
167 ga(i) = pm(36,mx)
168 pc(i) = pm(37,mx)
169 psh(i) = pm(88,mx)
170 ENDDO
171 ! pressure and sound speed derivatives
172 DO i=1,nel
173 IF (vnew(i) > zero) THEN
174 fac=one
175 fac1=one
176 IF(mu(i) > zero) THEN
177 xx= mu(i)/(one+mu(i))
178 ff=one+(one-half*g0(i))*mu(i)-half*ga(i)*mu2(i)
179 fg=one-(s1(i)-one+s2(i)*xx+s3(i)*xx*xx)*mu(i)
180 fac=ff/(fg*fg)
181 dff=one-half*g0(i)-ga(i)*mu(i)
182 dfg=one-s1(i)+xx*(-two*s2(i)+xx*(s2(i)-three*s3(i))+two*s3(i)*xx*xx)
183 fac1=fac*(one+mu(i)*(dff/ff-two*dfg/fg))
184 ENDIF
185 aa=fac*rho0(i)*cc(i)*cc(i)*mu(i)
186 bb=g0(i)+ga(i)*mu(i)
187 pnew(i)=max(aa+bb*espe(i),pc(i))*off(i)
188 dpdm(i)=fac1*rho0(i)*cc(i)*cc(i)+pnew(i)*df(i)*df(i)*bb+ga(i)*espe(i)
189 dpde(i)=bb
190 pnew(i)=pnew(i)-psh(i)
191 ENDIF
192 ENDDO
193
194 ENDIF
195
196 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21