1 /************************************************************************
2 *
3 * ***** *** ***
4 * * * * * *
5 * * *** ***
6 * * * * * *
7 * ***** *** ***
8 *
9 * A FREE Finite Elements Analysis Program in ANSI C for the Windows
10 * and UNIX OS.
11 *
12 * Composed and edited and copyright by
13 * Professor Dr.-Ing. Frank Rieg, University of Bayreuth, Germany
14 *
15 * eMail:
16 * frank.rieg@uni-bayreuth.de
17 * dr.frank.rieg@t-online.de
18 *
19 * V15.0 November 18, 2015
20 *
21 * Z88 should compile and run under any Windows OS and UNIX OS and
22 * GTK+.
23 *
24 * This program is free software; you can redistribute it and/or modify
25 * it under the terms of the GNU General Public License as published by
26 * the Free Software Foundation; either version 2, or (at your option)
27 * any later version.
28 *
29 * This program is distributed in the hope that it will be useful,
30 * but WITHOUT ANY WARRANTY; without even the implied warranty of
31 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 * GNU General Public License for more details.
33 *
34 * You should have received a copy of the GNU General Public License
35 * along with this program; see the file COPYING. If not, write to
36 * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
37 ***********************************************************************/
38 /***********************************************************************
39 * Functions zur Vergleichs-Spannungsberechnung
40 *
41 * torgh: gh fuer Torus 6,8,12 und 15
42 * sheigh: gh fuer Scheibe 3,7,11 und 14
43 * hexgh: gh fuer Hexaeder 1,10,Tetraeder 16,17, Schalen 21,22,24
44 * platgh: gh fuer Platte 18,19 und 20
45 *
46 * tornh: nh fuer Torus 6,8,12 und 15
47 * sheinh: nh fuer Scheibe 3,7,11 und 14
48 * hexnh: nh fuer Hexaeder 1,10,Tetraeder 16,17, Schalen 21,22,24
49 * platnh: nh fuer Platte 18,19 und 20
50 *
51 * torsh: sh fuer Torus 6,8,12 und 15
52 * sheish: sh fuer Scheibe 3,7,11 und 14
53 * hexsh: sh fuer Hexaeder 1,10,Tetraeder 16,17, Schalen 21,22,24
54 * platsh: sh fuer Platte 18,19 und 20
55 *
56 * 2.1.2010 Rieg
57 ***********************************************************************/
58
59 /***********************************************************************
60 * Fuer UNIX
61 ***********************************************************************/
62 #ifdef FR_UNIX
63 #include <z88r.h>
64 #include <math.h>
65 #endif
66
67 /***********************************************************************
68 * Fuer Windows
69 ***********************************************************************/
70 #ifdef FR_WIN
71 #include <z88r.h>
72 #include <math.h>
73 #endif
74
75 /***********************************************************************
76 * Fuer Windows und GTK+
77 ***********************************************************************/
78 #ifdef FR_GTKWIN
79 #include <z88r.h>
80 #include <math.h>
81 #endif
82
83 /***********************************************************************
84 * Gestaltsaenderungsenergie-Hypothese, v.Mises stresses
85 ***********************************************************************/
86 /***********************************************************************
87 * hier beginnt Function torgh
88 ***********************************************************************/
torgh(FR_DOUBLE sig[])89 FR_DOUBLE torgh(FR_DOUBLE sig[])
90 {
91 FR_DOUBLE fret;
92 fret= FR_SQRT(sig[1]*sig[1] + sig[4]*sig[4] + sig[2]*sig[2] -
93 (sig[1]*sig[4] + sig[4]*sig[2] + sig[1]*sig[2]) +
94 3.*sig[3]*sig[3]);
95 return(fret);
96 }
97
98 /***********************************************************************
99 * hier beginnt Function sheigh
100 ***********************************************************************/
sheigh(FR_DOUBLE sig[])101 FR_DOUBLE sheigh(FR_DOUBLE sig[])
102 {
103 FR_DOUBLE fret;
104 fret= FR_SQRT(sig[1]*sig[1] + sig[2]*sig[2] - sig[1]*sig[2] +
105 3.*sig[3]*sig[3]);
106 return(fret);
107 }
108
109 /***********************************************************************
110 * hier beginnt Function hexgh
111 ***********************************************************************/
hexgh(FR_DOUBLE sig[])112 FR_DOUBLE hexgh(FR_DOUBLE sig[])
113 {
114 FR_DOUBLE fret;
115 fret= FR_SQRT(sig[1]*sig[1] + sig[2]*sig[2] + sig[3]*sig[3] -
116 (sig[1]*sig[2] + sig[2]*sig[3] + sig[1]*sig[3]) +
117 3.*(sig[4]*sig[4] + sig[5]*sig[5] + sig[6]*sig[6]));
118 return(fret);
119 }
120
121 /***********************************************************************
122 * hier beginnt Function platgh
123 ***********************************************************************/
platgh(FR_DOUBLE rsig[])124 FR_DOUBLE platgh(FR_DOUBLE rsig[])
125 {
126 FR_DOUBLE fret;
127 fret= FR_SQRT(rsig[1]*rsig[1] + rsig[2]*rsig[2]
128 -rsig[1]*rsig[2] + 3.*rsig[3]*rsig[3]);
129 return(fret);
130 }
131
132 /***********************************************************************
133 * Normalspannungs-Hypothese, principal stresses
134 ***********************************************************************/
135 /***********************************************************************
136 * hier beginnt Function tornh
137 ***********************************************************************/
tornh(FR_DOUBLE sig[])138 FR_DOUBLE tornh(FR_DOUBLE sig[])
139 {
140 FR_DOUBLE fret,sig1,sig2,sig3,wur;
141
142 /*----------------------------------------------------------------------
143 * Hauptspannungen
144 *---------------------------------------------------------------------*/
145 wur= FR_SQRT(sig[1]*sig[1]-2*sig[1]*sig[2]+sig[2]*sig[2]+4*sig[3]*sig[3]);
146
147 sig1= FR_FABS(0.5*(-wur+sig[1]+sig[2]));
148 sig2= FR_FABS(0.5*( wur+sig[1]+sig[2]));
149 sig3= FR_FABS(sig[4]);
150
151 /*----------------------------------------------------------------------
152 * Spannungshypothese
153 *---------------------------------------------------------------------*/
154 fret= 0.;
155 if(sig1 >= sig2 && sig1 >= sig3) return(sig1);
156 if(sig2 >= sig1 && sig2 >= sig3) return(sig2);
157 if(sig3 >= sig1 && sig3 >= sig2) return(sig3);
158 return(fret);
159 }
160
161 /***********************************************************************
162 * hier beginnt Function sheinh
163 ***********************************************************************/
sheinh(FR_DOUBLE sig[])164 FR_DOUBLE sheinh(FR_DOUBLE sig[])
165 {
166 FR_DOUBLE sig1,sig2,wur;
167 FR_DOUBLE fret;
168
169 /*----------------------------------------------------------------------
170 * Hauptspannungen
171 *---------------------------------------------------------------------*/
172 wur= FR_SQRT(FR_POW((sig[1]-sig[2]),2)/2. + sig[3]*sig[3]);
173
174 sig1= FR_FABS((sig[1]+sig[2])/2. + wur);
175 sig2= FR_FABS((sig[1]+sig[2])/2. - wur);
176
177 /*----------------------------------------------------------------------
178 * Spannungshypothese
179 *---------------------------------------------------------------------*/
180 fret= 0.;
181 if(sig1 >= sig2) return(sig1);
182 if(sig2 >= sig1) return(sig2);
183 return(fret);
184 }
185
186 /***********************************************************************
187 * hier beginnt Function hexnh
188 ***********************************************************************/
hexnh(FR_DOUBLE sig[])189 FR_DOUBLE hexnh(FR_DOUBLE sig[])
190 {
191 extern FR_INT4 noci;
192
193 FR_DOUBLE J1,J2,J3,D,p,q,phi,y1,y2,y3,wup,sig1,sig2,sig3;
194 FR_DOUBLE fret=0;
195
196 /*----------------------------------------------------------------------
197 * die Invarianten berechnen
198 *---------------------------------------------------------------------*/
199 J1= sig[1]+sig[2]+sig[3];
200 J2= sig[1]*sig[2]+sig[1]*sig[3]+sig[2]*sig[3]-
201 sig[4]*sig[4]-sig[6]*sig[6]-sig[5]*sig[5];
202 J3= sig[1]*sig[2]*sig[3]-sig[1]*sig[5]*sig[5]+
203 sig[4]*sig[5]*sig[6]-sig[3]*sig[4]*sig[4]+
204 sig[6]*sig[4]*sig[5]-sig[2]*sig[6]*sig[6];
205
206 /*----------------------------------------------------------------------
207 * reduzierte Form der kubischen Gleichung y^3+py+q=0 loesen
208 *---------------------------------------------------------------------*/
209 p= J2- J1*J1/3.;
210 q= J1*J2/3. - 2*J1*J1*J1/27. - J3;
211
212 D= q*q/4. + p*p*p/27.; /* Diskriminante */
213
214 if(D >= 0) /* casus irreducibilis geht nicht */
215 {
216 noci= 1;
217 return(0.);
218 }
219
220 phi= FR_ACOS(-q/2./FR_SQRT(FR_POW((FR_FABS(p)/3.),3)));
221
222 wup= FR_SQRT(FR_FABS(p)/3.);
223 y1= 2*wup*FR_COS(phi/3.);
224 y2= -2*wup*FR_COS(phi/3.-1.0472);
225 y3= -2*wup*FR_COS(phi/3.+1.0472);
226
227 /*----------------------------------------------------------------------
228 * Hauptspannungen
229 *---------------------------------------------------------------------*/
230 sig1= FR_FABS(y1+J1/3.);
231 sig2= FR_FABS(y2+J1/3.);
232 sig3= FR_FABS(y3+J1/3.);
233
234 /*----------------------------------------------------------------------
235 * Spannungshypothese
236 *---------------------------------------------------------------------*/
237 fret= 0.;
238 if(sig1 >= sig2 && sig1 >= sig3) return(sig1);
239 if(sig2 >= sig1 && sig2 >= sig3) return(sig2);
240 if(sig3 >= sig1 && sig3 >= sig2) return(sig3);
241 return(fret);
242 }
243
244 /***********************************************************************
245 * hier beginnt Function platnh
246 ***********************************************************************/
platnh(FR_DOUBLE rsig[])247 FR_DOUBLE platnh(FR_DOUBLE rsig[])
248 {
249 FR_DOUBLE sig1,sig2,wur;
250 FR_DOUBLE fret;
251
252 /*----------------------------------------------------------------------
253 * Hauptspannungen
254 *---------------------------------------------------------------------*/
255 wur= FR_SQRT(FR_POW((rsig[1]-rsig[2]),2)/2. + rsig[3]*rsig[3]);
256
257 sig1= FR_FABS((rsig[1]+rsig[2])/2. + wur);
258 sig2= FR_FABS((rsig[1]+rsig[2])/2. - wur);
259
260 /*----------------------------------------------------------------------
261 * Spannungshypothese
262 *---------------------------------------------------------------------*/
263 fret= 0.;
264 if(sig1 >= sig2) return(sig1);
265 if(sig2 >= sig1) return(sig2);
266 return(fret);
267 }
268
269 /***********************************************************************
270 * Schubspannungs-Hypothese, Tresca stresses
271 ***********************************************************************/
272 /***********************************************************************
273 * hier beginnt Function torsh
274 ***********************************************************************/
torsh(FR_DOUBLE sig[])275 FR_DOUBLE torsh(FR_DOUBLE sig[])
276 {
277 FR_DOUBLE fret,sig1,sig2,sig3,wur;
278
279 /*----------------------------------------------------------------------
280 * Hauptspannungen
281 *---------------------------------------------------------------------*/
282 wur= FR_SQRT(sig[1]*sig[1]-2*sig[1]*sig[2]+sig[2]*sig[2]+4*sig[3]*sig[3]);
283
284 sig1= 0.5*(-wur+sig[1]+sig[2]);
285 sig2= 0.5*( wur+sig[1]+sig[2]);
286 sig3= sig[4];
287
288 /*----------------------------------------------------------------------
289 * Spannungshypothese
290 *---------------------------------------------------------------------*/
291 fret= 0.;
292 if(sig1 >= sig2 && sig2 >= sig3) return(sig1-sig3);
293 if(sig1 >= sig3 && sig3 >= sig2) return(sig1-sig2);
294 if(sig2 >= sig1 && sig1 >= sig3) return(sig2-sig3);
295 if(sig2 >= sig3 && sig3 >= sig1) return(sig2-sig1);
296 if(sig3 >= sig2 && sig2 >= sig1) return(sig3-sig1);
297 if(sig3 >= sig1 && sig1 >= sig2) return(sig3-sig2);
298 return(fret);
299 }
300
301 /***********************************************************************
302 * hier beginnt Function sheish
303 ***********************************************************************/
sheish(FR_DOUBLE sig[])304 FR_DOUBLE sheish(FR_DOUBLE sig[])
305 {
306 FR_DOUBLE sig1,sig2,wur;
307 FR_DOUBLE fret;
308
309 /*----------------------------------------------------------------------
310 * Hauptspannungen
311 *---------------------------------------------------------------------*/
312 wur= FR_SQRT(FR_POW((sig[1]-sig[2]),2)/2. + sig[3]*sig[3]);
313
314 sig1= (sig[1]+sig[2])/2. + wur;
315 sig2= (sig[1]+sig[2])/2. - wur;
316
317 /*----------------------------------------------------------------------
318 * Spannungshypothese
319 *---------------------------------------------------------------------*/
320 fret= 0.;
321 if(sig1 >= sig2 && sig2 >= 0.) return(FR_FABS(sig1));
322 if(sig1 <= sig2 && sig2 <= 0.) return(FR_FABS(sig1));
323 if(sig2 >= sig1 && sig1 >= 0.) return(FR_FABS(sig2));
324 if(sig2 <= sig1 && sig1 <= 0.) return(FR_FABS(sig2));
325 if(sig1 >= 0. && sig2 <= 0.) return(sig1-sig2);
326 if(sig2 >= 0. && sig1 <= 0.) return(sig2-sig1);
327 return(fret);
328 }
329
330 /***********************************************************************
331 * hier beginnt Function hexsh
332 ***********************************************************************/
hexsh(FR_DOUBLE sig[])333 FR_DOUBLE hexsh(FR_DOUBLE sig[])
334 {
335 extern FR_INT4 noci;
336
337 FR_DOUBLE J1,J2,J3,D,p,q,phi,y1,y2,y3,wup,sig1,sig2,sig3;
338 FR_DOUBLE fret=0;
339
340 /*----------------------------------------------------------------------
341 * die Invarianten berechnen
342 *---------------------------------------------------------------------*/
343 J1= sig[1]+sig[2]+sig[3];
344 J2= sig[1]*sig[2]+sig[1]*sig[3]+sig[2]*sig[3]-
345 sig[4]*sig[4]-sig[6]*sig[6]-sig[5]*sig[5];
346 J3= sig[1]*sig[2]*sig[3]-sig[1]*sig[5]*sig[5]+
347 sig[4]*sig[5]*sig[6]-sig[3]*sig[4]*sig[4]+
348 sig[6]*sig[4]*sig[5]-sig[2]*sig[6]*sig[6];
349
350 /*----------------------------------------------------------------------
351 * reduzierte Form der kubischen Gleichung y^3+py+q=0 loesen
352 *---------------------------------------------------------------------*/
353 p= J2- J1*J1/3.;
354 q= J1*J2/3. - 2*J1*J1*J1/27. - J3;
355
356 D= q*q/4. + p*p*p/27.; /* Diskriminante */
357
358 if(D >= 0) /* casus irreducibilis geht nicht */
359 {
360 noci= 1;
361 return(0.);
362 }
363
364 phi= FR_ACOS(-q/2./FR_SQRT(FR_POW((FR_FABS(p)/3.),3)));
365
366 wup= FR_SQRT(FR_FABS(p)/3.);
367 y1= 2*wup*FR_COS(phi/3.);
368 y2= -2*wup*FR_COS(phi/3.-1.0472);
369 y3= -2*wup*FR_COS(phi/3.+1.0472);
370
371 /*----------------------------------------------------------------------
372 * Hauptspannungen
373 *---------------------------------------------------------------------*/
374 sig1= y1+J1/3.;
375 sig2= y2+J1/3.;
376 sig3= y3+J1/3.;
377
378 /*----------------------------------------------------------------------
379 * Spannungshypothese
380 *---------------------------------------------------------------------*/
381 fret= 0.;
382 if(sig1 >= sig2 && sig2 >= sig3) return(sig1-sig3);
383 if(sig1 >= sig3 && sig3 >= sig2) return(sig1-sig2);
384 if(sig2 >= sig1 && sig1 >= sig3) return(sig2-sig3);
385 if(sig2 >= sig3 && sig3 >= sig1) return(sig2-sig1);
386 if(sig3 >= sig2 && sig2 >= sig1) return(sig3-sig1);
387 if(sig3 >= sig1 && sig1 >= sig2) return(sig3-sig2);
388 return(fret);
389 }
390
391 /***********************************************************************
392 * hier beginnt Function platsh
393 ***********************************************************************/
platsh(FR_DOUBLE rsig[])394 FR_DOUBLE platsh(FR_DOUBLE rsig[])
395 {
396 FR_DOUBLE sig1,sig2,wur;
397 FR_DOUBLE fret;
398
399 /*----------------------------------------------------------------------
400 * Hauptspannungen
401 *---------------------------------------------------------------------*/
402 wur= FR_SQRT(FR_POW((rsig[1]-rsig[2]),2)/2. + rsig[3]*rsig[3]);
403
404 sig1= (rsig[1]+rsig[2])/2. + wur;
405 sig2= (rsig[1]+rsig[2])/2. - wur;
406
407 /*----------------------------------------------------------------------
408 * Spannungshypothese
409 *---------------------------------------------------------------------*/
410 fret= 0.;
411 if(sig1 >= sig2 && sig2 >= 0.) return(FR_FABS(sig1));
412 if(sig1 <= sig2 && sig2 <= 0.) return(FR_FABS(sig1));
413 if(sig2 >= sig1 && sig1 >= 0.) return(FR_FABS(sig2));
414 if(sig2 <= sig1 && sig1 <= 0.) return(FR_FABS(sig2));
415 if(sig1 >= 0. && sig2 <= 0.) return(sig1-sig2);
416 if(sig2 >= 0. && sig1 <= 0.) return(sig2-sig1);
417 return(fret);
418 }
419