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 * V14.0 January 14, 2009
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 * Diese Compilereinheit BETI88 enthaelt:
40 * - beti88 fuer Balken Nr.25 in allg.Lage, Bernoulli + Timoshenko
41 * - sbet88, die dazugehoerige Spannungsroutine
42 * 4.1.2013 Rieg
43 * 27.02.2013 Wehmann, Aenderung zs fuer Rechtssystem
44 * 01.03.2013 Wehmann, Spannungsberechnung Timoshenko hinzugefuegt
45 * 27.08.2013 Rieg, umbenennen von vormals BETI88.C, das beti88 und
46 * sbet88 enthielt, in timo88 und stim88, damit ist
47 * die Namensgebung wieder konsequent zum Rest, d.h.
48 * NAME88 = Steifigkeitsmatrix Element
49 * BNAME88= Lastvektoren fuer das Element
50 * SNAM88 = Spannungsroutine Element
51 ***********************************************************************/
52
53 /***********************************************************************
54 * Fuer UNIX
55 ***********************************************************************/
56 #ifdef FR_UNIX
57 #include <z88r.h>
58 #include <stdio.h>
59 #endif
60
61 /***********************************************************************
62 * Fuer Windows
63 ***********************************************************************/
64 #ifdef FR_WIN
65 #include <z88r.h>
66 #include <stdio.h>
67 #endif
68
69 /***********************************************************************
70 * Fuer Windows & GTK+
71 ***********************************************************************/
72 #ifdef FR_GTKWIN
73 #include <z88r.h>
74 #include <stdio.h>
75 #endif
76
77 /***********************************************************************
78 * Leseformate
79 ***********************************************************************/
80 #define NL "\n"
81
82 #ifdef FR_XDOUB
83 #define P13E "%+#13.5lE "
84 #define P11E "%+#11.3lE"
85 #define P11EB "%+#11.3lE "
86 #endif
87
88 #ifdef FR_XQUAD
89 #define P13E "%+#13.5LE "
90 #define P11E "%+#11.3LE"
91 #define P11EB "%+#11.3LE "
92 #endif
93
94
95 /***********************************************************************
96 * hier beginnt Function stim88
97 ***********************************************************************/
stim88(void)98 int stim88(void)
99 {
100 extern FILE *fo3;
101
102 extern FR_DOUBLEAY tm;
103
104 extern FR_DOUBLE ul[];
105 extern FR_DOUBLE xk[],yk[],zk[];
106
107 extern FR_DOUBLE emode,rnuee,qparae,ezze,eyye,rite,wte;
108 extern FR_DOUBLE riyye, rizze;
109 extern FR_DOUBLE xkp,ykp,zkp,rkape;
110
111 extern FR_INT4 ifbetie,ifnili;
112
113 FR_DOUBLE ulc[13];
114
115 FR_DOUBLE a21x,a21y,a21z,rl21,a31x,a31y,a31z,rl31;
116 FR_DOUBLE azsx,azsy,azsz,rlzs,aysx,aysy,aysz,rlys;
117 FR_DOUBLE vxx,vxy,vxz,vyx,vyy,vyz,vzx,vzy,vzz;
118
119 FR_DOUBLE rlv,fac,qrl,sum;
120 FR_DOUBLE sigxx,tauxx,sigzz1,sigzz2,sigyy1,sigyy2;
121
122 FR_DOUBLE gmode;
123
124 FR_DOUBLE BeTiBzz, BeTiRzz, BeTiSzz;
125 FR_DOUBLE BeTiByy, BeTiRyy, BeTiSyy;
126
127 FR_INT4 i,k,inc;
128
129 /*----------------------------------------------------------------------
130 * globale Verschiebungen in lokale Verschiebungen umrechnen,
131 * Transformationsmatrix aufstellen
132 *---------------------------------------------------------------------*/
133 /*======================================================================
134 * Einheits-Vektor 2-1 spannt die lokale x-Achse x' auf
135 *=====================================================================*/
136 a21x= xk[2] - xk[1];
137 a21y= yk[2] - yk[1];
138 a21z= zk[2] - zk[1];
139 rl21 = FR_SQRT(a21x*a21x + a21y*a21y + a21z*a21z);
140
141 a21x/= rl21;
142 a21y/= rl21;
143 a21z/= rl21;
144
145 /*======================================================================
146 * Einheits-Vektor 3-1 fuer Kreuzprodukt ays,xkp/ykp/zkp=Kontrollpunkt 3
147 *=====================================================================*/
148 a31x= xkp-xk[1];
149 a31y= ykp-yk[1];
150 a31z= zkp-zk[1];
151 rl31= FR_SQRT(a31x*a31x + a31y*a31y + a31z*a31z);
152
153 a31x/= rl31;
154 a31y/= rl31;
155 a31z/= rl31;
156
157 /*======================================================================
158 * ays ist Kreuzprodukt aus a21 und a31 und bildet die lokale Achse y'
159 *=====================================================================*/
160 aysx= a21y*a31z - a21z*a31y;
161 aysy= a21z*a31x - a21x*a31z;
162 aysz= a21x*a31y - a21y*a31x;
163 rlys = FR_SQRT(aysx*aysx + aysy*aysy + aysz*aysz);
164
165 aysx/= rlys;
166 aysy/= rlys;
167 aysz/= rlys;
168
169 /*======================================================================
170 * azs ist Kreuzprodukt aus ays und a21 und bildet die lokale Achse z'
171 *=====================================================================*/
172 azsx= - aysy*a21z + aysz*a21y;
173 azsy= - aysz*a21x + aysx*a21z;
174 azsz= - aysx*a21y + aysy*a21x;
175 rlzs = FR_SQRT(azsx*azsx + azsy*azsy + azsz*azsz); /* Kontrolle */
176
177 /*======================================================================
178 * Richtungscosinus x'x, x'y und x'z
179 *=====================================================================*/
180 vxx=a21x;
181 vxy=a21y;
182 vxz=a21z;
183
184 /*======================================================================
185 * Richtungscosinus y'x, y'y und y'z
186 *=====================================================================*/
187 vyx=aysx;
188 vyy=aysy;
189 vyz=aysz;
190
191 /*======================================================================
192 * Richtungscosinus z'x, z'y und z'z
193 *=====================================================================*/
194 vzx=azsx;
195 vzy=azsy;
196 vzz=azsz;
197
198 /*======================================================================
199 * Transformationsmatrix tm
200 *=====================================================================*/
201 for(i = 1; i <= 144; i++)
202 tm[i]= 0.;
203
204 inc= 0;
205
206 for(i = 1; i <= 4; i++)
207 {
208 tm[inc*12+inc +1]= vxx;
209 tm[inc*12+inc +2]= vxy;
210 tm[inc*12+inc +3]= vxz;
211 tm[inc*12+inc+12+1]= vyx;
212 tm[inc*12+inc+12+2]= vyy;
213 tm[inc*12+inc+12+3]= vyz;
214 tm[inc*12+inc+24+1]= vzx;
215 tm[inc*12+inc+24+2]= vzy;
216 tm[inc*12+inc+24+3]= vzz;
217 inc+= 3;
218 }
219
220 /*----------------------------------------------------------------------
221 * die globalen Verschiebungen zwischensichern
222 *---------------------------------------------------------------------*/
223 for(i = 1; i <= 12; i++)
224 ulc[i]= ul[i];
225
226 /*----------------------------------------------------------------------
227 * Matrix-Vektorprodukt: tm * ulc(global) = ul(lokal)
228 *---------------------------------------------------------------------*/
229 for(i = 1;i <= 12;i++)
230 {
231 sum= 0.;
232 for(k = 1; k <= 12; k++)
233 {
234 sum+= tm[(i-1)*12+k] * ulc[k];
235 }
236 ul[i]= sum; /* ul ist jetzt lokal */
237 }
238
239 /*----------------------------------------------------------------------
240 * Spannungen berechnen
241 *---------------------------------------------------------------------*/
242 if(ifbetie == 0)
243 {
244 rlv= rl21+ul[7]-ul[1];
245 sigxx= emode*(rlv/rl21-1.);
246
247 tauxx= (ul[10]-ul[4])/rl21 * emode/(2.*(1.+rnuee))*rite/wte;
248
249 fac= emode*ezze;
250 qrl= rl21*rl21;
251 sigzz1= fac*2.*(3.*(ul[8]-ul[2])-rl21*(2.*ul[6]+ul[12]))/qrl;
252 sigzz2= sigzz1 + fac*6.*(2.*(ul[2]-ul[8])+rl21*( ul[12]+ul[6]))/qrl;
253
254 fac= emode*eyye;
255 sigyy1= fac*2.*(3.*(ul[9]-ul[3])+rl21*(2.*ul[5]+ul[11]))/qrl;
256 sigyy2= sigyy1 + fac*6.*(2.*(ul[3]-ul[9])-rl21*( ul[11]+ul[5]))/qrl;
257
258 if(ifnili == 0)
259 fprintf(fo3,NL P11EB P11EB P11EB P11EB P11EB P11E,
260 sigxx,tauxx,sigzz1,sigyy1,sigzz2,sigyy2);
261 }
262 else
263 {
264 rlv= rl21+ul[7]-ul[1];
265 sigxx= emode*(rlv/rl21-1.);
266
267 tauxx= (ul[10]-ul[4])/rl21 * emode/(2.*(1.+rnuee))*rite/wte;
268
269 gmode= emode/(2.*(1.+rnuee));
270
271 fac= ezze/rizze;
272 BeTiBzz= (6.0*emode*rizze/(rl21*rl21))
273 / (1.0+12.0*emode*rizze/(rkape*gmode*qparae*rl21*rl21));
274 BeTiSzz= ( 4.0*emode*rizze/rl21
275 + 12.0*emode*emode*rizze*rizze/(rkape*gmode*qparae*rl21*rl21*rl21))
276 / (1.0+12.0*emode*rizze/(rkape*gmode*qparae*rl21*rl21));
277 BeTiRzz= ( 2.0*emode*rizze/rl21
278 - 12.0*emode*emode*rizze*rizze/(rkape*gmode*qparae*rl21*rl21*rl21))
279 / (1.0+12.0*emode*rizze/(rkape*gmode*qparae*rl21*rl21));
280
281 sigzz1= fac*(BeTiBzz*ul[2]+BeTiSzz*ul[6]-BeTiBzz*ul[8]+BeTiRzz*ul[12]);
282 sigzz2= fac*(BeTiBzz*ul[2]+BeTiRzz*ul[6]-BeTiBzz*ul[8]+BeTiSzz*ul[12]);
283
284 fac= eyye/riyye;
285 BeTiByy= (6.0*emode*riyye/(rl21*rl21))
286 / (1.0+12.0*emode*riyye/(rkape*gmode*qparae*rl21*rl21));
287 BeTiSyy= ( 4.0*emode*riyye/rl21
288 + 12.0*emode*emode*riyye*riyye/(rkape*gmode*qparae*rl21*rl21*rl21))
289 / (1.0+12.0*emode*riyye/(rkape*gmode*qparae*rl21*rl21));
290 BeTiRyy= ( 2.0*emode*riyye/rl21
291 - 12.0*emode*emode*riyye*riyye/(rkape*gmode*qparae*rl21*rl21*rl21))
292 / (1.0+12.0*emode*riyye/(rkape*gmode*qparae*rl21*rl21));
293
294 sigyy1= fac*(-BeTiByy*ul[3]+BeTiSyy*ul[5]+BeTiByy*ul[9]+BeTiRyy*ul[11]);
295 sigyy2= fac*(-BeTiByy*ul[3]+BeTiRyy*ul[5]+BeTiByy*ul[9]+BeTiSyy*ul[11]);
296
297 if(ifnili == 0)
298 fprintf(fo3,NL P11EB P11EB P11EB P11EB P11EB P11E,
299 sigxx,tauxx,sigzz1,sigyy1,sigzz2,sigyy2);
300 }
301
302
303 return(0);
304 }
305
306