OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
czfintn.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| czfintn1 ../engine/source/elements/shell/coquez/czfintn.F
25!||--- called by ------------------------------------------------------
26!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
27!|| czforc3_crk ../engine/source/elements/xfem/czforc3_crk.F
28!||====================================================================
29 SUBROUTINE czfintn1(JFT ,JLT ,THK ,C1 ,AA ,VHG ,
30 2 X13 ,X24 ,Y13 ,Y24 ,Z1 ,MX23 ,
31 3 MX13 ,MX34 ,MY13 ,MY23 ,MY34 ,VGLAS,
32 4 VSTRE ,MSTRE ,VF ,VM ,FAC ,A11 ,
33 5 A12 ,G ,SHF ,SIGY ,OFF ,FAC1 ,
34 6 RHO ,AREA ,DT1 ,EINT ,AMU ,VHI ,
35 7 NPT ,IPARTC,EVIS ,KFTS ,GSR ,NEL ,
36 8 A11SR ,A12SR ,NUSR ,SHFSR,MTN ,FAC58)
37C elasto-plastic hourglass stresses+linear damping
38C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
39C CALCUL VF(3,4),VM(2,4) : FORCES INTERNES GENERALISEES
40C ENTREES : THK,1/AREA,PME,X13,X24,Y13,Y24,MX13,MX34,MY13,MY34
41C DHG(6) : INCREMENT DE LA DEF HOURGLASS GENERALISEE
42C VGLAS(12): DEF HOURGLASS GENERALISEE
43C VGLAS(1-6)->ETA; VGLAS(7-12)->KSI;
44C FAC=Et/E
45C SORTIES : VF,VM
46C VF(3,4),VM(2,4) :
47C J=1,2 STOCK LA PARTIE ANTISYM (F(1)=-F(3),F(2)=-F(4))
48C J=3,4 STOCK LA PARTIE SYM (F(1)=F(3),F(2)=F(4))
49C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
50#include "implicit_f.inc"
51#include "param_c.inc"
52#include "mvsiz_p.inc"
53C-----------------------------------------------
54C D U M M Y A R G U M E N T S
55C-----------------------------------------------
56 INTEGER JFT,JLT,NPT,IPARTC(*),NEL,MTN
57 my_real
58 . THK(*) ,C1(*) ,AA(*) ,VHG(MVSIZ,6),
59 . X13(*) ,X24(*) ,Y13(*) ,Y24(*) ,Z1(*),
60 . MX13(*),MX23(*),MX34(*),MY13(*),MY23(*),MY34(*),
61 . VGLAS(NEL,12),VSTRE(NEL,5),MSTRE(NEL,3),VF(MVSIZ,3,4),VM(MVSIZ,2,4),
62 . FAC(MVSIZ,2),A11(*) ,A12(*) ,G(*),SHF(*),VHI(MVSIZ,3),SIGY(*),
63 . fac1(*) ,rho(*) ,area(*), dt1,off(*),eint(jlt,2),amu(*) ,
64 . evis(npsav,*),
65 . gsr(*), a11sr(*), a12sr(*), nusr(*), shfsr(*),fac58(mvsiz,2)
66 INTEGER KFTS
67C-----------------------------------------------
68C L O C A L V A R I A B L E S
69C-----------------------------------------------
70 INTEGER I,MX, IC, II, J, JST(MVSIZ+1)
71 my_real
72 . STIER,C2,C3,C4,C5,C7,C8,HSURA,SXX,SYY,SVM,SVMN,
73 . ss1,ss2,ss3,sc6,sc5,cxx_k,cyy_k,hxx,hyy,sf1,sf2,b13,b24,
74 . cxx,cyy,y13s,x13s,y34s6,y23s5,x23s5,x34s6,hxx_k,hyy_k,
75 . bxx,byy,bxx_k,byy_k,
76 . c1m,c2m,c1b,c2b,cmm,cnn,ufac,tol,coef,eh1,eh2,
77 . cnnx,cnny,cnnx_k,cnny_k,cmmx,cmmy,cmmx_k,cmmy_k,sxy0,mxy0,
78 . svmc,sigy2,svm0,svmc_k,a_bk,a_b,a_k,undouzsr,aux,
79 . ssv0,ssv1,ssv2,ssv3,sc6_v,sc5_v,cxx_v,cyy_v,cxy_v,cxz_v,cyz_v,
80 . hxx_v,hyy_v,hxy_v,ss1_v,ss2_v,ss3_v,sf1_v,sf2_v,hvl,
81 . emy,ecx,ecy,esy,coef1,coefh,fbend,fbend_v,c12,
82 . c6(mvsiz),dglas(mvsiz,12),
83 . esx(mvsiz),tesy(mvsiz),emx(mvsiz),dhg(mvsiz,6),etmp(mvsiz,2)
84C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
85 c7=four_over_3
86 coefh= zep999
87 coef= zep85
88 stier = fivep333
89 tol= em18
90 fbend_v =threep464
91C
92 IF (npt==0) THEN
93 coef1=sixteen
94 ELSE
95 coef1=twenty5
96 ENDIF
97 IF (npt==1) THEN
98 fbend =zero
99 fbend_v =zero
100 ELSE
101 fbend =one_over_12
102 ENDIF
103 undouzsr=sqrt(one_over_12)
104
105 DO i=jft,jlt
106 dhg(i,1)=vhg(i,1)*dt1
107 dhg(i,2)=vhg(i,2)*dt1
108 dhg(i,3)=vhg(i,3)*dt1
109 dhg(i,4)=vhg(i,4)*dt1
110 dhg(i,5)=vhg(i,5)*dt1
111 dhg(i,6)=vhg(i,6)*dt1
112 END DO
113 IF(mtn==58) THEN
114 IF (npt==1) THEN
115 DO i=jft,jlt
116 c3=four*aa(i)
117 hxx=c3*my34(i)
118 hyy=c3*mx34(i)
119 hxx_k=c3*my23(i)
120 hyy_k=c3*mx23(i)
121 cxx=hxx*vglas(i,11)
122 cyy=hyy*vglas(i,12)
123 cxx_k=hxx_k*vglas(i,11)
124 cyy_k=hyy_k*vglas(i,12)
125 c2 = a11(i)*fac1(i)
126 c1m=c2*fac58(i,1)
127 c2m=c2*fac58(i,2)
128 c12 = em02*min(c1m,c2m)
129 dglas(i,1) =c1m*cxx-c12*cyy
130 dglas(i,2) =c2m*cyy-c12*cxx
131 dglas(i,7) =c1m*cxx_k-c12*cyy_k
132 dglas(i,8) =c2m*cyy_k-c12*cxx_k
133 END DO
134c if (mod(ncycle,10)==0) print *,'FAC58(I,1-2)=',FAC58(nel,1),FAC58(nel,2)
135 ELSE
136 DO i=jft,jlt
137 c3=four*aa(i)
138 hxx=c3*my34(i)
139 hyy=c3*mx34(i)
140 hxx_k=c3*my23(i)
141 hyy_k=c3*mx23(i)
142 cxx=hxx*dhg(i,1)
143 cyy=hyy*dhg(i,2)
144 cxx_k=hxx_k*dhg(i,1)
145 cyy_k=hyy_k*dhg(i,2)
146 bxx=hxx*dhg(i,3)
147 byy=hyy*dhg(i,4)
148 bxx_k=hxx_k*dhg(i,3)
149 byy_k=hyy_k*dhg(i,4)
150 c2 = a11(i)*fac1(i)
151 c1m=c2*fac58(i,1)
152 c2m=c2*fac58(i,2)
153 c12 = em02*min(c1m,c2m)
154C -----------for energy calculation-Mean_value--------
155 c6(i)=c1(i)*fbend
156 ss1= my34(i)*vglas(i,1) +my23(i)*vglas(i,7)
157 ss2= mx23(i)*vglas(i,8) +mx34(i)*vglas(i,2)
158 sf1= my34(i)*vglas(i,3) +my23(i)*vglas(i,9)
159 sf2= -mx23(i)*vglas(i,10)-mx34(i)*vglas(i,4)
160 sc5= my34(i)*vglas(i,5) +mx34(i)*vglas(i,6)
161 sc6= my23(i)*vglas(i,11)+mx23(i)*vglas(i,12)
162 c5=half*off(i)*thk(i)*c7
163 esx(i)= ss1*dhg(i,1)+ss2*dhg(i,2)
164 etmp(i,1) = c5*(esx(i) +
165 . fourth*(sc5*dhg(i,5)+sc6*dhg(i,6)))
166 emx(i)= (sf1*dhg(i,3)-sf2*dhg(i,4))*c6(i)
167 etmp(i,2) = c5*emx(i)
168 dglas(i,1) =c1m*cxx-c12*cyy
169 dglas(i,2) =c2m*cyy-c12*cxx
170 dglas(i,3) =c1m*bxx-c12*byy
171 dglas(i,4) =c2m*byy-c12*bxx
172 dglas(i,7) =c1m*cxx_k-c12*cyy_k
173 dglas(i,8) =c2m*cyy_k-c12*cxx_k
174 dglas(i,9) =c1m*bxx_k-c12*byy_k
175 dglas(i,10)=c2m*byy_k-c12*bxx_k
176 c2=fac1(i)*g(i)*shf(i)*one_over_64
177 dglas(i,5) =c2*hxx*dhg(i,5)
178 dglas(i,6) =c2*hyy*dhg(i,5)
179 dglas(i,11)=c2*hxx_k*dhg(i,6)
180 dglas(i,12)=c2*hyy_k*dhg(i,6)
181 END DO
182 END IF !(NPT==1) THEN
183 ELSE
184 DO i=jft,jlt
185 c3=four*aa(i)
186 hxx=c3*my34(i)
187 hyy=c3*mx34(i)
188 hxx_k=c3*my23(i)
189 hyy_k=c3*mx23(i)
190 cxx=hxx*dhg(i,1)
191 cyy=hyy*dhg(i,2)
192 cxx_k=hxx_k*dhg(i,1)
193 cyy_k=hyy_k*dhg(i,2)
194 bxx=hxx*dhg(i,3)
195 byy=hyy*dhg(i,4)
196 bxx_k=hxx_k*dhg(i,3)
197 byy_k=hyy_k*dhg(i,4)
198 c1m=a11(i)*fac1(i)
199 c2m=a12(i)*fac1(i)
200C -----------for energy calculation-Mean_value--------
201 c6(i)=c1(i)*fbend
202 ss1= my34(i)*vglas(i,1) +my23(i)*vglas(i,7)
203 ss2= mx23(i)*vglas(i,8) +mx34(i)*vglas(i,2)
204 sf1= my34(i)*vglas(i,3) +my23(i)*vglas(i,9)
205 sf2= -mx23(i)*vglas(i,10)-mx34(i)*vglas(i,4)
206 sc5= my34(i)*vglas(i,5) +mx34(i)*vglas(i,6)
207 sc6= my23(i)*vglas(i,11)+mx23(i)*vglas(i,12)
208 c5=half*off(i)*thk(i)*c7
209 esx(i)= ss1*dhg(i,1)+ss2*dhg(i,2)
210 etmp(i,1) = c5*(esx(i) +
211 . fourth*(sc5*dhg(i,5)+sc6*dhg(i,6)))
212 emx(i)= (sf1*dhg(i,3)-sf2*dhg(i,4))*c6(i)
213 etmp(i,2) = c5*emx(i)
214 dglas(i,1) =c1m*cxx-c2m*cyy
215 dglas(i,2) =c1m*cyy-c2m*cxx
216 dglas(i,3) =c1m*bxx-c2m*byy
217 dglas(i,4) =c1m*byy-c2m*bxx
218 dglas(i,7) =c1m*cxx_k-c2m*cyy_k
219 dglas(i,8) =c1m*cyy_k-c2m*cxx_k
220 dglas(i,9) =c1m*bxx_k-c2m*byy_k
221 dglas(i,10)=c1m*byy_k-c2m*bxx_k
222 c2=fac1(i)*g(i)*shf(i)*one_over_64
223 dglas(i,5) =c2*hxx*dhg(i,5)
224 dglas(i,6) =c2*hyy*dhg(i,5)
225 dglas(i,11)=c2*hxx_k*dhg(i,6)
226 dglas(i,12)=c2*hyy_k*dhg(i,6)
227 END DO
228 END IF !IF(MTN==58 )
229
230 IF(mtn==58 .AND.npt==1) THEN
231 DO i=jft,jlt
232C -----------for energy calculation-Mean_value--------
233 c8 = c7*off(i)
234 ss1= (my34(i)*dglas(i,1) +my23(i)*dglas(i,7))*c8
235 ss2= (mx23(i)*dglas(i,8) +mx34(i)*dglas(i,2))*c8
236 hsura=thk(i)*aa(i)
237 hvl = amu(i)*sqrt(rho(i)*area(i)*fac1(i))*off(i)
238 ssv0=my23(i)*my23(i)
239 ssv1=my34(i)*my34(i)
240 ssv2=mx23(i)*mx23(i)
241 ssv3=mx34(i)*mx34(i)
242 hxx_v= stier*(ssv1+ssv0)
243 hxy_v=-stier*(my34(i)*mx34(i)+my23(i)*mx23(i))
244 hyy_v= stier*(ssv2+ssv3)
245 aux=aa(i)*hvl
246 c1m = a11sr(i)*aux
247 c2m = a12sr(i)*aux
248 cxx_v=c1m*hxx_v
249 cyy_v=c1m*hyy_v
250 cxy_v=c2m*hxy_v
251 ss1_v=cxx_v*vhg(i,1)+cxy_v*vhg(i,2)
252 ss2_v=cyy_v*vhg(i,2)+cxy_v*vhg(i,1)
253C--------------------------------------------------
254 ss1=ss1+ss1_v
255 ss2=ss2+ss2_v
256 c2=fourth*thk(i)
257C--------NOEUDS 1,3-------
258 b13=(my13(i)*x24(i)-mx13(i)*y24(i))*hsura
259 vf(i,1,1)=vf(i,1,1)+b13*ss1
260 vf(i,1,3)=c2*ss1
261 vf(i,2,1)=vf(i,2,1)+b13*ss2
262 vf(i,2,3)=c2*ss2
263C--------NOEUDS 2,4-------
264 b24=(mx13(i)*y13(i)-my13(i)*x13(i))*hsura
265 vf(i,1,2)= vf(i,1,2)+b24*ss1
266 vf(i,1,4)=-vf(i,1,3)
267 vf(i,2,2)= vf(i,2,2)+b24*ss2
268 vf(i,2,4)=-vf(i,2,3)
269C
270 c2=z1(i)*hsura
271 vf(i,3,1)=vf(i,3,1)+c2*(ss1*y24(i)-ss2*x24(i))
272 vf(i,3,2)=vf(i,3,2)+c2*(-ss1*y13(i)+ss2*x13(i))
273C -----------for energy calculation-Mean_value--------
274 esy= ((ss1-ss1_v)*dhg(i,1)+(ss2-ss2_v)*dhg(i,2))*thk(i)
275 etmp(i,1) = esy
276 etmp(i,2) = zero
277C ----viscous energy is removed to hourglass energy----
278 tesy(i)= (ss1_v*dhg(i,1)+ss2_v*dhg(i,2))*thk(i)
279 ENDDO
280 ELSE
281 DO i=jft,jlt
282C ---------elstic increament----------
283 vglas(i,1) =vglas(i,1) +dglas(i,1)
284 vglas(i,2) =vglas(i,2) +dglas(i,2)
285 vglas(i,3) =vglas(i,3) +dglas(i,3)
286 vglas(i,4) =vglas(i,4) +dglas(i,4)
287 vglas(i,7) =vglas(i,7) +dglas(i,7)
288 vglas(i,8) =vglas(i,8) +dglas(i,8)
289 vglas(i,9) =vglas(i,9) +dglas(i,9)
290 vglas(i,10)=vglas(i,10)+dglas(i,10)
291C -----------/4 par rapport formulation-------
292 vglas(i,5) =vglas(i,5) +dglas(i,5)
293 vglas(i,6) =vglas(i,6) +dglas(i,6)
294 vglas(i,11)=vglas(i,11)+dglas(i,11)
295 vglas(i,12)=vglas(i,12)+dglas(i,12)
296C-------------enint------------------
297 eint(i,1) = eint(i,1)+etmp(i,1)
298 eint(i,2) = eint(i,2)+etmp(i,2)
299 END DO
300C
301 DO i=jft,jlt
302 IF (sigy(i)<zep9ep30) THEN
303 ufac=abs(min(fac(i,1),fac(i,2))-one)
304 sigy2 = sigy(i)*sigy(i)
305 svm =zero
306 sxy0=zero
307 IF (ufac<tol) THEN
308 sxy0= vstre(i,1)*vstre(i,1)+vstre(i,2)*vstre(i,2)
309 . -vstre(i,1)*vstre(i,2)+three*vstre(i,3)*vstre(i,3)
310 mxy0= mstre(i,1)*mstre(i,1)+mstre(i,2)*mstre(i,2)
311 . -mstre(i,1)*mstre(i,2)+three*mstre(i,3)*mstre(i,3)
312 cnn=coef
313 cmm=coef*thk(i)*one_over_16
314 cnnx=cnn*vglas(i,1)
315 cnny=cnn*vglas(i,2)
316 cnnx_k=cnn*vglas(i,7)
317 cnny_k=cnn*vglas(i,8)
318 cmmx=cmm*vglas(i,3)
319 cmmy=cmm*vglas(i,4)
320 cmmx_k=cmm*vglas(i,9)
321 cmmy_k=cmm*vglas(i,10)
322 sxy0= sxy0+cnnx*cnnx+cnny*cnny-cnnx*cnny
323 mxy0= mxy0+cmmx*cmmx+cmmy*cmmy-cmmx*cmmy
324 sxy0= sxy0+cnnx_k*cnnx_k+cnny_k*cnny_k-cnnx_k*cnny_k
325 mxy0= mxy0+cmmx_k*cmmx_k+cmmy_k*cmmy_k-cmmx_k*cmmy_k
326 sxy0= sxy0+abs(cnnx*(two*cnnx_k-cnny_k)
327 1 +cnny*(two*cnny_k-cnnx_k))
328 mxy0= mxy0+abs(cmmx*(two*cmmx_k-cmmy_k)
329 1 +cmmy*(two*cmmy_k-cmmx_k))
330 svm = sxy0+coef1*mxy0
331 ENDIF
332C-----------ELASTIC-PLASTIC global yield criterion------------
333 IF (ufac>=tol.OR.svm>sigy2) THEN
334 eh1=min(sxy0/max(sigy2,tol),one)
335 eh1=max(coefh*eh1,(one-fac(i,1)))
336 eh2=max(coefh,(one-fac(i,2)))
337 IF (esx(i)<zero) eh1=zero
338 IF (emx(i)<zero) eh2=zero
339C-------------MEMBRANE------------------------
340 vglas(i,1)=vglas(i,1)-eh1*dglas(i,1)
341 vglas(i,2)=vglas(i,2)-eh1*dglas(i,2)
342 vglas(i,7)=vglas(i,7)-eh1*dglas(i,7)
343 vglas(i,8)=vglas(i,8)-eh1*dglas(i,8)
344C-------------FLEXION------------------------
345 vglas(i,3) =vglas(i,3) -eh2*dglas(i,3)
346 vglas(i,4) =vglas(i,4) -eh2*dglas(i,4)
347 vglas(i,9) =vglas(i,9) -eh2*dglas(i,9)
348 vglas(i,10)=vglas(i,10)-eh2*dglas(i,10)
349 ENDIF
350 ENDIF
351 END DO
352
353 DO i=jft,jlt
354C -----------for energy calculation-Mean_value--------
355 c8 = c7*off(i)
356 ss1= (my34(i)*vglas(i,1) +my23(i)*vglas(i,7))*c8
357 ss2= (mx23(i)*vglas(i,8) +mx34(i)*vglas(i,2))*c8
358 sf1= (my34(i)*vglas(i,3) +my23(i)*vglas(i,9))*c8
359 sf2= -(mx23(i)*vglas(i,10)+mx34(i)*vglas(i,4))*c8
360 hsura=thk(i)*aa(i)
361 c2 = c8*thk(i)
362 sc5=(my34(i)*vglas(i,5) +mx34(i)*vglas(i,6) )*c2
363 sc6=(my23(i)*vglas(i,11)+mx23(i)*vglas(i,12))*c2
364 ss3=sc5+sc6
365 hvl = amu(i)*sqrt(rho(i)*area(i)*fac1(i))*off(i)
366 ssv0=my23(i)*my23(i)
367 ssv1=my34(i)*my34(i)
368 ssv2=mx23(i)*mx23(i)
369 ssv3=mx34(i)*mx34(i)
370 hxx_v= stier*(ssv1+ssv0)
371 hxy_v=-stier*(my34(i)*mx34(i)+my23(i)*mx23(i))
372 hyy_v= stier*(ssv2+ssv3)
373 c2 =hvl*gsr(i)*shfsr(i)*undouzsr
374 cxz_v=(ssv1+ssv3)*c2
375 cyz_v=(ssv2+ssv0)*c2
376 aux=aa(i)*hvl
377 c1m = a11sr(i)*aux
378 c2m = a12sr(i)*aux
379 cxx_v=c1m*hxx_v
380 cyy_v=c1m*hyy_v
381 cxy_v=c2m*hxy_v
382 ss1_v=cxx_v*vhg(i,1)+cxy_v*vhg(i,2)
383 ss2_v=cyy_v*vhg(i,2)+cxy_v*vhg(i,1)
384 sf1_v= (cxx_v*vhg(i,3)+cxy_v*vhg(i,4))*fbend_v
385 sf2_v=(-cyy_v*vhg(i,4)-cxy_v*vhg(i,3))*fbend_v
386 sc5_v=cxz_v*vhg(i,5)*hsura
387 sc6_v=cyz_v*vhg(i,6)*hsura
388 ss3_v = sc5_v+sc6_v
389C--------------------------------------------------
390 ss1=ss1+ss1_v
391 ss2=ss2+ss2_v
392 ss3=ss3+ss3_v
393 sc5=sc5+sc5_v
394 sc6=sc6+sc6_v
395 sf1=sf1+sf1_v
396 sf2=sf2+sf2_v
397 y13s= my13(i)*ss3
398 x13s= mx13(i)*ss3
399 y34s6=my34(i)*sc6
400 y23s5=my23(i)*sc5
401 x23s5=mx23(i)*sc5
402 x34s6=mx34(i)*sc6
403 c2=fourth*thk(i)
404C--------NOEUDS 1,3-------
405 b13=(my13(i)*x24(i)-mx13(i)*y24(i))*hsura
406 vf(i,1,1)=vf(i,1,1)+b13*ss1
407 vf(i,1,3)=c2*ss1
408 vf(i,2,1)=vf(i,2,1)+b13*ss2
409 vf(i,2,3)=c2*ss2
410 vf(i,3,3)=ss3
411C--------NOEUDS 2,4-------
412 b24=(mx13(i)*y13(i)-my13(i)*x13(i))*hsura
413 vf(i,1,2)= vf(i,1,2)+b24*ss1
414 vf(i,1,4)=-vf(i,1,3)
415 vf(i,2,2)= vf(i,2,2)+b24*ss2
416 vf(i,2,4)=-vf(i,2,3)
417 vf(i,3,4)=-vf(i,3,3)
418C--------NOEUDS 1,3-------
419 c3=c6(i)*b13
420 c4=c6(i)*c2
421 vm(i,1,1)=vm(i,1,1)+c3*sf2+y23s5+y34s6
422 vm(i,1,3)=vm(i,1,3)+c4*sf2-y13s
423 vm(i,2,1)=vm(i,2,1)+c3*sf1-x23s5-x34s6
424 vm(i,2,3)=vm(i,2,3)+c4*sf1+x13s
425C--------NOEUDS 2,4-------
426 c3=c6(i)*b24
427 vm(i,1,2)=vm(i,1,2)+c3*sf2+y23s5-y34s6
428 vm(i,1,4)=vm(i,1,4)-c4*sf2-y13s
429 vm(i,2,2)=vm(i,2,2)+c3*sf1-x23s5+x34s6
430 vm(i,2,4)=vm(i,2,4)-c4*sf1+x13s
431C
432 c2=z1(i)*hsura
433 vf(i,3,1)=vf(i,3,1)+c2*(ss1*y24(i)-ss2*x24(i))
434 vf(i,3,2)=vf(i,3,2)+c2*(-ss1*y13(i)+ss2*x13(i))
435C
436C -----------for energy calculation-Mean_value--------
437 esy= ((ss1-ss1_v)*dhg(i,1)+(ss2-ss2_v)*dhg(i,2))*thk(i)+
438 . fourth*((sc5-sc5_v)*dhg(i,5)+(sc6-sc6_v)*dhg(i,6))
439 etmp(i,1) = half*esy
440 emy= (sf1-sf1_v)*dhg(i,3)-(sf2-sf2_v)*dhg(i,4)
441 etmp(i,2) = half*c6(i)*emy*thk(i)
442C ----viscous energy is removed to hourglass energy----
443 tesy(i)= (ss1_v*dhg(i,1)+ss2_v*dhg(i,2))*thk(i)+
444 . (sf1_v*dhg(i,3)-sf2_v*dhg(i,4))*thk(i)*c6(i)+
445 . fourth*(sc5_v*dhg(i,5)+sc6_v*dhg(i,6))
446 ENDDO
447 END if!(MTN==58 .AND.NPT==1) THEN
448C
449 DO i=jft,jlt
450 eint(i,1) = eint(i,1)+etmp(i,1)
451 eint(i,2) = eint(i,2)+etmp(i,2)
452 END DO
453
454 ic=1
455 jst(ic)=jft
456 DO j=jft+1,jlt
457 IF (ipartc(j)/=ipartc(j-1)) THEN
458 ic=ic+1
459 jst(ic)=j
460 ENDIF
461 ENDDO
462 jst(ic+1)=jlt+1
463 IF(ic==1) THEN
464 mx = ipartc(jft)
465 DO i=jft,jlt
466 evis(8,mx)=evis(8,mx) + tesy(i)
467 ENDDO
468 ELSEIF(ic==2.AND.kfts>0)THEN
469 mx = ipartc(jft)
470 DO i=jft, kfts-1
471 evis(8,mx)=evis(8,mx) + tesy(i)
472 ENDDO
473 mx = ipartc(jlt)
474 DO i=kfts,jlt
475 evis(8,mx)=evis(8,mx) + tesy(i)
476 ENDDO
477
478 ELSE
479 DO ii=1,ic
480 mx=ipartc(jst(ii))
481 IF (jst(ii+1)-jst(ii)>15) THEN
482 DO j=jst(ii),jst(ii+1)-1
483 evis(8,mx)=evis(8,mx) + tesy(j)
484 ENDDO
485 ELSE
486 DO j=jst(ii),jst(ii+1)-1
487 evis(8,mx)=evis(8,mx) + tesy(j)
488 ENDDO
489 ENDIF
490 ENDDO
491 ENDIF
492 RETURN
493 END
494
495!||====================================================================
496!|| czfintnm ../engine/source/elements/shell/coquez/czfintn.F
497!||--- called by ------------------------------------------------------
498!|| czfintnm1 ../engine/source/elements/shell/coquez/czfintn.F
499!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
500!|| czforc3_crk ../engine/source/elements/xfem/czforc3_crk.F
501!||====================================================================
502 SUBROUTINE czfintnm(JFT ,JLT ,THK ,AA ,VHG ,
503 2 X13 ,X24 ,Y13 ,Y24 ,VF ,
504 3 G ,RHO ,AREA ,AMU ,DT1 ,
505 4 OFF ,IPARTC,EVIS,KFTS )
506C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
507C CALCUL viscous VF(3,*) due au shear when npt=1
508C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
509#include "implicit_f.inc"
510#include "param_c.inc"
511#include "mvsiz_p.inc"
512C-----------------------------------------------
513C D U M M Y A R G U M E N T S
514C-----------------------------------------------
515 INTEGER JFT,JLT,IPARTC(*)
516 my_real
517 . THK(*) ,AA(*) ,VHG(MVSIZ,6),
518 . X13(*) ,X24(*) ,Y13(*) ,Y24(*) ,
519 . VF(MVSIZ,3,4),G(*),
520 . RHO(*) ,AREA(*), OFF(*),AMU(*) ,
521 . EVIS(NPSAV,*),DT1
522 INTEGER KFTS
523C-----------------------------------------------
524C L O C A L V A R I A B L E S
525C-----------------------------------------------
526 INTEGER I,MX, IC, II, J, JST(MVSIZ+1)
527 my_real
528 . C2,HSURA,MX23,MY23,MX34,MY34,
529 . SC6_V,SC5_V,CXZ_V,CYZ_V,SS3_V,HVL,
530 . TESY(MVSIZ)
531C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
532 DO i=jft,jlt
533 mx23=half*(x13(i)-x24(i))
534 mx34=half*(x13(i)+x24(i))
535 my23=half*(y13(i)-y24(i))
536 my34=half*(y13(i)+y24(i))
537 c2 =one_over_12*g(i)*rho(i)*area(i)
538 hvl = twenty5*amu(i)*sqrt(c2)*off(i)
539 cxz_v=(my34*my34+mx34*mx34)*hvl
540 cyz_v=(my23*my23+mx23*mx23)*hvl
541 hsura=thk(i)*aa(i)
542 sc5_v=cxz_v*vhg(i,5)*hsura
543 sc6_v=cyz_v*vhg(i,6)*hsura
544 ss3_v = sc5_v+sc6_v
545C--------NOEUDS 1,3-------
546 vf(i,3,3)=vf(i,3,3)+ss3_v
547C--------NOEUDS 2,4-------
548 vf(i,3,4)=vf(i,3,4)-ss3_v
549 tesy(i)= (sc5_v*vhg(i,5)+sc6_v*vhg(i,6))*dt1
550 ENDDO
551C
552 mx = ipartc(jft)
553 DO i=jft,jlt
554 evis(8,mx)=evis(8,mx) + tesy(i)
555 ENDDO
556C
557 RETURN
558 END
559!||====================================================================
560!|| czfintnrz ../engine/source/elements/shell/coquez/czfintn.F
561!||--- called by ------------------------------------------------------
562!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
563!||====================================================================
564 SUBROUTINE czfintnrz(JFT ,JLT ,THK ,C1 ,AA ,VHG ,
565 2 X13 ,X24 ,Y13 ,Y24 ,Z1 ,MX23 ,
566 3 MX13 ,MX34 ,MY13 ,MY23 ,MY34 ,VGLAS,
567 4 VSTRE,MSTRE,VF ,VM ,FAC ,A11 ,
568 5 A12 ,G ,SHF ,SIGY ,OFF ,FAC1 ,
569 6 RHO ,AREA , DT1,EINT ,AMU ,VHI ,
570 7 NPT ,IPARTC,EVIS,KFTS ,GSR ,
571 8 A11SR,A12SR ,NUSR,SHFSR,BMKRZ,BMERZ,
572 9 VHGZK,VHGZE ,KRZ ,VMZ ,NEL )
573C elasto-plastic hourglass stresses+linear damping
574C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
575C CALCUL VF(3,4),VM(2,4) : FORCES INTERNES GENERALISEES
576C ENTREES : THK,1/AREA,PME,X13,X24,Y13,Y24,MX13,MX34,MY13,MY34
577C DHG(6) : INCREMENT DE LA DEF HOURGLASS GENERALISEE
578C VGLAS(12): DEF HOURGLASS GENERALISEE
579C VGLAS(1-6)->ETA; VGLAS(7-12)->KSI;
580C -shear: VGLAS(13-14)->ETA; VGLAS(15-16)->KSI;
581C -Sig_rz: VGLAS(17)->ETA; VGLAS(18)->KSI; VGLAS(19)->const
582C FAC=Et/E
583C SORTIES : VF,VM
584C VF(3,4),VM(2,4) :
585C J=1,2 STOCK LA PARTIE ANTISYM (F(1)=-F(3),F(2)=-F(4))
586C J=3,4 STOCK LA PARTIE SYM (F(1)=F(3),F(2)=F(4))
587C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
588#include "implicit_f.inc"
589#include "param_c.inc"
590#include "mvsiz_p.inc"
591C-----------------------------------------------
592C D U M M Y A R G U M E N T S
593C-----------------------------------------------
594 INTEGER JFT,JLT,NPT,IPARTC(*),NEL
595 my_real
596 . THK(*) ,C1(*) ,AA(*) ,VHG(MVSIZ,6),
597 . X13(*) ,X24(*) ,Y13(*) ,Y24(*) ,Z1(*),
598 . MX13(*),MX23(*),MX34(*),MY13(*),MY23(*),MY34(*),
599 . VGLAS(NEL,19),VSTRE(NEL,5),MSTRE(NEL,3),VF(MVSIZ,3,4),VM(MVSIZ,2,4),
600 . FAC(MVSIZ,2),A11(*) ,A12(*) ,G(*),SHF(*),VHI(MVSIZ,3),SIGY(*),
601 . FAC1(*) ,RHO(*) ,AREA(*), DT1,OFF(*),EINT(JLT,2),AMU(*) ,
602 . EVIS(NPSAV,*),VMZ(MVSIZ,4),
603 . VHGZK(MVSIZ,5),VHGZE(MVSIZ,5),KRZ(*),
604 . BM0RZ(MVSIZ,4,4),BMKRZ(MVSIZ,4,4),BMERZ(MVSIZ,4,4),
605 . GSR(*), A11SR(*), A12SR(*), NUSR(*), SHFSR(*)
606 INTEGER KFTS
607C-----------------------------------------------
608C L O C A L V A R I A B L E S
609C-----------------------------------------------
610 INTEGER I,MX, IC, II, J, JST(MVSIZ+1)
611 my_real
612 . stier,c2,c3,c4,c5,c7,c8,hsura,sxx,syy,svm,svmn,
613 . ss1,ss2,ss3,sc6,sc5,cxx_k,cyy_k,hxx,hyy,sf1,sf2,b13,b24,
614 . cxx,cyy,y13s,x13s,y34s6,y23s5,x23s5,x34s6,hxx_k,hyy_k,
615 . bxx,byy,bxx_k,byy_k,
616 . c1m,c2m,c1b,c2b,cmm,cnn,ufac,tol,coef,eh1,eh2,
617 . cnnx,cnny,cnnx_k,cnny_k,cmmx,cmmy,cmmx_k,cmmy_k,sxy0,mxy0,
618 . svmc,sigy2,svm0,svmc_k,a_bk,a_b,a_k,undouzsr,aux,
619 . ssv0,ssv1,ssv2,ssv3,sc6_v,sc5_v,cxx_v,cyy_v,cxy_v,cxz_v,cyz_v,
620 . hxx_v,hyy_v,hxy_v,ss1_v,ss2_v,ss3_v,sf1_v,sf2_v,hvl,crz,
621 . emy,ecx,ecy,esy,coef1,coefh,fbend,fbend_v,
622 . cnnxy,cnnxy_k,cmmxy,cmmxy_k,ssz1,ssz2,c3m,erz,
623 . c6(mvsiz), esx(mvsiz), tesy(mvsiz), emx(mvsiz),
624 . dglas(mvsiz,18),
625 . dhg(mvsiz,6), dhgzk(mvsiz,5),dhgze(mvsiz,5), etmp(mvsiz,2)
626C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
627 c7=four_over_3
628 coefh= zep999
629 coef= zep85
630 stier = fivep333
631 tol= em18
632 fbend_v =threep464
633 IF (npt==0) THEN
634 coef1=sixteen
635 ELSE
636 coef1=twenty5
637 ENDIF
638 IF (npt==1) THEN
639 fbend =zero
640 fbend_v =zero
641 ELSE
642 fbend =one_over_12
643 ENDIF
644 undouzsr=sqrt(one_over_12)
645C
646 DO i=jft,jlt
647 dhg(i,1)=vhg(i,1)*dt1
648 dhg(i,2)=vhg(i,2)*dt1
649 dhg(i,3)=vhg(i,3)*dt1
650 dhg(i,4)=vhg(i,4)*dt1
651 dhg(i,5)=vhg(i,5)*dt1
652 dhg(i,6)=vhg(i,6)*dt1
653 dhgzk(i,1)=vhgzk(i,1)*dt1
654 dhgzk(i,2)=vhgzk(i,2)*dt1
655 dhgzk(i,3)=vhgzk(i,3)*dt1
656 dhgzk(i,4)=vhgzk(i,4)*dt1
657 dhgzk(i,5)=vhgzk(i,5)*dt1
658 dhgze(i,1)=vhgze(i,1)*dt1
659 dhgze(i,2)=vhgze(i,2)*dt1
660 dhgze(i,3)=vhgze(i,3)*dt1
661 dhgze(i,4)=vhgze(i,4)*dt1
662 dhgze(i,5)=vhgze(i,5)*dt1
663 c3=four*aa(i)
664 hxx=c3*my34(i)
665 hyy=c3*mx34(i)
666 hxx_k=c3*my23(i)
667 hyy_k=c3*mx23(i)
668 cxx=hxx*dhg(i,1)
669 cyy=hyy*dhg(i,2)
670 cxx_k=hxx_k*dhg(i,1)
671 cyy_k=hyy_k*dhg(i,2)
672 bxx=hxx*dhg(i,3)
673 byy=hyy*dhg(i,4)
674 bxx_k=hxx_k*dhg(i,3)
675 byy_k=hyy_k*dhg(i,4)
676 c1m=a11(i)*fac1(i)
677 c2m=a12(i)*fac1(i)
678 c3m=g(i)*fac1(i)
679 crz=half*krz(i)*fac1(i)
680C -----------for energy calculation-Mean_value--------
681 c6(i)=c1(i)*fbend
682 ss1= my34(i)*vglas(i,1) +my23(i)*vglas(i,7)
683 ss2= mx23(i)*vglas(i,8) +mx34(i)*vglas(i,2)
684 ss1 = ss1 - mx34(i)*vglas(i,13)+mx23(i)*vglas(i,15)
685 . +half*(mx34(i)*vglas(i,17)-mx23(i)*vglas(i,18))
686 ss2 = ss2 + my34(i)*vglas(i,13) -my23(i)*vglas(i,15)
687 . +half*(my34(i)*vglas(i,17)-my23(i)*vglas(i,18))
688C-------add Srz---------
689 c8 = c7*off(i)
690 c3= half*c8
691 ss1= ss1 + c3*(mx34(i)*vglas(i,17)-mx23(i)*vglas(i,18))
692 ss2= ss2 + c3*(my34(i)*vglas(i,17)-my23(i)*vglas(i,18))
693 sf1= my34(i)*vglas(i,3) +my23(i)*vglas(i,9)
694 sf2= -mx23(i)*vglas(i,10)-mx34(i)*vglas(i,4)
695 sf1=sf1+(-mx34(i)*vglas(i,14)+mx23(i)*vglas(i,16))
696 sf2=sf2-(my34(i)*vglas(i,14)-my23(i)*vglas(i,16))
697 sc5= my34(i)*vglas(i,5) +mx34(i)*vglas(i,6)
698 sc6= my23(i)*vglas(i,11)+mx23(i)*vglas(i,12)
699 erz = vglas(i,1)*dhgze(i,1)-vglas(i,2)*dhgze(i,2)-
700 . vglas(i,7)*dhgzk(i,1)+vglas(i,8)*dhgzk(i,2)+
701 . vglas(i,13)*dhgze(i,3)+vglas(i,15)*dhgzk(i,3)+
702 . half*(vglas(i,17)*dhgze(i,4)+vglas(i,18)*dhgzk(i,4))
703 c5=half*off(i)*thk(i)*c7
704 esx(i)= ss1*dhg(i,1)+ss2*dhg(i,2)+erz*area(i)*fourth
705 etmp(i,1) = c5*(esx(i) +
706 . fourth*(sc5*dhg(i,5)+sc6*dhg(i,6)))
707 emx(i)= (sf1*dhg(i,3)-sf2*dhg(i,4))*c6(i)
708 etmp(i,2) = c5*emx(i)
709C
710 dglas(i,1) =c1m*cxx-c2m*cyy+c1m*dhgze(i,1)+c2m*dhgze(i,2)
711 dglas(i,2) =c1m*cyy-c2m*cxx-c1m*dhgze(i,2)-c2m*dhgze(i,1)
712 dglas(i,3) =c1m*bxx-c2m*byy
713 dglas(i,4) =c1m*byy-c2m*bxx
714 dglas(i,7) =c1m*cxx_k-c2m*cyy_k-c1m*dhgzk(i,1)-c2m*dhgzk(i,2)
715 dglas(i,8) =c1m*cyy_k-c2m*cxx_k+c1m*dhgzk(i,2)+c2m*dhgzk(i,1)
716 dglas(i,9) =c1m*bxx_k-c2m*byy_k
717 dglas(i,10)=c1m*byy_k-c2m*bxx_k
718 c2=fac1(i)*g(i)*shf(i)*one_over_64
719 dglas(i,5) =c2*hxx*dhg(i,5)
720 dglas(i,6) =c2*hyy*dhg(i,5)
721 dglas(i,11)=c2*hxx_k*dhg(i,6)
722 dglas(i,12)=c2*hyy_k*dhg(i,6)
723C --------add shear due to rz dof--no sign change like 2,4,7,9-------
724 dglas(i,13)=c3m*(hxx*dhg(i,2)-hyy*dhg(i,1)+dhgze(i,3))
725 dglas(i,14)=c3m*(hxx*dhg(i,4)-hyy*dhg(i,3))
726 dglas(i,15)=c3m*(-hxx_k*dhg(i,2)+hyy_k*dhg(i,1)+dhgzk(i,3))
727 dglas(i,16)=c3m*(-hxx_k*dhg(i,4)+hyy_k*dhg(i,3))
728C --------add srz dof--------
729 dglas(i,17)=crz*dhgze(i,5)
730 dglas(i,18)=crz*dhgzk(i,5)
731 END DO
732C
733 DO i=jft,jlt
734C ---------elstic increament----------
735 vglas(i,1) =vglas(i,1) +dglas(i,1)
736 vglas(i,2) =vglas(i,2) +dglas(i,2)
737 vglas(i,3) =vglas(i,3) +dglas(i,3)
738 vglas(i,4) =vglas(i,4) +dglas(i,4)
739 vglas(i,5) =vglas(i,5) +dglas(i,5)
740 vglas(i,6) =vglas(i,6) +dglas(i,6)
741 vglas(i,7) =vglas(i,7) +dglas(i,7)
742 vglas(i,8) =vglas(i,8) +dglas(i,8)
743 vglas(i,9) =vglas(i,9) +dglas(i,9)
744 vglas(i,10)=vglas(i,10)+dglas(i,10)
745 vglas(i,11)=vglas(i,11)+dglas(i,11)
746 vglas(i,12)=vglas(i,12)+dglas(i,12)
747 vglas(i,13)=vglas(i,13)+dglas(i,13)
748 vglas(i,14)=vglas(i,14)+dglas(i,14)
749 vglas(i,15)=vglas(i,15)+dglas(i,15)
750 vglas(i,16)=vglas(i,16)+dglas(i,16)
751 vglas(i,17)=vglas(i,17)+dglas(i,17)
752 vglas(i,18)=vglas(i,18)+dglas(i,18)
753C ---------enint----------
754 eint(i,1) = eint(i,1)+etmp(i,1)
755 eint(i,2) = eint(i,2)+etmp(i,2)
756 END DO
757C
758 DO i=jft,jlt
759 IF (sigy(i)<zep9ep30) THEN
760 ufac=abs(min(fac(i,1),fac(i,2))-one)
761 sigy2 = sigy(i)*sigy(i)
762 svm =zero
763 sxy0=zero
764 IF (ufac<tol) THEN
765 sxy0= vstre(i,1)*vstre(i,1)+vstre(i,2)*vstre(i,2)
766 . -vstre(i,1)*vstre(i,2)+three*vstre(i,3)*vstre(i,3)
767 mxy0= mstre(i,1)*mstre(i,1)+mstre(i,2)*mstre(i,2)
768 . -mstre(i,1)*mstre(i,2)+three*mstre(i,3)*mstre(i,3)
769 cnn=coef
770 cmm=coef*thk(i)*one_over_16
771 cnnx=cnn*vglas(i,1)
772 cnny=cnn*vglas(i,2)
773 cnnxy=cnn*vglas(i,13)
774 cnnx_k=cnn*vglas(i,7)
775 cnny_k=cnn*vglas(i,8)
776 cnnxy_k=cnn*vglas(i,15)
777 cmmx=cmm*vglas(i,3)
778 cmmy=cmm*vglas(i,4)
779 cmmxy=cmm*vglas(i,14)
780 cmmx_k=cmm*vglas(i,9)
781 cmmy_k=cmm*vglas(i,10)
782 cmmxy_k=cmm*vglas(i,16)
783C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
784 sxy0= sxy0+cnnx*cnnx+cnny*cnny-cnnx*cnny+three*cnnxy*cnnxy
785 mxy0= mxy0+cmmx*cmmx+cmmy*cmmy-cmmx*cmmy+three*cmmxy*cmmxy
786 sxy0= sxy0+cnnx_k*cnnx_k+cnny_k*cnny_k-cnnx_k*cnny_k+
787 1 three*cnnxy_k*cnnxy_k
788 mxy0= mxy0+cmmx_k*cmmx_k+cmmy_k*cmmy_k-cmmx_k*cmmy_k+
789 1 three*cmmxy_k*cmmxy_k
790 sxy0= sxy0+abs(cnnx*(two*cnnx_k-cnny_k)
791 1 +cnny*(two*cnny_k-cnnx_k))
792 mxy0= mxy0+abs(cmmx*(two*cmmx_k-cmmy_k)
793 1 +cmmy*(two*cmmy_k-cmmx_k))
794 svm = sxy0+coef1*mxy0
795 ENDIF
796 IF (ufac>=tol.OR.svm>sigy2) THEN
797 eh1=min(sxy0/max(sigy2,tol),one)
798 eh1=max(coefh*eh1,(one-fac(i,1)))
799 eh2=max(coefh,(one-fac(i,2)))
800 IF (esx(i)<zero) eh1=zero
801 IF (emx(i)<zero) eh2=zero
802C-------------MEMBRANE------------------------
803 vglas(i,1)=vglas(i,1)-eh1*dglas(i,1)
804 vglas(i,2)=vglas(i,2)-eh1*dglas(i,2)
805 vglas(i,7)=vglas(i,7)-eh1*dglas(i,7)
806 vglas(i,8)=vglas(i,8)-eh1*dglas(i,8)
807 vglas(i,13)=vglas(i,13)-eh1*dglas(i,13)
808 vglas(i,15)=vglas(i,15)-eh1*dglas(i,15)
809C-------------FLEXION------------------------
810 vglas(i,3) =vglas(i,3) -eh2*dglas(i,3)
811 vglas(i,4) =vglas(i,4) -eh2*dglas(i,4)
812 vglas(i,9) =vglas(i,9) -eh2*dglas(i,9)
813 vglas(i,10)=vglas(i,10)-eh2*dglas(i,10)
814 vglas(i,14)=vglas(i,14)-eh1*dglas(i,14)
815 vglas(i,16)=vglas(i,16)-eh1*dglas(i,16)
816 ENDIF
817 ENDIF
818 END DO
819 DO i=jft,jlt
820C -----------for energy calculation-Mean_value--------
821 c8 = c7*off(i)
822 ss1= (my34(i)*vglas(i,1) +my23(i)*vglas(i,7))*c8
823 ss2= (mx23(i)*vglas(i,8) +mx34(i)*vglas(i,2))*c8
824C-------add shear
825 ss1 = ss1 + c8*(-mx34(i)*vglas(i,13)+mx23(i)*vglas(i,15))
826 ss2 = ss2 + c8*( my34(i)*vglas(i,13) -my23(i)*vglas(i,15))
827C-------add Srz---------
828 c3= half*c8
829 ss1= ss1 + c3*(mx34(i)*vglas(i,17)-mx23(i)*vglas(i,18))
830 ss2= ss2 + c3*(my34(i)*vglas(i,17)-my23(i)*vglas(i,18))
831C-------save for Fz--and add coupling terms(warped)------
832 ssz1 = ss1
833 ssz2 = ss2
834 sf1= (my34(i)*vglas(i,3) +my23(i)*vglas(i,9))*c8
835 sf2= -(mx23(i)*vglas(i,10)+mx34(i)*vglas(i,4))*c8
836C-------add shear SF1->y;SF2->x
837 sf1=sf1+(-mx34(i)*vglas(i,14)+mx23(i)*vglas(i,16))*c8
838 sf2=sf2-(my34(i)*vglas(i,14)-my23(i)*vglas(i,16))*c8
839 hsura=thk(i)*aa(i)
840 c2 = c8*thk(i)
841 sc5=(my34(i)*vglas(i,5) +mx34(i)*vglas(i,6) )*c2
842 sc6=(my23(i)*vglas(i,11)+mx23(i)*vglas(i,12))*c2
843 ss3=sc5+sc6
844 hvl = amu(i)*sqrt(rho(i)*area(i)*fac1(i))*off(i)
845 ssv0=my23(i)*my23(i)
846 ssv1=my34(i)*my34(i)
847 ssv2=mx23(i)*mx23(i)
848 ssv3=mx34(i)*mx34(i)
849 hxx_v= stier*(ssv1+ssv0)
850 hxy_v=-stier*(my34(i)*mx34(i)+my23(i)*mx23(i))
851 hyy_v= stier*(ssv2+ssv3)
852 c2 =hvl*gsr(i)*shfsr(i)*undouzsr
853 cxz_v=(ssv1+ssv3)*c2
854 cyz_v=(ssv2+ssv0)*c2
855 aux=aa(i)*hvl
856 c1m = a11sr(i)*aux
857 c2m = a12sr(i)*aux
858 cxx_v=c1m*hxx_v
859 cyy_v=c1m*hyy_v
860 cxy_v=c2m*hxy_v
861 ss1_v=cxx_v*vhg(i,1)+cxy_v*vhg(i,2)
862 ss2_v=cyy_v*vhg(i,2)+cxy_v*vhg(i,1)
863 sf1_v= (cxx_v*vhg(i,3)+cxy_v*vhg(i,4))*fbend_v
864 sf2_v=(-cyy_v*vhg(i,4)-cxy_v*vhg(i,3))*fbend_v
865 sc5_v=cxz_v*vhg(i,5)*hsura
866 sc6_v=cyz_v*vhg(i,6)*hsura
867 ss3_v = sc5_v+sc6_v
868C--------------------------------------------------
869 ss1=ss1+ss1_v
870 ss2=ss2+ss2_v
871 ss3=ss3+ss3_v
872 sc5=sc5+sc5_v
873 sc6=sc6+sc6_v
874 sf1=sf1+sf1_v
875 sf2=sf2+sf2_v
876 y13s= my13(i)*ss3
877 x13s= mx13(i)*ss3
878 y34s6=my34(i)*sc6
879 y23s5=my23(i)*sc5
880 x23s5=mx23(i)*sc5
881 x34s6=mx34(i)*sc6
882 c2=fourth*thk(i)
883C--------NOEUDS 1,3-------
884 b13=(my13(i)*x24(i)-mx13(i)*y24(i))*hsura
885 vf(i,1,1)=vf(i,1,1)+b13*ss1
886 vf(i,1,3)=c2*ss1
887 vf(i,2,1)=vf(i,2,1)+b13*ss2
888 vf(i,2,3)=c2*ss2
889 vf(i,3,3)=ss3
890C--------NOEUDS 2,4-------
891 b24=(mx13(i)*y13(i)-my13(i)*x13(i))*hsura
892 vf(i,1,2)= vf(i,1,2)+b24*ss1
893 vf(i,1,4)=-vf(i,1,3)
894 vf(i,2,2)= vf(i,2,2)+b24*ss2
895 vf(i,2,4)=-vf(i,2,3)
896 vf(i,3,4)=-vf(i,3,3)
897C--------NOEUDS 1,3-------
898 c3=c6(i)*b13
899 c4=c6(i)*c2
900 vm(i,1,1)=vm(i,1,1)+c3*sf2+y23s5+y34s6
901 vm(i,1,3)=vm(i,1,3)+c4*sf2-y13s
902 vm(i,2,1)=vm(i,2,1)+c3*sf1-x23s5-x34s6
903 vm(i,2,3)=vm(i,2,3)+c4*sf1+x13s
904C--------NOEUDS 2,4-------
905 c3=c6(i)*b24
906 vm(i,1,2)=vm(i,1,2)+c3*sf2+y23s5-y34s6
907 vm(i,1,4)=vm(i,1,4)-c4*sf2-y13s
908 vm(i,2,2)=vm(i,2,2)+c3*sf1-x23s5+x34s6
909 vm(i,2,4)=vm(i,2,4)-c4*sf1+x13s
910C
911 c2=z1(i)*hsura
912 vf(i,3,1)=vf(i,3,1)+c2*(ssz1*y24(i)-ssz2*x24(i))
913 vf(i,3,2)=vf(i,3,2)+c2*(-ssz1*y13(i)+ssz2*x13(i))
914C-----------non const mz----17,18--
915 c2=thk(i)*one_over_6*off(i)
916 c3=thk(i)*third*off(i)
917C
918 j=1
919 vmz(i,j)= vmz(i,j)
920 . +c2*(vglas(i,17)*bmerz(i,4,j)+vglas(i,18)*bmkrz(i,4,j))
921 . +c3*(bmerz(i,1,j)*vglas(i,1)-bmkrz(i,1,j)*vglas(i,7)
922 . -bmerz(i,2,j)*vglas(i,2)+bmkrz(i,2,j)*vglas(i,8)
923 . +bmerz(i,3,j)*vglas(i,13)+bmkrz(i,3,j)*vglas(i,15))
924 j=2
925 vmz(i,j)= vmz(i,j)
926 . +c2*(vglas(i,17)*bmerz(i,4,j)+vglas(i,18)*bmkrz(i,4,j))
927 . +c3*(bmerz(i,1,j)*vglas(i,1)-bmkrz(i,1,j)*vglas(i,7)
928 . -bmerz(i,2,j)*vglas(i,2)+bmkrz(i,2,j)*vglas(i,8)
929 . +bmerz(i,3,j)*vglas(i,13)+bmkrz(i,3,j)*vglas(i,15))
930 j=3
931 vmz(i,j)= vmz(i,j)
932 . +c2*(vglas(i,17)*bmerz(i,4,j)+vglas(i,18)*bmkrz(i,4,j))
933 . +c3*(bmerz(i,1,j)*vglas(i,1)-bmkrz(i,1,j)*vglas(i,7)
934 . -bmerz(i,2,j)*vglas(i,2)+bmkrz(i,2,j)*vglas(i,8)
935 . +bmerz(i,3,j)*vglas(i,13)+bmkrz(i,3,j)*vglas(i,15))
936 j=4
937 vmz(i,j)= vmz(i,j)
938 . +c2*(vglas(i,17)*bmerz(i,4,j)+vglas(i,18)*bmkrz(i,4,j))
939 . +c3*(bmerz(i,1,j)*vglas(i,1)-bmkrz(i,1,j)*vglas(i,7)
940 . -bmerz(i,2,j)*vglas(i,2)+bmkrz(i,2,j)*vglas(i,8)
941 . +bmerz(i,3,j)*vglas(i,13)+bmkrz(i,3,j)*vglas(i,15))
942C
943 erz = vglas(i,1)*dhgze(i,1)-vglas(i,2)*dhgze(i,2)-
944 . vglas(i,7)*dhgzk(i,1)+vglas(i,8)*dhgzk(i,2)+
945 . vglas(i,13)*dhgze(i,3)+vglas(i,15)*dhgzk(i,3)+
946 . half*(vglas(i,17)*dhgze(i,4)+vglas(i,18)*dhgzk(i,4))
947 c3 = third*area(i)
948C -----------for energy calculation-Mean_value--------
949 esy= ((ss1-ss1_v)*dhg(i,1)+(ss2-ss2_v)*dhg(i,2)+c3*erz)*thk(i)+
950 . fourth*((sc5-sc5_v)*dhg(i,5)+(sc6-sc6_v)*dhg(i,6))
951 etmp(i,1) = half*esy*off(i)
952 emy= (sf1-sf1_v)*dhg(i,3)-(sf2-sf2_v)*dhg(i,4)
953 etmp(i,2) = half*c6(i)*emy*thk(i)*off(i)
954C ----viscous energy is removed to hourglass energy----
955 tesy(i)= (ss1_v*dhg(i,1)+ss2_v*dhg(i,2))*thk(i)+
956 . (sf1_v*dhg(i,3)-sf2_v*dhg(i,4))*thk(i)*c6(i)+
957 . sc5_v*dhg(i,5)+sc6_v*dhg(i,6)
958 ENDDO
959C
960 DO i=jft,jlt
961 eint(i,1) = eint(i,1)+etmp(i,1)
962 eint(i,2) = eint(i,2)+etmp(i,2)
963 END DO
964C
965 ic=1
966 jst(ic)=jft
967 DO j=jft+1,jlt
968 IF (ipartc(j)/=ipartc(j-1)) THEN
969 ic=ic+1
970 jst(ic)=j
971 ENDIF
972 ENDDO
973 jst(ic+1)=jlt+1
974
975 IF(ic==1) THEN
976 mx = ipartc(jft)
977 DO i=jft,jlt
978 evis(8,mx)=evis(8,mx) + tesy(i)
979 ENDDO
980 ELSEIF(ic==2.AND.kfts>0)THEN
981 mx = ipartc(jft)
982 DO i=jft, kfts-1
983 evis(8,mx)=evis(8,mx) + tesy(i)
984 ENDDO
985 mx = ipartc(jlt)
986 DO i=kfts,jlt
987 evis(8,mx)=evis(8,mx) + tesy(i)
988 ENDDO
989 ELSE
990 DO ii=1,ic
991 mx=ipartc(jst(ii))
992 IF (jst(ii+1)-jst(ii)>15) THEN
993 DO j=jst(ii),jst(ii+1)-1
994 evis(8,mx)=evis(8,mx) + tesy(j)
995 ENDDO
996 ELSE
997 DO j=jst(ii),jst(ii+1)-1
998 evis(8,mx)=evis(8,mx) + tesy(j)
999 ENDDO
1000 ENDIF
1001 ENDDO
1002 ENDIF
1003 RETURN
1004 END
1005!||====================================================================
1006!|| czfintn_or ../engine/source/elements/shell/coquez/czfintn.F
1007!||--- called by ------------------------------------------------------
1008!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
1009!|| czforc3_crk ../engine/source/elements/xfem/czforc3_crk.F
1010!||--- calls -----------------------------------------------------
1011!|| czehourv ../engine/source/elements/shell/coquez/czfintn.F
1012!|| czhourep_or ../engine/source/elements/shell/coquez/czfintn.F
1013!||====================================================================
1014 SUBROUTINE czfintn_or(JFT ,JLT ,THK ,C1 ,AA ,VHG ,
1015 2 X13 ,X24 ,Y13 ,Y24 ,Z1 ,MX23 ,
1016 3 MX13 ,MX34 ,MY13 ,MY23 ,MY34 ,VGLAS,
1017 4 VSTRE,MSTRE,VF ,VM ,FAC ,A11 ,
1018 5 A12 ,G ,GS ,SIGY ,OFF ,FAC1 ,
1019 6 RHO ,AREA , DT1 ,EINT ,AMU ,VHI ,
1020 7 NPT ,IPARTC,EVIS,KFTS ,GSR ,
1021 8 A11SR,A12SR ,NUSR,SHFSR,IORTH,HM ,
1022 9 HF ,HC ,HMFOR,MTN ,NEL)
1023C orthotropic elasto-plastic hourglass stresses+linear damping
1024C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1025C CALCUL VF(3,4),VM(2,4) : FORCES INTERNES GENERALISEES
1026C ENTREES : THK,1/AREA,PME,X13,X24,Y13,Y24,MX13,MX34,MY13,MY34
1027C DHG(6) : INCREMENT DE LA DEF HOURGLASS GENERALISEE
1028C VGLAS(12): DEF HOURGLASS GENERALISEE
1029C VGLAS(1-6)->ETA; VGLAS(7-12)->KSI;
1030C FAC=Et/E
1031C SORTIES : VF,VM
1032C VF(3,4),VM(2,4) :
1033C J=1,2 STOCK LA PARTIE ANTISYM (F(1)=-F(3),F(2)=-F(4))
1034C J=3,4 STOCK LA PARTIE SYM (F(1)=F(3),F(2)=F(4))
1035C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1036#include "implicit_f.inc"
1037#include "param_c.inc"
1038#include "mvsiz_p.inc"
1039
1040C-----------------------------------------------
1041C D U M M Y A R G U M E N T S
1042C-----------------------------------------------
1043 INTEGER JFT,JLT,NPT,IPARTC(*),IORTH,MTN,NEL
1044 my_real
1045 . THK(*) ,C1(*) ,AA(*) ,VHG(MVSIZ,6),
1046 . X13(*) ,X24(*) ,Y13(*) ,Y24(*) ,Z1(*),
1047 . MX13(*),MX23(*),MX34(*),MY13(*),MY23(*),MY34(*),
1048 . VGLAS(NEL,12),VSTRE(NEL,5),MSTRE(NEL,3),VF(MVSIZ,3,4),VM(MVSIZ,2,4),
1049 . FAC(MVSIZ,2),A11(*) ,A12(*) ,G(*),GS(*),VHI(MVSIZ,3),SIGY(*),
1050 . FAC1(*) ,RHO(*) ,AREA(*), DT1,OFF(*),EINT(JLT,2),AMU(*) ,
1051 . EVIS(NPSAV,*),HM(MVSIZ,6),HF(MVSIZ,6),HC(MVSIZ,2),HMFOR(MVSIZ,6),
1052 . GSR(*), A11SR(*), A12SR(*), NUSR(*), SHFSR(*)
1053 INTEGER KFTS
1054C-----------------------------------------------
1055C L O C A L V A R I A B L E S
1056C-----------------------------------------------
1057 INTEGER I,MX, IC, II, J
1058 my_real
1059 . STIER,C2,C3,C4,C7,B13,B24,C3M,
1060 . Y13S,X13S,Y34S6,Y23S5,X23S5,X34S6,
1061 . DGLAS(MVSIZ,12),
1062 . C1B,C2B,CMM,CNN,A_BK,A_B,A_K,UNDOUZSR,AUX,
1063 . SSV0,SSV1,SSV2,SSV3,SC6_V,SC5_V,CXX_V,CYY_V,CXY_V,CXZ_V,CYZ_V,
1064 . HXX_V,HYY_V,HXY_V,SS1_V,SS2_V,SS3_V,SF1_V,SF2_V,HVL,CRZ,
1065 . EMY,ECX,ECY,ESY,FBEND,FBEND_V,ERZ,FACM(MVSIZ),FACF(MVSIZ)
1066 my_real
1067 . C5(MVSIZ),CXX_K(MVSIZ),CYY_K(MVSIZ),HXX(MVSIZ),HYY(MVSIZ),
1068 . CXX(MVSIZ),CYY(MVSIZ),HXX_K(MVSIZ),HYY_K(MVSIZ),
1069 . bxx(mvsiz),byy(mvsiz),bxx_k(mvsiz),byy_k(mvsiz),
1070 . c1m(mvsiz),c2m(mvsiz),c6(mvsiz),tesy(mvsiz),esx0(mvsiz),
1071 . ss1(mvsiz),ss2(mvsiz),ss3(mvsiz),sc6(mvsiz),sc5(mvsiz),
1072 . sf1(mvsiz),sf2(mvsiz),hsura(mvsiz),c8(mvsiz),ssz1(mvsiz),
1073 . ssz2(mvsiz),esx(mvsiz),emx(mvsiz),cm(mvsiz,3,3),cf(mvsiz,3,3),
1074 . cmf(mvsiz,3,3),cmh(mvsiz),cfh(mvsiz),cmh_k(mvsiz),cfh_k(mvsiz)
1075C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1076 c7=four_over_3
1077 stier = fivep333
1078 IF (npt==1) THEN
1079 fbend =zero
1080 fbend_v =zero
1081 ELSE
1082 fbend =one_over_12
1083 fbend_v =threep464
1084 ENDIF
1085 undouzsr=sqrt(one_over_12)
1086C------HXX: ETA_i*Y_i,HYY: ETA_i*X_i; HXX_K: Ksi_i*Y_i,HYY_K: Ksi_i*X_i;
1087 DO i=jft,jlt
1088 c3=four*aa(i)
1089 hxx(i)=c3*my34(i)
1090 hyy(i)=c3*mx34(i)
1091 hxx_k(i)=c3*my23(i)
1092 hyy_k(i)=c3*mx23(i)
1093 cxx(i)=hxx(i)*vhg(i,1)
1094 cyy(i)=hyy(i)*vhg(i,2)
1095 cxx_k(i)=hxx_k(i)*vhg(i,1)
1096 cyy_k(i)=hyy_k(i)*vhg(i,2)
1097 bxx(i)=hxx(i)*vhg(i,3)
1098 byy(i)=hyy(i)*vhg(i,4)
1099 bxx_k(i)=hxx_k(i)*vhg(i,3)
1100 byy_k(i)=hyy_k(i)*vhg(i,4)
1101 c5(i)=half*off(i)*thk(i)*c7*dt1
1102 c6(i)=c1(i)*fbend
1103 hsura(i)=thk(i)*aa(i)
1104 END DO
1105C -----------for energy calculation-Mean_value--------
1106 DO i=jft,jlt
1107 ss1(i)= my34(i)*vglas(i,1) +my23(i)*vglas(i,7)
1108 ss2(i)= mx23(i)*vglas(i,8) +mx34(i)*vglas(i,2)
1109 sf1(i)= my34(i)*vglas(i,3) +my23(i)*vglas(i,9)
1110 sf2(i)= -mx23(i)*vglas(i,10)-mx34(i)*vglas(i,4)
1111 sc5(i)= my34(i)*vglas(i,5) +mx34(i)*vglas(i,6)
1112 sc6(i)= my23(i)*vglas(i,11)+mx23(i)*vglas(i,12)
1113 esx0(i)= ss1(i)*vhg(i,1)+ss2(i)*vhg(i,2)
1114 esx(i)= esx0(i)+
1115 . fourth*(sc5(i)*vhg(i,5)+sc6(i)*vhg(i,6))
1116 emx(i)= (sf1(i)*vhg(i,3)-sf2(i)*vhg(i,4))*c6(i)
1117 END DO
1118C -----------for energy calculation-Mean_value--VGLAS(t)------
1119 DO i=jft,jlt
1120 eint(i,1) = eint(i,1)+c5(i)*esx(i)
1121 eint(i,2) = eint(i,2)+c5(i)*emx(i)
1122 END DO
1123C -----------elastic increment
1124 IF (iorth>0) THEN
1125 IF ((mtn==19).OR.(mtn==119)) THEN
1126 DO i=jft,jlt
1127 c2=fac1(i)*dt1
1128 facm(i)=fac(i,1)*c2
1129 facf(i)=fac(i,2)*c2
1130 END DO
1131 ELSE
1132C---------coef of integration should be removed, because has been done inside HF
1133C---------have to be done here for VGLAS
1134 DO i=jft,jlt
1135 c2=fac1(i)*dt1
1136 facm(i)=c2
1137 facf(i)=c2*twelve
1138 END DO
1139 END IF !(MTN==19) THEN
1140C-------Modified modulas for damping---
1141 DO i=jft,jlt
1142 gsr(i) =sqrt(hm(i,4))
1143 a11sr(i)=sqrt(half*(hm(i,1)+hm(i,2)))
1144 a12sr(i)=sqrt(max(em20,hm(i,3)))
1145cc A12SR(I)=SQRT(HM(3,I))
1146 END DO
1147 DO i=jft,jlt
1148 cm(i,1,1)=facm(i)*hm(i,1)
1149 cm(i,2,2)=facm(i)*hm(i,2)
1150 cm(i,1,2)=facm(i)*hm(i,3)
1151 cm(i,3,3)=facm(i)*hm(i,4)
1152 cm(i,1,3)=facm(i)*hm(i,5)
1153 cm(i,2,3)=facm(i)*hm(i,6)
1154 cf(i,1,1)=facf(i)*hf(i,1)
1155 cf(i,2,2)=facf(i)*hf(i,2)
1156 cf(i,1,2)=facf(i)*hf(i,3)
1157 cf(i,3,3)=facf(i)*hf(i,4)
1158 cf(i,1,3)=facf(i)*hf(i,5)
1159 cf(i,2,3)=facf(i)*hf(i,6)
1160 END DO
1161C VGLAS(1-6)->ETA (x:+ y-); VGLAS(7-12)->KSI (x:- y+);
1162C-----part iso-commun
1163 DO i=jft,jlt
1164C----- x_ETA, signe +
1165 dglas(i,1) =cm(i,1,1)*cxx(i)-cm(i,1,2)*cyy(i)
1166C----- y_ETA, signe -
1167 dglas(i,2) =cm(i,2,2)*cyy(i)-cm(i,1,2)*cxx(i)
1168C----- x_ETA, signe+ bending
1169 dglas(i,3) =cf(i,1,1)*bxx(i)-cf(i,1,2)*byy(i)
1170C----- y_ETA, signe - bending
1171 dglas(i,4) =cf(i,2,2)*byy(i)-cf(i,1,2)*bxx(i)
1172C----- x_Ksi, signe -
1173 dglas(i,7) =cm(i,1,1)*cxx_k(i)-cm(i,1,2)*cyy_k(i)
1174C----- y_Ksi, signe +
1175 dglas(i,8) =cm(i,2,2)*cyy_k(i)-cm(i,1,2)*cxx_k(i)
1176C----- x_Ksi, signe - bending
1177 dglas(i,9) =cf(i,1,1)*bxx_k(i)-cf(i,1,2)*byy_k(i)
1178C----- y_Ksi, signe + bending
1179 dglas(i,10)=cf(i,2,2)*byy_k(i)-cf(i,1,2)*bxx_k(i)
1180C------------
1181 c4=fac1(i)*one_over_64*dt1
1182 c2=c4*hc(i,1)
1183 c3=c4*hc(i,2)
1184 dglas(i,5) =c2*hxx(i)*vhg(i,5)
1185 dglas(i,6) =c2*hyy(i)*vhg(i,5)
1186 dglas(i,11)=c3*hxx_k(i)*vhg(i,6)
1187 dglas(i,12)=c3*hyy_k(i)*vhg(i,6)
1188 END DO
1189C-----part C(1,3),C(2,3)
1190C-------finally after crash-box tets, constant in-plane shear strain and stress
1191C--------- like QBAT is used to avoid instability
1192 DO i=jft,jlt
1193 cmh(i) = zero
1194 cmh_k(i) = zero
1195 cfh(i) = zero
1196 cfh_k(i) = zero
1197 END DO
1198C-----part CMF
1199 IF (iorth>1) THEN
1200 DO i=jft,jlt
1201 c2=em01*fac1(i)*dt1
1202 cmf(i,1,1)=c2*hmfor(i,1)
1203 cmf(i,2,2)=c2*hmfor(i,2)
1204 cmf(i,1,2)=c2*hmfor(i,3)
1205 cmf(i,3,3)=c2*hmfor(i,4)
1206 cmf(i,1,3)=c2*hmfor(i,5)
1207 cmf(i,2,3)=c2*hmfor(i,6)
1208C
1209 dglas(i,1) =dglas(i,1)+cmf(i,1,1)*bxx(i)-cmf(i,1,2)*byy(i)
1210 + -cmf(i,1,3)*cfh(i)
1211 dglas(i,2) =dglas(i,2)+cmf(i,2,2)*byy(i)-cmf(i,1,2)*bxx(i)
1212 + +cmf(i,2,3)*cfh(i)
1213 dglas(i,3) =dglas(i,3)+cmf(i,1,1)*cxx(i)-cmf(i,1,2)*cyy(i)
1214 + -cmf(i,1,3)*cmh(i)
1215 dglas(i,4) =dglas(i,4)+cmf(i,2,2)*cyy(i)-cmf(i,1,2)*cxx(i)
1216 + +cmf(i,2,3)*cmh(i)
1217 dglas(i,7) =dglas(i,7)+cmf(i,1,1)*bxx_k(i)-cmf(i,1,2)*byy_k(i)
1218 + -cmf(i,1,3)*cfh_k(i)
1219 dglas(i,8) =dglas(i,8)+cmf(i,2,2)*byy_k(i)-cmf(i,1,2)*bxx_k(i)
1220 + +cmf(i,2,3)*cfh_k(i)
1221 dglas(i,9) =dglas(i,9)+cmf(i,1,1)*cxx_k(i)-cmf(i,1,2)*cyy_k(i)
1222 + -cmf(i,1,3)*cmh_k(i)
1223 dglas(i,10)=dglas(i,10)+cmf(i,2,2)*cyy_k(i)-cmf(i,1,2)*cxx_k(i)
1224 + +cmf(i,2,3)*cmh_k(i)
1225 END DO
1226 END IF !(IORTH>1) THEN
1227 ELSE
1228 DO i=jft,jlt
1229 c1m(i)=a11(i)*fac1(i)*dt1
1230 c2m(i)=a12(i)*fac1(i)*dt1
1231 END DO
1232 DO i=jft,jlt
1233 dglas(i,1) =c1m(i)*cxx(i)-c2m(i)*cyy(i)
1234 dglas(i,2) =c1m(i)*cyy(i)-c2m(i)*cxx(i)
1235 dglas(i,3) =c1m(i)*bxx(i)-c2m(i)*byy(i)
1236 dglas(i,4) =c1m(i)*byy(i)-c2m(i)*bxx(i)
1237 dglas(i,7) =c1m(i)*cxx_k(i)-c2m(i)*cyy_k(i)
1238 dglas(i,8) =c1m(i)*cyy_k(i)-c2m(i)*cxx_k(i)
1239 dglas(i,9) =c1m(i)*bxx_k(i)-c2m(i)*byy_k(i)
1240 dglas(i,10)=c1m(i)*byy_k(i)-c2m(i)*bxx_k(i)
1241C------------
1242 c2=fac1(i)*gs(i)*one_over_64*dt1
1243 dglas(i,5) =c2*hxx(i)*vhg(i,5)
1244 dglas(i,6) =c2*hyy(i)*vhg(i,5)
1245 dglas(i,11)=c2*hxx_k(i)*vhg(i,6)
1246 dglas(i,12)=c2*hyy_k(i)*vhg(i,6)
1247 END DO
1248 END IF !(IORTH>0) THEN
1249 DO i=jft,jlt
1250 DO j=1,12
1251 vglas(i,j) =vglas(i,j) +dglas(i,j)
1252 END DO
1253 END DO
1254C-----------elasto-plastic----
1255 CALL czhourep_or(jft ,jlt ,sigy ,vstre ,mstre ,
1256 1 thk ,fac ,esx0 ,emx ,npt ,
1257 2 dglas ,vglas,nel )
1258C-----------to use the same coef C5
1259 DO i=jft,jlt
1260 ss1(i)= my34(i)*vglas(i,1) +my23(i)*vglas(i,7)
1261 ss2(i)= mx23(i)*vglas(i,8) +mx34(i)*vglas(i,2)
1262 sf1(i)= my34(i)*vglas(i,3) +my23(i)*vglas(i,9)
1263 sf2(i)= -(mx23(i)*vglas(i,10)+mx34(i)*vglas(i,4))
1264 sc5(i)=my34(i)*vglas(i,5) +mx34(i)*vglas(i,6)
1265 sc6(i)=my23(i)*vglas(i,11)+mx23(i)*vglas(i,12)
1266 END DO !I=JFT,JLT
1267C -------for energy calculation-Mean_value-VGLAS(t+dt)-------
1268 DO i=jft,jlt
1269 esx(i)= ss1(i)*vhg(i,1)+ss2(i)*vhg(i,2)+
1270 . fourth*(sc5(i)*vhg(i,5)+sc6(i)*vhg(i,6))
1271 emx(i)= (sf1(i)*vhg(i,3)-sf2(i)*vhg(i,4))*c6(i)
1272 END DO
1273 DO i=jft,jlt
1274 eint(i,1) = eint(i,1)+c5(i)*esx(i)
1275 eint(i,2) = eint(i,2)+c5(i)*emx(i)
1276 END DO !I=JFT,JLT
1277 DO i=jft,jlt
1278 c8(i) = c7*off(i)
1279 ss1(i)= ss1(i)*c8(i)
1280 ss2(i)= ss2(i)*c8(i)
1281 sf1(i)= sf1(i)*c8(i)
1282 sf2(i)= sf2(i)*c8(i)
1283 c2 = c8(i)*thk(i)
1284 sc5(i)=sc5(i)*c2
1285 sc6(i)=sc6(i)*c2
1286 ss3(i)=sc5(i)+sc6(i)
1287 END DO !I=JFT,JLT
1288C--------Fint :terms for linear damping
1289 DO i=jft,jlt
1290C-------save for Fz--and add coupling terms(warped)------
1291 ssz1(i) = ss1(i)
1292 ssz2(i) = ss2(i)
1293 hvl = amu(i)*sqrt(rho(i)*area(i)*fac1(i))*off(i)
1294 ssv0=my23(i)*my23(i)
1295 ssv1=my34(i)*my34(i)
1296 ssv2=mx23(i)*mx23(i)
1297 ssv3=mx34(i)*mx34(i)
1298 hxx_v= stier*(ssv1+ssv0)
1299 hxy_v=-stier*(my34(i)*mx34(i)+my23(i)*mx23(i))
1300 hyy_v= stier*(ssv2+ssv3)
1301 c2 =hvl*gsr(i)*shfsr(i)*undouzsr
1302 cxz_v=(ssv1+ssv3)*c2
1303 cyz_v=(ssv2+ssv0)*c2
1304C
1305 aux=aa(i)*hvl
1306 c1m(i) = a11sr(i)*aux
1307 c2m(i) = a12sr(i)*aux
1308C
1309 cxx_v=c1m(i)*hxx_v
1310 cyy_v=c1m(i)*hyy_v
1311 cxy_v=c2m(i)*hxy_v
1312 ss1_v=cxx_v*vhg(i,1)+cxy_v*vhg(i,2)
1313 ss2_v=cyy_v*vhg(i,2)+cxy_v*vhg(i,1)
1314 sf1_v= (cxx_v*vhg(i,3)+cxy_v*vhg(i,4))*fbend_v
1315 sf2_v=(-cyy_v*vhg(i,4)-cxy_v*vhg(i,3))*fbend_v
1316 sc5_v=cxz_v*vhg(i,5)*hsura(i)
1317 sc6_v=cyz_v*vhg(i,6)*hsura(i)
1318 ss3_v = sc5_v+sc6_v
1319C--------------------------------------------------
1320 ss1(i)=ss1(i)+ss1_v
1321 ss2(i)=ss2(i)+ss2_v
1322 ss3(i)=ss3(i)+ss3_v
1323 sc5(i)=sc5(i)+sc5_v
1324 sc6(i)=sc6(i)+sc6_v
1325 sf1(i)=sf1(i)+sf1_v
1326 sf2(i)=sf2(i)+sf2_v
1327 tesy(i)= ((ss1_v*vhg(i,1)+ss2_v*vhg(i,2))*thk(i)+
1328 . (sf1_v*vhg(i,3)-sf2_v*vhg(i,4))*thk(i)*c6(i) +
1329 . (sc5_v*vhg(i,5)+sc6_v*vhg(i,6))*fourth)*dt1
1330 END DO !I=JFT,JLT
1331C--------Fint
1332 DO i=jft,jlt
1333 y13s= my13(i)*ss3(i)
1334 x13s= mx13(i)*ss3(i)
1335 y34s6=my34(i)*sc6(i)
1336 y23s5=my23(i)*sc5(i)
1337 x23s5=mx23(i)*sc5(i)
1338 x34s6=mx34(i)*sc6(i)
1339 c2=fourth*thk(i)
1340C--------NOEUDS 1,3-------
1341 b13=(my13(i)*x24(i)-mx13(i)*y24(i))*hsura(i)
1342 vf(i,1,1)=vf(i,1,1)+b13*ss1(i)
1343 vf(i,1,3)=c2*ss1(i)
1344 vf(i,2,1)=vf(i,2,1)+b13*ss2(i)
1345 vf(i,2,3)=c2*ss2(i)
1346 vf(i,3,3)=ss3(i)
1347C--------NOEUDS 2,4-------
1348 b24=(mx13(i)*y13(i)-my13(i)*x13(i))*hsura(i)
1349 vf(i,1,2)= vf(i,1,2)+b24*ss1(i)
1350 vf(i,1,4)=-vf(i,1,3)
1351 vf(i,2,2)= vf(i,2,2)+b24*ss2(i)
1352 vf(i,2,4)=-vf(i,2,3)
1353 vf(i,3,4)=-vf(i,3,3)
1354C--------NOEUDS 1,3-------
1355 c3=c6(i)*b13
1356 c4=c6(i)*c2
1357 vm(i,1,1)=vm(i,1,1)+c3*sf2(i)+y23s5+y34s6
1358 vm(i,1,3)=vm(i,1,3)+c4*sf2(i)-y13s
1359 vm(i,2,1)=vm(i,2,1)+c3*sf1(i)-x23s5-x34s6
1360 vm(i,2,3)=vm(i,2,3)+c4*sf1(i)+x13s
1361C--------NOEUDS 2,4-------
1362 c3=c6(i)*b24
1363 vm(i,1,2)=vm(i,1,2)+c3*sf2(i)+y23s5-y34s6
1364 vm(i,1,4)=vm(i,1,4)-c4*sf2(i)-y13s
1365 vm(i,2,2)=vm(i,2,2)+c3*sf1(i)-x23s5+x34s6
1366 vm(i,2,4)=vm(i,2,4)-c4*sf1(i)+x13s
1367C
1368 c2=z1(i)*hsura(i)
1369 vf(i,3,1)=vf(i,3,1)+c2*(ssz1(i)*y24(i)-ssz2(i)*x24(i))
1370 vf(i,3,2)=vf(i,3,2)+c2*(-ssz1(i)*y13(i)+ssz2(i)*x13(i))
1371 END DO !I=JFT,JLT
1372C
1373 CALL czehourv(jft ,jlt ,ipartc,evis,kfts ,tesy )
1374 RETURN
1375 END
1376!||====================================================================
1377!|| czehourv ../engine/source/elements/shell/coquez/czfintn.F
1378!||--- called by ------------------------------------------------------
1379!|| czfintn_or ../engine/source/elements/shell/coquez/czfintn.F
1380!|| czfintnrz_or ../engine/source/elements/shell/coquez/czfintn.F
1381!||====================================================================
1382 SUBROUTINE czehourv(JFT ,JLT ,IPARTC,EVIS,KFTS ,
1383 8 TESY )
1384C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1385#include "implicit_f.inc"
1386#include "param_c.inc"
1387#include "mvsiz_p.inc"
1388C-----------------------------------------------
1389C D U M M Y A R G U M E N T S
1390C-----------------------------------------------
1391 INTEGER JFT,JLT,IPARTC(*)
1392 my_real
1393 . evis(npsav,*),tesy(*)
1394 INTEGER KFTS
1395C-----------------------------------------------
1396C L O C A L V A R I A B L E S
1397C-----------------------------------------------
1398 INTEGER I,MX, IC, II, J, JST(MVSIZ+1)
1399C-----------------------------------------------
1400 ic=1
1401 jst(ic)=jft
1402 DO j=jft+1,jlt
1403 IF (ipartc(j)/=ipartc(j-1)) THEN
1404 ic=ic+1
1405 jst(ic)=j
1406 ENDIF
1407 ENDDO
1408 jst(ic+1)=jlt+1
1409 IF(ic==1) THEN
1410 mx = ipartc(jft)
1411 DO i=jft,jlt
1412 evis(8,mx)=evis(8,mx) + tesy(i)
1413 ENDDO
1414
1415 ELSEIF(ic==2.AND.kfts>0)THEN
1416 mx = ipartc(jft)
1417 DO i=jft, kfts-1
1418 evis(8,mx)=evis(8,mx) + tesy(i)
1419 ENDDO
1420 mx = ipartc(jlt)
1421 DO i=kfts,jlt
1422 evis(8,mx)=evis(8,mx) + tesy(i)
1423 ENDDO
1424
1425 ELSE
1426C
1427 DO ii=1,ic
1428 mx=ipartc(jst(ii))
1429 IF (jst(ii+1)-jst(ii)>15) THEN
1430 DO j=jst(ii),jst(ii+1)-1
1431 evis(8,mx)=evis(8,mx) + tesy(j)
1432 ENDDO
1433 ELSE
1434 DO j=jst(ii),jst(ii+1)-1
1435 evis(8,mx)=evis(8,mx) + tesy(j)
1436 ENDDO
1437 ENDIF
1438 ENDDO
1439 ENDIF
1440C
1441 RETURN
1442 END
1443!||====================================================================
1444!|| czfintnrz_or ../engine/source/elements/shell/coquez/czfintn.F
1445!||--- called by ------------------------------------------------------
1446!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
1447!|| czforc3_crk ../engine/source/elements/xfem/czforc3_crk.F
1448!||--- calls -----------------------------------------------------
1449!|| czehourv ../engine/source/elements/shell/coquez/czfintn.F
1450!|| czhoureprz_or ../engine/source/elements/shell/coquez/czfintn.F
1451!||====================================================================
1452 SUBROUTINE czfintnrz_or(JFT ,JLT ,THK ,C1 ,AA ,VHG ,
1453 2 X13 ,X24 ,Y13 ,Y24 ,Z1 ,MX23 ,
1454 3 MX13 ,MX34 ,MY13 ,MY23 ,MY34 ,VGLAS,
1455 4 VSTRE,MSTRE,VF ,VM ,FAC ,A11 ,
1456 5 A12 ,G ,GS ,SIGY ,OFF ,FAC1 ,
1457 6 RHO ,AREA , DT1,EINT ,AMU ,VHI ,
1458 7 NPT ,IPARTC,EVIS,KFTS ,GSR ,
1459 8 A11SR,A12SR ,NUSR,SHFSR,BMKRZ,BMERZ,
1460 9 VHGZK,VHGZE ,KRZ ,VMZ ,IORTH,HM ,
1461 9 HF ,HC ,HMFOR,MTN ,NEL)
1462C orthotropic elasto-plastic hourglass stresses+linear damping with drilling dof
1463C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1464C CALCUL VF(3,4),VM(2,4) : FORCES INTERNES GENERALISEES
1465C ENTREES : THK,1/AREA,PME,X13,X24,Y13,Y24,MX13,MX34,MY13,MY34
1466C DHG(6) : INCREMENT DE LA DEF HOURGLASS GENERALISEE
1467C VGLAS(12): DEF HOURGLASS GENERALISEE
1468C VGLAS(1-6)->ETA; VGLAS(7-12)->KSI;
1469C -shear: VGLAS(13-14)->ETA; VGLAS(15-16)->KSI;
1470C -Sig_rz: VGLAS(17)->ETA; VGLAS(18)->KSI; VGLAS(19)->const
1471C FAC=Et/E
1472C SORTIES : VF,VM
1473C VF(3,4),VM(2,4) :
1474C J=1,2 STOCK LA PARTIE ANTISYM (F(1)=-F(3),F(2)=-F(4))
1475C J=3,4 STOCK LA PARTIE SYM (F(1)=F(3),F(2)=F(4))
1476C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1477#include "implicit_f.inc"
1478#include "param_c.inc"
1479#include "mvsiz_p.inc"
1480
1481C-----------------------------------------------
1482C D U M M Y A R G U M E N T S
1483C-----------------------------------------------
1484 INTEGER JFT,JLT,NPT,IPARTC(*),IORTH,MTN,NEL
1485 my_real
1486 . THK(*) ,C1(*) ,AA(*) ,VHG(MVSIZ,6),
1487 . X13(*) ,X24(*) ,Y13(*) ,Y24(*) ,Z1(*),
1488 . mx13(*),mx23(*),mx34(*),my13(*),my23(*),my34(*),
1489 . vglas(nel,19),vstre(nel,5),mstre(nel,3),vf(mvsiz,3,4),vm(mvsiz,2,4),
1490 . fac(mvsiz,2),a11(*) ,a12(*) ,g(*),gs(*),vhi(mvsiz,3),sigy(*),
1491 . fac1(*) ,rho(*) ,area(*), dt1,off(*),eint(jlt,2),amu(*) ,
1492 . evis(npsav,*),vmz(mvsiz,4),
1493 . vhgzk(mvsiz,5),vhgze(mvsiz,5),krz(*),
1494 . bm0rz(mvsiz,4,4),bmkrz(mvsiz,4,4),bmerz(mvsiz,4,4),
1495 . gsr(*), a11sr(*), a12sr(*), nusr(*), shfsr(*),
1496 . hm(mvsiz,6),hf(mvsiz,6),hc(mvsiz,2),hmfor(mvsiz,6)
1497 INTEGER KFTS
1498C-----------------------------------------------
1499C L O C A L V A R I A B L E S
1500C-----------------------------------------------
1501 INTEGER I,MX, IC, II, J
1502 my_real
1503 . STIER,C2,C3,C4,C7,B13,B24,C3M,
1504 . Y13S,X13S,Y34S6,Y23S5,X23S5,X34S6,
1505 . DGLAS(MVSIZ,18),FACM(MVSIZ),FACF(MVSIZ),
1506 . C1B,C2B,CMM,CNN,A_BK,A_B,A_K,UNDOUZSR,AUX,
1507 . SSV0,SSV1,SSV2,SSV3,SC6_V,SC5_V,CXX_V,CYY_V,CXY_V,CXZ_V,CYZ_V,
1508 . HXX_V,HYY_V,HXY_V,SS1_V,SS2_V,SS3_V,SF1_V,SF2_V,HVL,CRZ,
1509 . EMY,ECX,ECY,ESY,FBEND,FBEND_V,ERZ
1510 my_real
1511 . C5(MVSIZ),CXX_K(MVSIZ),CYY_K(MVSIZ),HXX(MVSIZ),HYY(MVSIZ),
1512 . CXX(MVSIZ),CYY(MVSIZ),HXX_K(MVSIZ),HYY_K(MVSIZ),
1513 . BXX(MVSIZ),BYY(MVSIZ),BXX_K(MVSIZ),BYY_K(MVSIZ),
1514 . C1M(MVSIZ),C2M(MVSIZ),C6(MVSIZ),TESY(MVSIZ),ESX0(MVSIZ),
1515 . SS1(MVSIZ),SS2(MVSIZ),SS3(MVSIZ),SC6(MVSIZ),SC5(MVSIZ),
1516 . SF1(MVSIZ),SF2(MVSIZ),HSURA(MVSIZ),C8(MVSIZ),SSZ1(MVSIZ),
1517 . SSZ2(MVSIZ),ESX(MVSIZ),EMX(MVSIZ),CM(MVSIZ,3,3),CF(MVSIZ,3,3),
1518 . CMF(MVSIZ,3,3),CMH(MVSIZ),CFH(MVSIZ),CMH_K(MVSIZ),CFH_K(MVSIZ)
1519C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1520 c7=four_over_3
1521 stier = fivep333
1522 IF (npt==1) THEN
1523 fbend =zero
1524 fbend_v =zero
1525 ELSE
1526 fbend =one_over_12
1527 fbend_v =threep464
1528 ENDIF
1529 undouzsr=sqrt(one_over_12)
1530C
1531 DO i=jft,jlt
1532 c3=four*aa(i)
1533 hxx(i)=c3*my34(i)
1534 hyy(i)=c3*mx34(i)
1535 hxx_k(i)=c3*my23(i)
1536 hyy_k(i)=c3*mx23(i)
1537 cxx(i)=hxx(i)*vhg(i,1)
1538 cyy(i)=hyy(i)*vhg(i,2)
1539 cxx_k(i)=hxx_k(i)*vhg(i,1)
1540 cyy_k(i)=hyy_k(i)*vhg(i,2)
1541 bxx(i)=hxx(i)*vhg(i,3)
1542 byy(i)=hyy(i)*vhg(i,4)
1543 bxx_k(i)=hxx_k(i)*vhg(i,3)
1544 byy_k(i)=hyy_k(i)*vhg(i,4)
1545 c5(i)=half*off(i)*thk(i)*c7*dt1
1546 c6(i)=c1(i)*fbend
1547 hsura(i)=thk(i)*aa(i)
1548 END DO
1549C -----------for energy calculation-Mean_value--------
1550 DO i=jft,jlt
1551 ss1(i)= my34(i)*vglas(i,1) +my23(i)*vglas(i,7)
1552 ss2(i)= mx23(i)*vglas(i,8) +mx34(i)*vglas(i,2)
1553 sf1(i)= my34(i)*vglas(i,3) +my23(i)*vglas(i,9)
1554 sf2(i)= -mx23(i)*vglas(i,10)-mx34(i)*vglas(i,4)
1555 sc5(i)= my34(i)*vglas(i,5) +mx34(i)*vglas(i,6)
1556 sc6(i)= my23(i)*vglas(i,11)+mx23(i)*vglas(i,12)
1557 esx0(i)= ss1(i)*vhg(i,1)+ss2(i)*vhg(i,2)
1558 esx(i)= esx0(i)+
1559 . fourth*(sc5(i)*vhg(i,5)+sc6(i)*vhg(i,6))
1560 emx(i)= (sf1(i)*vhg(i,3)-sf2(i)*vhg(i,4))*c6(i)
1561 END DO
1562C
1563 DO i=jft,jlt
1564 ss1(i)=-mx34(i)*vglas(i,13)+mx23(i)*vglas(i,15)+
1565 + half*(mx34(i)*vglas(i,17)-mx23(i)*vglas(i,18))
1566 ss2(i)=my34(i)*vglas(i,13) -my23(i)*vglas(i,15)+
1567 + half*(my34(i)*vglas(i,17)-my23(i)*vglas(i,18))
1568 sf1(i)=(-mx34(i)*vglas(i,14)+mx23(i)*vglas(i,16))
1569 sf2(i)=-(my34(i)*vglas(i,14)-my23(i)*vglas(i,16))
1570 erz = vglas(i,1)*vhgze(i,1)-vglas(i,2)*vhgze(i,2)-
1571 . vglas(i,7)*vhgzk(i,1)+vglas(i,8)*vhgzk(i,2)+
1572 . vglas(i,13)*vhgze(i,3)+vglas(i,15)*vhgzk(i,3)+
1573 . half*(vglas(i,17)*vhgze(i,4)+vglas(i,18)*vhgzk(i,4))
1574C--------
1575 esx(i)=esx(i)+
1576 + ss1(i)*vhg(i,1)+ss2(i)*vhg(i,2)+erz*area(i)*fourth
1577 emx(i)=emx(i)+(sf1(i)*vhg(i,3)-sf2(i)*vhg(i,4))*c6(i)
1578 END DO
1579C -----------for energy calculation-Mean_value--VGLAS(t)------
1580 DO i=jft,jlt
1581 eint(i,1) = eint(i,1)+c5(i)*esx(i)
1582 eint(i,2) = eint(i,2)+c5(i)*emx(i)
1583 END DO
1584C -----------elastic increment
1585 IF (iorth>0) THEN
1586 IF ((mtn==19).OR.(mtn==119)) THEN
1587 DO i=jft,jlt
1588 c2=fac1(i)*dt1
1589 facm(i)=fac(i,1)*c2
1590 facf(i)=fac(i,2)*c2
1591 END DO
1592 ELSE
1593 DO i=jft,jlt
1594 c2=fac1(i)*dt1
1595 facm(i)=c2
1596 facf(i)=c2*twelve
1597 END DO
1598 END IF !(MTN==19) THEN
1599C-------Modified modulas for damping---
1600 DO i=jft,jlt
1601 gsr(i) =sqrt(hm(i,4))
1602 a11sr(i)=sqrt(half*(hm(i,1)+hm(i,2)))
1603 a12sr(i)=sqrt(hm(i,3))
1604 END DO
1605 DO i=jft,jlt
1606 cm(i,1,1)=facm(i)*hm(i,1)
1607 cm(i,2,2)=facm(i)*hm(i,2)
1608 cm(i,1,2)=facm(i)*hm(i,3)
1609 cm(i,3,3)=facm(i)*hm(i,4)
1610 cm(i,1,3)=facm(i)*hm(i,5)
1611 cm(i,2,3)=facm(i)*hm(i,6)
1612 cf(i,1,1)=facf(i)*hf(i,1)
1613 cf(i,2,2)=facf(i)*hf(i,2)
1614 cf(i,1,2)=facf(i)*hf(i,3)
1615 cf(i,3,3)=facf(i)*hf(i,4)
1616 cf(i,1,3)=facf(i)*hf(i,5)
1617 cf(i,2,3)=facf(i)*hf(i,6)
1618 END DO
1619 DO i=jft,jlt
1620 cmh(i) = hyy(i)*vhg(i,1)-hxx(i)*vhg(i,2)
1621 cfh(i) = hyy(i)*vhg(i,3)-hxx(i)*vhg(i,4)
1622 cmh_k(i) = hyy_k(i)*vhg(i,1)-hxx_k(i)*vhg(i,2)
1623 cfh_k(i) = hyy_k(i)*vhg(i,3)-hxx_k(i)*vhg(i,4)
1624 END DO
1625C VGLAS(1-6)->ETA (x:+ y-); VGLAS(7-12)->KSI (x:- y+);
1626C-----part iso-commun
1627 DO i=jft,jlt
1628C----- x_ETA, signe +
1629 dglas(i,1) =cm(i,1,1)*cxx(i)-cm(i,1,2)*cyy(i)
1630C----- y_ETA, signe -
1631 dglas(i,2) =cm(i,2,2)*cyy(i)-cm(i,1,2)*cxx(i)
1632C----- x_ETA, signe + bending
1633 dglas(i,3) =cf(i,1,1)*bxx(i)-cf(i,1,2)*byy(i)
1634C----- y_ETA, signe - bending
1635 dglas(i,4) =cf(i,2,2)*byy(i)-cf(i,1,2)*bxx(i)
1636C----- x_Ksi, signe -
1637 dglas(i,7) =cm(i,1,1)*cxx_k(i)-cm(i,1,2)*cyy_k(i)
1638C----- y_Ksi, signe +
1639 dglas(i,8) =cm(i,2,2)*cyy_k(i)-cm(i,1,2)*cxx_k(i)
1640C----- x_Ksi, signe - bending
1641 dglas(i,9) =cf(i,1,1)*bxx_k(i)-cf(i,1,2)*byy_k(i)
1642C----- y_Ksi, signe + bending
1643 dglas(i,10)=cf(i,2,2)*byy_k(i)-cf(i,1,2)*bxx_k(i)
1644C------------
1645 c4=fac1(i)*one_over_64*dt1
1646 c2=c4*hc(i,1)
1647 c3=c4*hc(i,2)
1648 dglas(i,5) =c2*hxx(i)*vhg(i,5)
1649 dglas(i,6) =c2*hyy(i)*vhg(i,5)
1650 dglas(i,11)=c3*hxx_k(i)*vhg(i,6)
1651 dglas(i,12)=c3*hyy_k(i)*vhg(i,6)
1652 c3m=cm(i,3,3)
1653 crz=half*krz(i)*fac1(i)*dt1
1654C --------add shear due to rz dof--no sign change-------
1655C----- xy_ETA, signe +
1656 dglas(i,13)=c3m*(-cmh(i)+vhgze(i,3))
1657C----- xy_ETA, signe + bending
1658 dglas(i,14)=-c3m*cfh(i)
1659C----- xy_Ksi, signe +
1660 dglas(i,15)=c3m*(cmh_k(i)+vhgzk(i,3))
1661C----- xy_Ksi, signe + bending
1662 dglas(i,16)=c3m*cfh_k(i)
1663C --------add srz dof--------
1664 dglas(i,17)=crz*vhgze(i,5)
1665 dglas(i,18)=crz*vhgzk(i,5)
1666 END DO
1667C-----part C(1,3),C(2,3)
1668 DO i=jft,jlt
1669C
1670 dglas(i,1) =dglas(i,1)-cm(i,1,3)*cmh(i)
1671 dglas(i,2) =dglas(i,2)+cm(i,2,3)*cmh(i)
1672 dglas(i,3) =dglas(i,3)-cf(i,1,3)*cfh(i)
1673 dglas(i,4) =dglas(i,4)+cf(i,2,3)*cfh(i)
1674 dglas(i,7) =dglas(i,7)-cm(i,1,3)*cmh_k(i)
1675 dglas(i,8) =dglas(i,8)+cm(i,2,3)*cmh_k(i)
1676 dglas(i,9) =dglas(i,9)-cf(i,1,3)*cfh_k(i)
1677 dglas(i,10)=dglas(i,10)+cf(i,2,3)*cfh_k(i)
1678C --------add shear due to rz dof--no sign change-------
1679 dglas(i,13)=dglas(i,13)+
1680 + cm(i,1,3)*hxx(i)*vhg(i,1)-cm(i,2,3)*hyy(i)*vhg(i,2)
1681 dglas(i,14)=dglas(i,14)+
1682 + cf(i,1,3)*hxx(i)*vhg(i,3)-cf(i,2,3)*hyy(i)*vhg(i,4)
1683 dglas(i,15)=dglas(i,15)
1684 + -cm(i,1,3)*hxx_k(i)*vhg(i,1)+cm(i,2,3)*hyy_k(i)*vhg(i,2)
1685 dglas(i,16)=dglas(i,16)
1686 + -cf(i,1,3)*hxx_k(i)*vhg(i,3)+cf(i,2,3)*hyy_k(i)*vhg(i,4)
1687 END DO
1688C-----part CMF
1689 IF (iorth>1) THEN
1690 DO i=jft,jlt
1691 c2=fac1(i)*dt1
1692 cmf(i,1,1)=c2*hmfor(i,1)
1693 cmf(i,2,2)=c2*hmfor(i,2)
1694 cmf(i,1,2)=c2*hmfor(i,3)
1695 cmf(i,3,3)=c2*hmfor(i,4)
1696 cmf(i,1,3)=c2*hmfor(i,5)
1697 cmf(i,2,3)=c2*hmfor(i,6)
1698C
1699 dglas(i,1) =dglas(i,1)+cmf(i,1,1)*bxx(i)-cmf(i,1,2)*byy(i)
1700 + -cmf(i,1,3)*cfh(i)
1701 dglas(i,2) =dglas(i,2)+cmf(i,2,2)*byy(i)-cmf(i,1,2)*bxx(i)
1702 + +cmf(i,2,3)*cfh(i)
1703 dglas(i,3) =dglas(i,3)+cmf(i,1,1)*cxx(i)-cmf(i,1,2)*cyy(i)
1704 + -cmf(i,1,3)*cmh(i)
1705 dglas(i,4) =dglas(i,4)+cmf(i,2,2)*cyy(i)-cmf(i,1,2)*cxx(i)
1706 + +cmf(i,2,3)*cmh(i)
1707 dglas(i,7) =dglas(i,7)+cmf(i,1,1)*bxx_k(i)-cmf(i,1,2)*byy_k(i)
1708 + -cmf(i,1,3)*cfh_k(i)
1709 dglas(i,8) =dglas(i,8)+cmf(i,2,2)*byy_k(i)-cmf(i,1,2)*bxx_k(i)
1710 + +cmf(i,2,3)*cfh_k(i)
1711 dglas(i,9) =dglas(i,9)+cmf(i,1,1)*cxx_k(i)-cmf(i,1,2)*cyy_k(i)
1712 + -cmf(i,1,3)*cmh_k(i)
1713 dglas(i,10)=dglas(i,10)+cmf(i,2,2)*cyy_k(i)-cmf(i,1,2)*cxx_k(i)
1714 + +cmf(i,2,3)*cmh_k(i)
1715 dglas(i,13)=dglas(i,13)+cmf(i,3,3)*(-cfh(i))+
1716 + cmf(i,1,3)*hxx(i)*vhg(i,3)-cmf(i,2,3)*hyy(i)*vhg(i,4)
1717 dglas(i,14)=dglas(i,14)-cmf(i,3,3)*cmh(i)+
1718 + cmf(i,1,3)*hxx(i)*vhg(i,1)-cmf(i,2,3)*hyy(i)*vhg(i,2)
1719 dglas(i,15)=dglas(i,15)+cmf(i,3,3)*cfh_k(i)
1720 + -cmf(i,1,3)*hxx_k(i)*vhg(i,3)+cmf(i,2,3)*hyy_k(i)*vhg(i,4)
1721 dglas(i,16)=dglas(i,16)+cmf(i,3,3)*cmh_k(i)
1722 + -cmf(i,1,3)*hxx_k(i)*vhg(i,1)+cmf(i,2,3)*hyy_k(i)*vhg(i,2)
1723 END DO
1724 END IF !(IORTH>1) THEN
1725 ELSE
1726 DO i=jft,jlt
1727 c1m(i)=a11(i)*fac1(i)*dt1
1728 c2m(i)=a12(i)*fac1(i)*dt1
1729 END DO
1730 DO i=jft,jlt
1731 dglas(i,1) =c1m(i)*cxx(i)-c2m(i)*cyy(i)+
1732 + c1m(i)*vhgze(i,1)+c2m(i)*vhgze(i,2)
1733 dglas(i,2) =c1m(i)*cyy(i)-c2m(i)*cxx(i)
1734 + -c1m(i)*vhgze(i,2)-c2m(i)*vhgze(i,1)
1735 dglas(i,3) =c1m(i)*bxx(i)-c2m(i)*byy(i)
1736 dglas(i,4) =c1m(i)*byy(i)-c2m(i)*bxx(i)
1737 dglas(i,7) =c1m(i)*cxx_k(i)-c2m(i)*cyy_k(i)
1738 + -c1m(i)*vhgzk(i,1)-c2m(i)*vhgzk(i,2)
1739 dglas(i,8) =c1m(i)*cyy_k(i)-c2m(i)*cxx_k(i)+
1740 + c1m(i)*vhgzk(i,2)+c2m(i)*vhgzk(i,1)
1741 dglas(i,9) =c1m(i)*bxx_k(i)-c2m(i)*byy_k(i)
1742 dglas(i,10)=c1m(i)*byy_k(i)-c2m(i)*bxx_k(i)
1743C------------
1744 c2=fac1(i)*gs(i)*one_over_64*dt1
1745 dglas(i,5) =c2*hxx(i)*vhg(i,5)
1746 dglas(i,6) =c2*hyy(i)*vhg(i,5)
1747 dglas(i,11)=c2*hxx_k(i)*vhg(i,6)
1748 dglas(i,12)=c2*hyy_k(i)*vhg(i,6)
1749C
1750 c3m=g(i)*fac1(i)*dt1
1751 crz=half*krz(i)*fac1(i)*dt1
1752C --------add shear due to rz dof--no sign change like 2,4,7,9-------
1753 dglas(i,13)=c3m*(hxx(i)*vhg(i,2)-hyy(i)*vhg(i,1)+vhgze(i,3))
1754 dglas(i,14)=c3m*(hxx(i)*vhg(i,4)-hyy(i)*vhg(i,3))
1755 dglas(i,15)=c3m*(-hxx_k(i)*vhg(i,2)+hyy_k(i)*vhg(i,1)+
1756 + vhgzk(i,3))
1757 dglas(i,16)=c3m*(-hxx_k(i)*vhg(i,4)+hyy_k(i)*vhg(i,3))
1758C --------add srz dof--------
1759 dglas(i,17)=crz*vhgze(i,5)
1760 dglas(i,18)=crz*vhgzk(i,5)
1761 END DO
1762 END IF !(IORTH>0) THEN
1763 DO i=jft,jlt
1764 DO j=1,18
1765 vglas(i,j) =vglas(i,j) +dglas(i,j)
1766 END DO
1767 END DO
1768C-----------elasto-plastic----
1769 CALL czhoureprz_or(jft ,jlt ,sigy ,vstre ,mstre ,
1770 1 thk ,fac ,esx0 ,emx ,npt ,
1771 2 dglas ,vglas ,nel )
1772C -------for energy calculation-Mean_value-VGLAS(t+dt)-------
1773C-----------to use the same coef: C5
1774 DO i=jft,jlt
1775 ss1(i)= my34(i)*vglas(i,1) +my23(i)*vglas(i,7)
1776 ss2(i)= mx23(i)*vglas(i,8) +mx34(i)*vglas(i,2)
1777 sf1(i)= my34(i)*vglas(i,3) +my23(i)*vglas(i,9)
1778 sf2(i)= -(mx23(i)*vglas(i,10)+mx34(i)*vglas(i,4))
1779 sc5(i)=my34(i)*vglas(i,5) +mx34(i)*vglas(i,6)
1780 sc6(i)=my23(i)*vglas(i,11)+mx23(i)*vglas(i,12)
1781 END DO !I=JFT,JLT
1782C
1783 DO i=jft,jlt
1784 ss1(i)= ss1(i) -mx34(i)*vglas(i,13)+mx23(i)*vglas(i,15)+
1785 + half*(mx34(i)*vglas(i,17)-mx23(i)*vglas(i,18))
1786 ss2(i)= ss2(i) +my34(i)*vglas(i,13) -my23(i)*vglas(i,15)+
1787 + half*(my34(i)*vglas(i,17)-my23(i)*vglas(i,18))
1788 sf1(i)=sf1(i)-mx34(i)*vglas(i,14)+mx23(i)*vglas(i,16)
1789 sf2(i)=sf2(i)-my34(i)*vglas(i,14)+my23(i)*vglas(i,16)
1790 erz = vglas(i,1)*vhgze(i,1)-vglas(i,2)*vhgze(i,2)-
1791 . vglas(i,7)*vhgzk(i,1)+vglas(i,8)*vhgzk(i,2)+
1792 . vglas(i,13)*vhgze(i,3)+vglas(i,15)*vhgzk(i,3)+
1793 . half*(vglas(i,17)*vhgze(i,4)+vglas(i,18)*vhgzk(i,4))
1794C--------
1795 esx(i)=ss1(i)*vhg(i,1)+ss2(i)*vhg(i,2)+erz*area(i)*fourth
1796 + + fourth*(sc5(i)*vhg(i,5)+sc6(i)*vhg(i,6))
1797 emx(i)=(sf1(i)*vhg(i,3)-sf2(i)*vhg(i,4))*c6(i)
1798 END DO
1799 DO i=jft,jlt
1800 eint(i,1) = eint(i,1)+c5(i)*esx(i)
1801 eint(i,2) = eint(i,2)+c5(i)*emx(i)
1802 END DO
1803 DO i=jft,jlt
1804 c8(i) = c7*off(i)
1805 ss1(i)= ss1(i)*c8(i)
1806 ss2(i)= ss2(i)*c8(i)
1807 sf1(i)= sf1(i)*c8(i)
1808 sf2(i)= sf2(i)*c8(i)
1809 c2 = c8(i)*thk(i)
1810 sc5(i)=sc5(i)*c2
1811 sc6(i)=sc6(i)*c2
1812 ss3(i)=sc5(i)+sc6(i)
1813 END DO !I=JFT,JLT
1814C-----------non const mz----17,18--
1815 DO i=jft,jlt
1816 c2=thk(i)*one_over_6*off(i)
1817 c3=thk(i)*third*off(i)
1818 DO j=1,4
1819 vmz(i,j)= vmz(i,j)
1820 . +c2*(vglas(i,17)*bmerz(i,4,j)+vglas(i,18)*bmkrz(i,4,j))
1821 . +c3*(bmerz(i,1,j)*vglas(i,1)-bmkrz(i,1,j)*vglas(i,7)
1822 . -bmerz(i,2,j)*vglas(i,2)+bmkrz(i,2,j)*vglas(i,8)
1823 . +bmerz(i,3,j)*vglas(i,13)+bmkrz(i,3,j)*vglas(i,15))
1824 END DO
1825 ENDDO
1826C--------Fint :terms for linear damping
1827 DO i=jft,jlt
1828C-------save for Fz--and add coupling terms(warped)------
1829 ssz1(i) = ss1(i)
1830 ssz2(i) = ss2(i)
1831 hvl = amu(i)*sqrt(rho(i)*area(i)*fac1(i))*off(i)
1832 ssv0=my23(i)*my23(i)
1833 ssv1=my34(i)*my34(i)
1834 ssv2=mx23(i)*mx23(i)
1835 ssv3=mx34(i)*mx34(i)
1836 hxx_v= stier*(ssv1+ssv0)
1837 hxy_v=-stier*(my34(i)*mx34(i)+my23(i)*mx23(i))
1838 hyy_v= stier*(ssv2+ssv3)
1839 c2 =hvl*gsr(i)*shfsr(i)*undouzsr
1840 cxz_v=(ssv1+ssv3)*c2
1841 cyz_v=(ssv2+ssv0)*c2
1842C
1843 aux=aa(i)*hvl
1844 c1m(i) = a11sr(i)*aux
1845 c2m(i) = a12sr(i)*aux
1846C
1847 cxx_v=c1m(i)*hxx_v
1848 cyy_v=c1m(i)*hyy_v
1849 cxy_v=c2m(i)*hxy_v
1850 ss1_v=cxx_v*vhg(i,1)+cxy_v*vhg(i,2)
1851 ss2_v=cyy_v*vhg(i,2)+cxy_v*vhg(i,1)
1852 sf1_v= (cxx_v*vhg(i,3)+cxy_v*vhg(i,4))*fbend_v
1853 sf2_v=(-cyy_v*vhg(i,4)-cxy_v*vhg(i,3))*fbend_v
1854 sc5_v=cxz_v*vhg(i,5)*hsura(i)
1855 sc6_v=cyz_v*vhg(i,6)*hsura(i)
1856 ss3_v = sc5_v+sc6_v
1857C--------------------------------------------------
1858 ss1(i)=ss1(i)+ss1_v
1859 ss2(i)=ss2(i)+ss2_v
1860 ss3(i)=ss3(i)+ss3_v
1861 sc5(i)=sc5(i)+sc5_v
1862 sc6(i)=sc6(i)+sc6_v
1863 sf1(i)=sf1(i)+sf1_v
1864 sf2(i)=sf2(i)+sf2_v
1865 tesy(i)= ((ss1_v*vhg(i,1)+ss2_v*vhg(i,2))*thk(i)+
1866 . (sf1_v*vhg(i,3)-sf2_v*vhg(i,4))*thk(i)*c6(i) +
1867 . sc5_v*vhg(i,5)+sc6_v*vhg(i,6))*dt1
1868 END DO !I=JFT,JLT
1869C--------Fint
1870 DO i=jft,jlt
1871 y13s= my13(i)*ss3(i)
1872 x13s= mx13(i)*ss3(i)
1873 y34s6=my34(i)*sc6(i)
1874 y23s5=my23(i)*sc5(i)
1875 x23s5=mx23(i)*sc5(i)
1876 x34s6=mx34(i)*sc6(i)
1877 c2=fourth*thk(i)
1878C--------NOEUDS 1,3-------
1879 b13=(my13(i)*x24(i)-mx13(i)*y24(i))*hsura(i)
1880 vf(i,1,1)=vf(i,1,1)+b13*ss1(i)
1881 vf(i,1,3)=c2*ss1(i)
1882 vf(i,2,1)=vf(i,2,1)+b13*ss2(i)
1883 vf(i,2,3)=c2*ss2(i)
1884 vf(i,3,3)=ss3(i)
1885C--------NOEUDS 2,4-------
1886 b24=(mx13(i)*y13(i)-my13(i)*x13(i))*hsura(i)
1887 vf(i,1,2)= vf(i,1,2)+b24*ss1(i)
1888 vf(i,1,4)=-vf(i,1,3)
1889 vf(i,2,2)= vf(i,2,2)+b24*ss2(i)
1890 vf(i,2,4)=-vf(i,2,3)
1891 vf(i,3,4)=-vf(i,3,3)
1892C--------NOEUDS 1,3-------
1893 c3=c6(i)*b13
1894 c4=c6(i)*c2
1895 vm(i,1,1)=vm(i,1,1)+c3*sf2(i)+y23s5+y34s6
1896 vm(i,1,3)=vm(i,1,3)+c4*sf2(i)-y13s
1897 vm(i,2,1)=vm(i,2,1)+c3*sf1(i)-x23s5-x34s6
1898 vm(i,2,3)=vm(i,2,3)+c4*sf1(i)+x13s
1899C--------NOEUDS 2,4-------
1900 c3=c6(i)*b24
1901 vm(i,1,2)=vm(i,1,2)+c3*sf2(i)+y23s5-y34s6
1902 vm(i,1,4)=vm(i,1,4)-c4*sf2(i)-y13s
1903 vm(i,2,2)=vm(i,2,2)+c3*sf1(i)-x23s5+x34s6
1904 vm(i,2,4)=vm(i,2,4)-c4*sf1(i)+x13s
1905C
1906 c2=z1(i)*hsura(i)
1907 vf(i,3,1)=vf(i,3,1)+c2*(ssz1(i)*y24(i)-ssz2(i)*x24(i))
1908 vf(i,3,2)=vf(i,3,2)+c2*(-ssz1(i)*y13(i)+ssz2(i)*x13(i))
1909 END DO !I=JFT,JLT
1910C
1911 CALL czehourv(jft ,jlt ,ipartc,evis,kfts ,tesy )
1912 RETURN
1913 END
1914!||====================================================================
1915!|| czhourep_or ../engine/source/elements/shell/coquez/czfintn.F
1916!||--- called by ------------------------------------------------------
1917!|| czfintn_or ../engine/source/elements/shell/coquez/czfintn.F
1918!||====================================================================
1919 SUBROUTINE czhourep_or(JFT ,JLT ,SIGY ,VSTRE ,MSTRE ,
1920 1 THK ,FAC ,ESX ,EMX ,NPT ,
1921 2 DGLAS ,VGLAS,NEL )
1922C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1923#include "implicit_f.inc"
1924#include "mvsiz_p.inc"
1925C-----------------------------------------------
1926C D U M M Y A R G U M E N T S
1927C-----------------------------------------------
1928 INTEGER JFT ,JLT,NPT,NEL
1929 my_real
1930 . thk(*) ,vglas(nel,12),vstre(nel,5),mstre(nel,3),
1931 . fac(mvsiz,2),sigy(*),esx(*),emx(*),dglas(mvsiz,12)
1932C-----------------------------------------------
1933C L O C A L V A R I A B L E S
1934C-----------------------------------------------
1935 INTEGER I,MX, IC, II, J,NPT1
1936 my_real
1937 . svm,svmn,c1b,c2b,cmm,cnn,ufac,tol,coef,eh1,eh2,
1938 . cnnx,cnny,cnnx_k,cnny_k,cmmx,cmmy,cmmx_k,cmmy_k,sxy0,mxy0,
1939 . svmc,sigy2,svm0,svmc_k,a_bk,a_b,a_k,undouzsr,aux,
1940 . coef1,coefh,sigy0,sxyh,mxyh,sigyh,facm,facf,fact
1941C-----------------------------------------------
1942 coefh= zep99
1943 coef= zep85
1944 tol= em18
1945 IF (npt==0) THEN
1946 coef1=sixteen
1947 ELSE
1948 coef1=twenty5
1949 ENDIF
1950C---- FAC(1,I)= mean value of Et/G; FAC(2,I)= min(Et_ip/G ), negative value when failure
1951C---- en Elastic FAC(1,I)=0;FAC(2,I)= 1
1952 DO i=jft,jlt
1953 IF (esx(i)<zero.AND.emx(i)<zero) cycle
1954 IF (sigy(i)<zep9ep30) THEN
1955 ufac=abs(min(fac(i,1),abs(fac(i,2)))-one)
1956C---------elastic inc on quadrature-----------
1957 IF (ufac<tol) THEN
1958 sigy2 = sigy(i)*sigy(i)
1959 sxy0= vstre(i,1)*vstre(i,1)+vstre(i,2)*vstre(i,2)
1960 . -vstre(i,1)*vstre(i,2)+three*vstre(i,3)*vstre(i,3)
1961 mxy0= mstre(i,1)*mstre(i,1)+mstre(i,2)*mstre(i,2)
1962 . -mstre(i,1)*mstre(i,2)+three*mstre(i,3)*mstre(i,3)
1963 sigy0 = sxy0+coef1*mxy0
1964C------------due to different criteria, SIGY0: Von misese
1965 sigy2 = max(sigy2,sigy0)
1966 cnn=coef
1967 cmm=coef*thk(i)*one_over_16
1968 cnnx=cnn*vglas(i,1)
1969 cnny=cnn*vglas(i,2)
1970 cnnx_k=cnn*vglas(i,7)
1971 cnny_k=cnn*vglas(i,8)
1972 cmmx=cmm*vglas(i,3)
1973 cmmy=cmm*vglas(i,4)
1974 cmmx_k=cmm*vglas(i,9)
1975 cmmy_k=cmm*vglas(i,10)
1976 sxyh= cnnx*cnnx+cnny*cnny-cnnx*cnny
1977 mxyh= cmmx*cmmx+cmmy*cmmy-cmmx*cmmy
1978 sxyh= sxyh+cnnx_k*cnnx_k+cnny_k*cnny_k-cnnx_k*cnny_k
1979 mxyh= mxyh+cmmx_k*cmmx_k+cmmy_k*cmmy_k-cmmx_k*cmmy_k
1980 sxyh= sxyh+abs(cnnx*(two*cnnx_k-cnny_k)
1981 1 +cnny*(two*cnny_k-cnnx_k))
1982 mxyh= mxyh+abs(cmmx*(two*cmmx_k-cmmy_k)
1983 1 +cmmy*(two*cmmy_k-cmmx_k))
1984 sigyh = sxyh+coef1*mxyh
1985C
1986 svm=sigy0+sigyh
1987 eh1=zero
1988 eh2=zero
1989C
1990 IF (svm > sigy2) THEN
1991 fact = svm/max(sigy2,tol) -one
1992 eh1 = min(fact,coefh)
1993 eh2 = eh1
1994 ENDIF
1995 ELSE
1996 eh1=one-fac(i,1)
1997 eh2=one-abs(fac(i,2))
1998 END IF
1999 IF (esx(i)<zero) eh1=zero
2000 IF (emx(i)<zero) eh2=zero
2001C-------------MEMBRANE------------------------
2002 vglas(i,1)=vglas(i,1)-eh1*dglas(i,1)
2003 vglas(i,2)=vglas(i,2)-eh1*dglas(i,2)
2004 vglas(i,7)=vglas(i,7)-eh1*dglas(i,7)
2005 vglas(i,8)=vglas(i,8)-eh1*dglas(i,8)
2006C-------------FLEXION------------------------
2007 vglas(i,3) =vglas(i,3) -eh2*dglas(i,3)
2008 vglas(i,4) =vglas(i,4) -eh2*dglas(i,4)
2009 vglas(i,9) =vglas(i,9) -eh2*dglas(i,9)
2010 vglas(i,10)=vglas(i,10)-eh2*dglas(i,10)
2011 ENDIF
2012 END DO !I=JFT,JLT
2013C-----taking into account failing in the layers(law25 first)----
2014 npt1=min(8,npt)
2015 facm = zep85 +zep015*(npt1-2)
2016 facf = zep7 +zep015*(npt1-2)
2017 fact = zep9999
2018 DO i=jft,jlt
2019 IF (fac(i,2)<zero) THEN
2020 vglas(i,1)=facm*vglas(i,1)
2021 vglas(i,2)=facm*vglas(i,2)
2022 vglas(i,7)=facm*vglas(i,7)
2023 vglas(i,8)=facm*vglas(i,8)
2024 vglas(i,5)=fact*vglas(i,5)
2025 vglas(i,6)=fact*vglas(i,6)
2026 vglas(i,11)=fact*vglas(i,11)
2027 vglas(i,12)=fact*vglas(i,12)
2028 vglas(i,3)=facf*vglas(i,3)
2029 vglas(i,4)=facf*vglas(i,4)
2030 vglas(i,9)=facf*vglas(i,9)
2031 vglas(i,10)=facf*vglas(i,10)
2032 ENDIF
2033 END DO
2034C
2035 RETURN
2036 END
2037!||====================================================================
2038!|| czhoureprz_or ../engine/source/elements/shell/coquez/czfintn.f
2039!||--- called by ------------------------------------------------------
2040!|| czfintnrz_or ../engine/source/elements/shell/coquez/czfintn.F
2041!||====================================================================
2042 SUBROUTINE czhoureprz_or(JFT ,JLT ,SIGY ,VSTRE ,MSTRE ,
2043 1 THK ,FAC ,ESX ,EMX ,NPT ,
2044 2 DGLAS ,VGLAS ,NEL )
2045C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2046#include "implicit_f.inc"
2047#include "mvsiz_p.inc"
2048C-----------------------------------------------
2049C D U M M Y A R G U M E N T S
2050C-----------------------------------------------
2051 INTEGER JFT ,JLT,NPT,NEL
2052 my_real
2053 . thk(*) ,vglas(nel,19),vstre(nel,5),mstre(nel,3),
2054 . fac(mvsiz,2),sigy(*),esx(*),emx(*),dglas(mvsiz,18)
2055C-----------------------------------------------
2056C L O C A L V A R I A B L E S
2057C-----------------------------------------------
2058 INTEGER I,MX, IC, II, J,NPT1
2059 my_real
2060 . svm,svmn,c1b,c2b,cmm,cnn,ufac,tol,coef,eh1,eh2,
2061 . cnnx,cnny,cnnx_k,cnny_k,cmmx,cmmy,cmmx_k,cmmy_k,sxy0,mxy0,
2062 . svmc,sigy2,svm0,svmc_k,a_bk,a_b,a_k,undouzsr,aux,
2063 . coef1,coefh,cnnxy,cnnxy_k,cmmxy,cmmxy_k,facm,facf,fact,sigy0
2064C-----------------------------------------------
2065 coefh= zep99
2066 coef= zep85
2067 tol= em18
2068 IF (npt==0) THEN
2069 coef1=sixteen
2070 ELSE
2071 coef1=twenty5
2072 ENDIF
2073 DO i=jft,jlt
2074 IF (esx(i)<zero.AND.emx(i)<zero) cycle
2075 IF (sigy(i)<zep9ep30) THEN
2076 ufac=abs(min(fac(i,1),abs(fac(i,2)))-one)
2077 sigy2 = sigy(i)*sigy(i)
2078 svm =zero
2079 sxy0=zero
2080 IF (ufac<tol) THEN
2081c CYCLE
2082 sxy0= vstre(i,1)*vstre(i,1)+vstre(i,2)*vstre(i,2)
2083 . -vstre(i,1)*vstre(i,2)+three*vstre(i,3)*vstre(i,3)
2084 mxy0= mstre(i,1)*mstre(i,1)+mstre(i,2)*mstre(i,2)
2085 . -mstre(i,1)*mstre(i,2)+three*mstre(i,3)*mstre(i,3)
2086 sigy0 = sxy0+coef1*mxy0
2087 sigy2 = max(sigy2,sigy0)
2088 cnn=coef
2089 cmm=coef*thk(i)*one_over_16
2090 cnnx=cnn*vglas(i,1)
2091 cnny=cnn*vglas(i,2)
2092 cnnxy=cnn*vglas(i,13)
2093 cnnx_k=cnn*vglas(i,7)
2094 cnny_k=cnn*vglas(i,8)
2095 cnnxy_k=cnn*vglas(i,15)
2096 cmmx=cmm*vglas(i,3)
2097 cmmy=cmm*vglas(i,4)
2098 cmmxy=cmm*vglas(i,14)
2099 cmmx_k=cmm*vglas(i,9)
2100 cmmy_k=cmm*vglas(i,10)
2101 cmmxy_k=cmm*vglas(i,16)
2102C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2103 sxy0= sxy0+cnnx*cnnx+cnny*cnny-cnnx*cnny+three*cnnxy*cnnxy
2104 mxy0= mxy0+cmmx*cmmx+cmmy*cmmy-cmmx*cmmy+three*cmmxy*cmmxy
2105 sxy0= sxy0+cnnx_k*cnnx_k+cnny_k*cnny_k-cnnx_k*cnny_k+
2106 1 three*cnnxy_k*cnnxy_k
2107 mxy0= mxy0+cmmx_k*cmmx_k+cmmy_k*cmmy_k-cmmx_k*cmmy_k+
2108 1 three*cmmxy_k*cmmxy_k
2109 sxy0= sxy0+abs(cnnx*(two*cnnx_k-cnny_k)
2110 1 +cnny*(two*cnny_k-cnnx_k))
2111 mxy0= mxy0+abs(cmmx*(two*cmmx_k-cmmy_k)
2112 1 +cmmy*(two*cmmy_k-cmmx_k))
2113 svm = sxy0+coef1*mxy0
2114 eh1 = zero
2115 eh2 = zero
2116 IF (svm > sigy2) THEN
2117 fact = svm/max(sigy2,tol) -one
2118 eh1 = min(fact,coefh)
2119 eh2 = eh1
2120 ENDIF
2121 ELSE
2122 eh1=one-fac(i,1)
2123 eh2=one-abs(fac(i,2))
2124 END IF
2125 IF (esx(i)<zero) eh1=zero
2126 IF (emx(i)<zero) eh2=zero
2127C-------------MEMBRANE------------------------
2128 vglas(i,1)=vglas(i,1)-eh1*dglas(i,1)
2129 vglas(i,2)=vglas(i,2)-eh1*dglas(i,2)
2130 vglas(i,7)=vglas(i,7)-eh1*dglas(i,7)
2131 vglas(i,8)=vglas(i,8)-eh1*dglas(i,8)
2132 vglas(i,13)=vglas(i,13)-eh1*dglas(i,13)
2133 vglas(i,15)=vglas(i,15)-eh1*dglas(i,15)
2134C-------------FLEXION------------------------
2135 vglas(i,3) =vglas(i,3) -eh2*dglas(i,3)
2136 vglas(i,4) =vglas(i,4) -eh2*dglas(i,4)
2137 vglas(i,9) =vglas(i,9) -eh2*dglas(i,9)
2138 vglas(i,10)=vglas(i,10)-eh2*dglas(i,10)
2139 vglas(i,14)=vglas(i,14)-eh1*dglas(i,14)
2140 vglas(i,16)=vglas(i,16)-eh1*dglas(i,16)
2141 ENDIF
2142 END DO !I=JFT,JLT
2143C-----taking into account failing in the layers(law25 first)----
2144 npt1=min(8,npt)
2145 facm = zep85 +zep015*(npt1-2)
2146 facf = zep7 +zep015*(npt1-2)
2147 fact = zep9999
2148 DO i=jft,jlt
2149 IF (fac(i,2)<zero) THEN
2150 vglas(i,1)=facm*vglas(i,1)
2151 vglas(i,2)=facm*vglas(i,2)
2152 vglas(i,7)=facm*vglas(i,7)
2153 vglas(i,8)=facm*vglas(i,8)
2154 vglas(i,13)=facm*vglas(i,13)
2155 vglas(i,15)=facm*vglas(i,15)
2156 vglas(i,5)=fact*vglas(i,5)
2157 vglas(i,6)=fact*vglas(i,6)
2158 vglas(i,11)=fact*vglas(i,11)
2159 vglas(i,12)=fact*vglas(i,12)
2160 vglas(i,3)=facf*vglas(i,3)
2161 vglas(i,4)=facf*vglas(i,4)
2162 vglas(i,9)=facf*vglas(i,9)
2163 vglas(i,10)=facf*vglas(i,10)
2164 vglas(i,14)=facf*vglas(i,14)
2165 vglas(i,16)=facf*vglas(i,16)
2166 ENDIF
2167 END DO
2168C
2169 RETURN
2170 END
2171!||====================================================================
2172!|| czfintnm1 ../engine/source/elements/shell/coquez/czfintn.F
2173!||--- called by ------------------------------------------------------
2174!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
2175!||--- calls -----------------------------------------------------
2176!|| czfintnm ../engine/source/elements/shell/coquez/czfintn.F
2177!||====================================================================
2178 SUBROUTINE czfintnm1(JFT ,JLT ,THK ,AA ,VHG ,
2179 2 X13 ,X24 ,Y13 ,Y24 ,VF ,
2180 3 MX13 ,MX23 ,MY13 ,MY23 ,
2181 4 G ,RHO ,AREA ,AMU ,DT1 ,
2182 5 V13 ,V24 ,VHI ,VGLAS,
2183 6 OFF ,IPARTC,EVIS,KFTS ,NEL )
2184C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2185C CALCUL viscous VF(3,*) due au shear when npt=1
2186C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2187#include "implicit_f.inc"
2188#include "mvsiz_p.inc"
2189#include "param_c.inc"
2190C-----------------------------------------------
2191C D U M M Y A R G U M E N T S
2192C-----------------------------------------------
2193 INTEGER JFT,JLT,IPARTC(*),NEL
2194 my_real
2195 . thk(*) ,aa(*) ,vhg(mvsiz,6),
2196 . x13(*) ,x24(*) ,y13(*) ,y24(*) ,
2197 . mx13(*),mx23(*),my13(*),my23(*),
2198 . vf(mvsiz,3,4),g(*),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
2199 . rho(*) ,area(*), off(*),amu(*) ,
2200 . evis(npsav,*),dt1,vglas(nel,12)
2201 INTEGER KFTS
2202C-----------------------------------------------
2203C L O C A L V A R I A B L E S
2204C-----------------------------------------------
2205 INTEGER I,MX, IC, II, J, JST(MVSIZ+1)
2206 my_real
2207 . c2,hsura,mx34,my34,fv,gv(mvsiz),
2208 . sc6_v,sc5_v,cxz_v,cyz_v,ss3_v,hvl,hvz(mvsiz),
2209 . gama1(mvsiz), gama2(mvsiz), gama3(mvsiz), gama4(mvsiz),
2210 . gamavz(mvsiz),px1(mvsiz), px2(mvsiz), py1(mvsiz), py2(mvsiz),
2211 . px1v, px2v, py1v, py2v, tesy(mvsiz),hx(mvsiz),hy(mvsiz),hxa,hya,fac
2212C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2213c gama_i=1/4(h_i-hx*bx_i-hy*byi); bx_i=(y24 -y13 -y24 y13)/A;by_i=(-x24 x13 x24 -x13)/A;
2214 DO i=jft,jlt
2215 px1(i) = y24(i)
2216 px2(i) = -y13(i)
2217 py1(i) = -x24(i)
2218 py2(i) = x13(i)
2219C----1/4*hx
2220 hx(i) = half*(-x13(i)+x24(i))+(mx13(i)-mx23(i))
2221 hy(i) = half*(-y13(i)+y24(i))+(my13(i)-my23(i))
2222 ENDDO
2223c gama_i*vz_i = vhi-hx*(y24*v13-y13*v24)-hy*(-x24*v13+x13*v24)
2224 DO i=jft,jlt
2225 gamavz(i) = vhi(i,3)-hx(i)*(px1(i)*v13(i,3)+px2(i)*v24(i,3))
2226 . -hy(i)*(py1(i)*v13(i,3)+py2(i)*v24(i,3))
2227 ENDDO
2228 DO i=jft,jlt
2229 hxa =hx(i)*aa(i)
2230 hya =hy(i)*aa(i)
2231 px1v = px1(i)*hxa
2232 px2v = px2(i)*hxa
2233 py1v = py1(i)*hya
2234 py2v = py2(i)*hya
2235C J=1,2 STOCK LA PARTIE ANTISYM (F(1)=-F(3),F(2)=-F(4))
2236C J=3,4 STOCK LA PARTIE SYM (F(1)=F(3),F(2)=F(4))
2237 gama1(i)= - px1v-py1v
2238 gama3(i)= fourth
2239 gama2(i)= - px2v-py2v
2240 gama4(i)= -fourth
2241 ENDDO
2242C----fixing fac as Q4, can be changed w/ G(i)
2243 fac = 0.83*em02
2244C--- same as Q4, no more damping but elastic w/ FAC*G
2245 DO i=jft,jlt
2246 hvl=fac*g(i)*thk(i)*off(i)*dt1
2247 vglas(i,10)=vglas(i,10)+hvl*gamavz(i)
2248 c2 = four*vglas(i,10)
2249 vf(i,3,1)=vf(i,3,1)+c2*gama1(i)
2250 vf(i,3,2)=vf(i,3,2)+c2*gama2(i)
2251 vf(i,3,3)=vf(i,3,3)+c2*gama3(i)
2252 vf(i,3,4)=vf(i,3,4)+c2*gama4(i)
2253C
2254 tesy(i) = area(i)*vglas(i,10)*dt1*(
2255 . (gama1(i)+gama2(i))*v13(i,3)+(gama3(i)+gama4(i))*v24(i,3))
2256 ENDDO
2257 mx = ipartc(jft)
2258 DO i=jft,jlt
2259 evis(8,mx)=evis(8,mx) + tesy(i)
2260 ENDDO
2261C---- damping
2262c FAC=0.01
2263cc FAC=0.5
2264c DO I=JFT,JLT
2265c HVZ(I)=FOURTH*FAC*RHO(I)*THK(I)*SQRT(AREA(I))
2266c ENDDO
2267cC---- mode Z1+Z2-(Z3+Z4)
2268c DO I=JFT,JLT
2269c FV=HVZ(I)*(V13(I,3)+V24(I,3))
2270c VF(I,3,1)=VF(I,3,1)+FV
2271c VF(I,3,2)=VF(I,3,2)+FV
2272c ENDDO
2273cC---- mode Z1+Z4-(Z2+Z3)
2274c DO I=JFT,JLT
2275c FV=HVZ(I)*(V13(I,3)-V24(I,3))
2276c VF(I,3,1)=VF(I,3,1)+FV
2277c VF(I,3,2)=VF(I,3,2)-FV
2278c ENDDO
2279
2280 gv(jft:jlt) = 0.001*g(jft:jlt)
2281 CALL czfintnm(jft ,jlt ,thk ,aa ,vhg ,
2282 2 x13 ,x24 ,y13 ,y24 ,vf ,
2283 3 gv ,rho ,area ,amu ,dt1 ,
2284 4 off ,ipartc,evis,kfts )
2285C
2286 RETURN
2287 END
my_real
#define my_real
Definition cppsort.cpp:32
czfintn1
subroutine czfintn1(jft, jlt, thk, c1, aa, vhg, x13, x24, y13, y24, z1, mx23, mx13, mx34, my13, my23, my34, vglas, vstre, mstre, vf, vm, fac, a11, a12, g, shf, sigy, off, fac1, rho, area, dt1, eint, amu, vhi, npt, ipartc, evis, kfts, gsr, nel, a11sr, a12sr, nusr, shfsr, mtn, fac58)
Definition czfintn.F:37
czhourep_or
subroutine czhourep_or(jft, jlt, sigy, vstre, mstre, thk, fac, esx, emx, npt, dglas, vglas, nel)
Definition czfintn.F:1922
czehourv
subroutine czehourv(jft, jlt, ipartc, evis, kfts, tesy)
Definition czfintn.F:1384
czfintn_or
subroutine czfintn_or(jft, jlt, thk, c1, aa, vhg, x13, x24, y13, y24, z1, mx23, mx13, mx34, my13, my23, my34, vglas, vstre, mstre, vf, vm, fac, a11, a12, g, gs, sigy, off, fac1, rho, area, dt1, eint, amu, vhi, npt, ipartc, evis, kfts, gsr, a11sr, a12sr, nusr, shfsr, iorth, hm, hf, hc, hmfor, mtn, nel)
Definition czfintn.F:1023
czfintnm1
subroutine czfintnm1(jft, jlt, thk, aa, vhg, x13, x24, y13, y24, vf, mx13, mx23, my13, my23, g, rho, area, amu, dt1, v13, v24, vhi, vglas, off, ipartc, evis, kfts, nel)
Definition czfintn.F:2184
czhoureprz_or
subroutine czhoureprz_or(jft, jlt, sigy, vstre, mstre, thk, fac, esx, emx, npt, dglas, vglas, nel)
Definition czfintn.F:2045
czfintnm
subroutine czfintnm(jft, jlt, thk, aa, vhg, x13, x24, y13, y24, vf, g, rho, area, amu, dt1, off, ipartc, evis, kfts)
Definition czfintn.F:506
czfintnrz
subroutine czfintnrz(jft, jlt, thk, c1, aa, vhg, x13, x24, y13, y24, z1, mx23, mx13, mx34, my13, my23, my34, vglas, vstre, mstre, vf, vm, fac, a11, a12, g, shf, sigy, off, fac1, rho, area, dt1, eint, amu, vhi, npt, ipartc, evis, kfts, gsr, a11sr, a12sr, nusr, shfsr, bmkrz, bmerz, vhgzk, vhgze, krz, vmz, nel)
Definition czfintn.F:573
czfintnrz_or
subroutine czfintnrz_or(jft, jlt, thk, c1, aa, vhg, x13, x24, y13, y24, z1, mx23, mx13, mx34, my13, my23, my34, vglas, vstre, mstre, vf, vm, fac, a11, a12, g, gs, sigy, off, fac1, rho, area, dt1, eint, amu, vhi, npt, ipartc, evis, kfts, gsr, a11sr, a12sr, nusr, shfsr, bmkrz, bmerz, vhgzk, vhgze, krz, vmz, iorth, hm, hf, hc, hmfor, mtn, nel)
Definition czfintn.F:1462
if
if(complex_arithmetic) id
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21