Télécharger con_calc.procedur

Retour à la liste

Numérotation des lignes :

  1. * CON_CALC PROCEDUR BP208322 16/09/01 21:15:02 9010
  2. *=======================================================
  3. *
  4. * CON_CALC
  5. * ANALyse : calcul de la REPONSE_BALOURD NONlineaire par CONTINUATION
  6. * calcule le Résidu R, la raideur K=dR/dU, les dérivées dR/dt et dK/dt
  7. * creation : BP 28/01/2014
  8. *
  9. *=======================================================
  10. *
  11. ************************************************************************
  12. *
  13. * CREATION : BP208322 28/01/2014
  14. *
  15. * OBJET : Procedure de calcul de :
  16. * 1. du Residu = F^ext - F^int
  17. * 2. de la Raideur K = dResidu / du
  18. * 3. de dt * la dérivée du Residu par rapport a t
  19. * 4. de dt * la dérivée de la Raideur par rapport a t
  20. *
  21. * ENTREE :
  22. * TAB1 = TABLE
  23. * . 'MODELE'
  24. * . 'CARACTERISTIQUES'
  25. * . 'RIGIDITE_CONSTANTE'
  26. * . 'BLOCAGES_MECANIQUES'
  27. * . 'CHARGEMENT'
  28. * . 'TEMPS_CALCULES' = plage indicative des temps a calculer et
  29. * de la finesse des pas
  30. * . 'GRANDS_DEPLACEMENTS'
  31. * . 'K_SIGMA'
  32. * . 'MAXITERATION'
  33. * . 'PROCEDURE_CHAR_MECA' = vrai si chargement dependant de la confi
  34. * -guration (pressions suiveuses par ex.)
  35. * . 'WTABLE' = table avec infos de calcul interne
  36. * . 'CONTRAINTES' = contraintes de Cauchy a t_n+1
  37. * LMOT1 = LISTMOTS de mots-clés définissant les actions a effectuer
  38. * a choisir parmi { RESI(DU), RIGI, DRES(IDU), DRIG(I), ...}
  39. *
  40. * SORTIE :
  41. * TAB1 . 'WTABLE'
  42. * . 'RESIDU' = Residu = F^ext - F^int
  43. * . 'RAIDEUR' = K = dResidu/du
  44. * . 'DERIVEE_RESIDU' = dResidu/dt * dt
  45. * . 'DERIVEE_RAIDEUR' = dK/dt * dt
  46. *
  47. * REM :
  48. * on suppose etre sur la config finale n+1
  49. *
  50. ************************************************************************
  51.  
  52. DEBPROC CON_CALC TAB1*'TABLE' LMOT1*'LISTMOTS'
  53. time*'FLOTTANT' dtime*'FLOTTANT';
  54.  
  55.  
  56. ************************************************************************
  57. * *
  58. * VERIFICATION ou simple recup DES DONNEES D'ENTREE *
  59. * *
  60. ************************************************************************
  61. * fldebug = (vale debug) >eg 1;
  62. fldebug = faux;
  63.  
  64. SI (EXIS TAB1 'WTABLE');
  65. WTAB = TAB1 . 'WTABLE';
  66. SINON;
  67. MESS 'IL MANQUE LA WTABLE'; ERRE 21;
  68. FINSI;
  69.  
  70. FLCHAR1 = exis TAB1 'CHARGEMENT';
  71. si (FLCHAR1);
  72. CHAR1 = TAB1 . 'CHARGEMENT';
  73. finsi;
  74.  
  75. FLMOD1 = exis TAB1 'MODELE';
  76. si (FLMOD1);
  77. MOD1 = TAB1 . 'MODELE';
  78. MAT1 = TAB1 . 'CARACTERISTIQUES';
  79. SIG1 = WTAB . 'CONTRAINTES' ;
  80. finsi;
  81. DEP1 = WTAB . 'DEPLACEMENTS';
  82.  
  83. * prendre les matrices constantes ou HBM ?
  84. FLHBM = WTAB . 'HBM';
  85. * Alternating Frequency / Time ?
  86. IAFT = WTAB . 'AFT';
  87. * si AFT, Jacobienne calculee par TFR(Knl*gamma+Cnl*w*gamma°)
  88. * ou par differentiation numerique ?
  89. IDIFF = WTAB . 'DIFFERENTIATION';
  90. * dans le cas dynamique frequentiel, on a besoin de la frequence w
  91. * est-elle une inconnue ou fixee ?
  92. IFREQ = WTAB . 'FREQ_INCONNUE';
  93. si IFREQ; w = EXTR DEP1 'FREQ' WTAB . 'PT_FREQ';
  94. WTAB . 'FREQ' = w;
  95. sinon; w = IPOL TAB1 . 'FREQUENCE' time;
  96. finsi;
  97.  
  98. * initialisations
  99. TNL1 = mot 'NONDEFINI';
  100. Fext1 = mot 'NONDEFINI';
  101. TVARI1 = mot 'NONDEFINI';
  102.  
  103. SI FLHBM;
  104. morigi = mot 'RIGIDITE_HBM';
  105. momass = mot 'MASSE_HBM';
  106. moamor = mot 'AMORTISSEMENT_HBM';
  107. mocori = mot 'CORIOLIS_HBM';
  108. mokcen = mot 'CENTRIFUGE_HBM';
  109. mobloc = mot 'BLOCAGES_HBM';
  110. SINON;
  111. morigi = mot 'RIGIDITE_CONSTANTE';
  112. momass = mot 'MASSE_CONSTANTE';
  113. moamor = mot 'AMORTISSEMENT_CONSTANT';
  114. mobloc = mot 'BLOCAGES_MECANIQUES';
  115. FINSI;
  116.  
  117. KBLO1 = TAB1 . mobloc;
  118.  
  119.  
  120. ************************************************************************
  121. * *
  122. * CALCUL DU RESIDU *
  123. * *
  124. ************************************************************************
  125.  
  126. *** calcul des efforts Fext + Fnl ***
  127.  
  128. SI ((EXIS LMOT1 'RESI') ou (EXIS LMOT1 'DRES') ou IFREQ);
  129.  
  130. * Forces exterieure dues au chargement
  131. si (FLCHAR1); Fext1 = TIRE CHAR1 time;
  132. sinon; Fext1 = VIDE 'CHPOINT'/'DISCRET';
  133. finsi;
  134. * mess '_Fext1='; list resum Fext1;
  135.  
  136. * Autres forces (types suiveuses par ex ou F^NL perso)
  137. SI (TAB1 . 'PROCEDURE_CHAR_MECA');
  138. SI (EGA IAFT 0); TNL1 = CHARMECA TAB1 time; FINSI;
  139. SI (EGA IAFT 1); TNL1 = DFT TAB1 time; FINSI;
  140. SI (EGA IAFT 2); TNL1 = AFT TAB1 time LMOT1; FINSI;
  141. FNL1 = TNL1 . 'ADDI_SECOND';
  142. * mess 'concalc -> AFT 'time'=> FNL1='; list resum FNL1;
  143. Fext1 = Fext1 + FNL1;
  144. FINSI;
  145.  
  146. * Forces non-lineaires des structures variables BRASERO
  147. SI (EXIS TAB1 'STRUCT_VARIABLE');
  148. * attention a l'unite de w !
  149. TVARI1 = BRASE403 TAB1 DEP1 w LMOT1;
  150. FVARI1 = TVARI1 . 'ADDI_SECOND';
  151. Fext1 = Fext1 + FVARI1;
  152. FINSI;
  153.  
  154. FINSI;
  155.  
  156. *** calcul des effort internes ***
  157.  
  158. SI ((EXIS LMOT1 'RESI') ou (EXIS LMOT1 'DRES'));
  159.  
  160. * Forces internes
  161. * ! on suppose deja etre sur la config deformee !
  162. SI (FLMOD1); Fint1 = BSIG SIG1 MOD1 MAT1;
  163. SINON; Fint1 = VIDE 'CHPOINT'/'DISCRET';
  164. FINSI;
  165. * mess '_Fint1='; list resum Fint1;
  166.  
  167. * rappel : en frequentiel, K^dyn = K + t*C + t^2*M
  168. FLdyn1 = faux;
  169. SI (EXIS TAB1 morigi);
  170. Kdyn1 = TAB1 . morigi;
  171. * mess '_on prend K^Harm';
  172. FLdyn1 = vrai;
  173. SINON;
  174. Kdyn1 = VIDE 'RIGIDITE';
  175. FINSI;
  176. SI (EXIS TAB1 moamor);
  177. * attention a l'unite de w !
  178. Kdyn1 = Kdyn1 et (w * TAB1 . moamor);
  179. * mess '_on prend t * C^Harm';
  180. FLdyn1 = vrai;
  181. FINSI;
  182. SI (EXIS TAB1 momass);
  183. * attention a l'unite de w !
  184. Kdyn1 = Kdyn1 et ((w**2) * TAB1 . momass);
  185. * mess '_on prend t^2 * M^Harm';
  186. FLdyn1 = vrai;
  187. FINSI;
  188. SI (FLdyn1);
  189. * Fdyn1 = Kdyn1 * (DEP1 enle 'LX');
  190. Fdyn1 = Kdyn1 * DEP1;
  191. * mess '_Fdyn1='; list resum Fdyn1;
  192. Fint1 = Fint1 + Fdyn1;
  193. FINSI;
  194. *dl + prise en compte des reactions
  195. Freac1 = REAC KBLO1 DEP1;
  196. WTAB . 'REACTION' = Freac1;
  197. Fint1 = Fint1 - Freac1;
  198. *dl + :fin de l'ajout
  199.  
  200. * calcul et stockage effectif du Residu
  201. *bp-test SI (EXIS LMOT1 'RESI');
  202. Res1 = Fext1 - Fint1;
  203. * on enleve les deplacements imposés deja vérifiés
  204. Res1 = Res1 - ((KBLO1*DEP1) EXCO 'FLX' 'NOID' 'FLX');
  205. * mess '_Res1='; list resum Res1;
  206. WTAB . 'RESIDU' = Res1;
  207. *bp-test FINSI;
  208.  
  209. * on va calculer ici une norme pour tester la convergence
  210. * si (ega NF2 'NONDEFINI');
  211. * xden1 = maxi ((PSCA Res1 Res1 motCf motCf)**0.5);
  212. * xden2 = maxi ((PSCA Fext1 Fext1 motCf motCf)**0.5);
  213. * xden3 = maxi ((PSCA Fint1 Fint1 motCf motCf)**0.5);
  214. * mess 'xdeno = ' xden1 xden2 xden3;
  215. * * si (NF2 < 1.E-30); mess 'def de NF2 a revoir'; erre 21; finsi;
  216. * * WTAB . 'XDENO' = NF2;
  217. * finsi;
  218.  
  219. FINSI;
  220.  
  221.  
  222. ************************************************************************
  223. * *
  224. * CALCUL DE LA RAIDEUR (= - dResidu/dU) *
  225. * *
  226. ************************************************************************
  227. SI (EXIS LMOT1 'RIGI');
  228.  
  229. * Raideur elastique du modele
  230. * ! on suppose deja etre sur la config deformee !
  231. SI (FLMOD1); Ktot1 = RIGI MOD1 MAT1;
  232. SINON; Ktot1 = VIDE 'RIGIDITE';
  233. FINSI;
  234. * relations et blocages
  235. Ktot1 = Ktot1 et KBLO1;
  236.  
  237. * Ksigma
  238. SI (TAB1 . 'K_SIGMA');
  239. Ksig1 = KSIG SIG1 MOD1 MAT1;
  240. Ktot1 = Ktot1 et Ksig1;
  241. FINSI;
  242.  
  243. * Autres raideurs (dues au forces suiveuses ou autre K^NL perso)
  244. SI (TAB1 . 'PROCEDURE_CHAR_MECA');
  245. si (neg (type TNL1) 'TABLE');
  246. SI (EGA IAFT 0); TNL1 = CHARMECA TAB1 time; FINSI;
  247. SI (EGA IAFT 1); TNL1 = DFT TAB1 time; FINSI;
  248. SI (EGA IAFT 2); TNL1 = AFT TAB1 time LMOT1; FINSI;
  249. finsi;
  250. KNL1 = TNL1 . 'ADDI_MATRICE';
  251. Ktot1 = Ktot1 et KNL1;
  252. SINON;
  253. KNL1 = VIDE 'RIGIDITE';
  254. FINSI;
  255.  
  256. * raideurs non-lineaires des structures variables BRASERO
  257. SI (EXIS TAB1 'STRUCT_VARIABLE');
  258. si (neg (type TVARI1) 'TABLE');
  259. * attention a l'unite de w !
  260. TVARI1 = BRASE403 TAB1 DEP1 w LMOT1;
  261. finsi;
  262. KVARI1 = TVARI1 . 'ADDI_MATRICE';
  263. Ktot1 = Ktot1 et KVARI1;
  264. WTAB . 'RAIDEUR_NL' = KVARI1;
  265. FINSI;
  266.  
  267. * rappel : en frequentiel, K^dyn = K + w*C + w^2*M
  268. SI (EXIS TAB1 morigi);
  269. Ktot1 = Ktot1 et TAB1 . morigi;
  270. FINSI;
  271. SI (EXIS TAB1 moamor);
  272. * attention a l'unite de w !
  273. Ktot1 = Ktot1 et (w * TAB1 . moamor);
  274. FINSI;
  275. SI (EXIS TAB1 momass);
  276. * attention a l'unite de w !
  277. Ktot1 = Ktot1 et ((w**2) * TAB1 . momass);
  278. FINSI;
  279.  
  280. * cas ou w est une inconnue indpte : il faut ajouter -dR/dw a -dR/dU
  281. SI IFREQ;
  282. SI (TAB1 . 'PROCEDURE_CHAR_MECA');
  283. * on perturbe de 1.E-6
  284. dw1c = 1.E-6;
  285. w1c = w + dw1c;
  286. WTAB . 'FREQ' = w1c;
  287. * Fext est considéré indpt de w, mais pas FNL
  288. SI (EGA IAFT 0); TNL1c = CHARMECA TAB1 time; FINSI;
  289. * SI (EGA IAFT 1); TNL1c = DFT TAB1 time; FINSI;
  290. SI (EGA IAFT 2); TNL1c = AFT TAB1 time; FINSI;
  291. FNL1c = TNL1c . 'ADDI_SECOND';
  292. * -dR/dw = ...
  293. dRdw = (-1./dw1c) * (FNL1c - FNL1) ;
  294. dRdw = MANU 'RIGI' dRdw 'COLO' 'QUEL' WTAB . 'PT_FREQ' 'FREQ';
  295. Ktot1 = Ktot1 et dRdw;
  296. * on remet en l'etat initial
  297. WTAB . 'FREQ' = w;
  298. FINSI;
  299. * sans oublier la partie constante
  300. FLdyn1 = faux;
  301. SI (EXIS TAB1 moamor);
  302. dKdyn1 = TAB1 . moamor;
  303. FLdyn1 = vrai;
  304. SINON;
  305. dKdyn1 = VIDE 'RIGIDITE';
  306. FINSI;
  307. SI (EXIS TAB1 momass);
  308. dKdyn1 = dKdyn1 et ((2. * w) * TAB1 . momass);
  309. FLdyn1 = vrai;
  310. FINSI;
  311. SI FLdyn1;
  312. dZUdw = dKdyn1 * DEP1;
  313. dZUdw = MANU 'RIGI' dZUdw 'COLO' 'QUEL' WTAB . 'PT_FREQ' 'FREQ';
  314. Ktot1 = Ktot1 et dZUdw;
  315. FINSI;
  316. FINSI;
  317.  
  318. * stockage effectif de la raideur
  319. WTAB . 'RAIDEUR' = Ktot1;
  320.  
  321. FINSI;
  322.  
  323.  
  324. ************************************************************************
  325. * *
  326. * CALCUL DE LA DERIVEE DU RESIDU avec t *
  327. * *
  328. ************************************************************************
  329. SI (EXIS LMOT1 'DRES');
  330.  
  331. * REM : quand on ne sait pas faire autrement, on calcule F(t+dt1b)-F(t)
  332. * dt1b = 0.5*dtime;
  333. * dt1b = 0.01*dtime;
  334. si(dtime > time); dt1b = 1.E-2 * dtime;
  335. sinon ; dt1b = 1.E-6 * time;
  336. finsi;
  337. time1b = time + dt1b;
  338. si IFREQ;
  339. * cas ou w est une inconnue indpte: w ne varie pas explictement avec t
  340. w1b = w;
  341. sinon;
  342. w1b = IPOL TAB1 . 'FREQUENCE' time1b;
  343. si (neg dt1b 0.); dw1b = (w1b - w) / dt1b;
  344. sinon; dw1b = 0.;
  345. finsi;
  346. finsi;
  347.  
  348. * Variation des Forces exterieures dues au chargement
  349. si (FLCHAR1); Fext1b = TIRE CHAR1 time1b;
  350. sinon; Fext1b = VIDE 'CHPOINT'/'DISCRET';
  351. finsi;
  352.  
  353. * Variation des Autres forces (types suiveuses par ex ou F^NL perso)
  354. * dependent-elles explicitement de time?
  355. SI (TAB1 . 'PROCEDURE_CHAR_MECA');
  356. SI (EGA IAFT 0); TNL1b = CHARMECA TAB1 time1b; FINSI;
  357. * SI (EGA IAFT 1); TNL1b = DFT TAB1 time1b; FINSI;
  358. LMOT1b = MOTS 'DRES';
  359. SI (EGA IAFT 2); TNL1b = AFT TAB1 time1b LMOT1b; FINSI;
  360. FNL1b = TNL1b . 'ADDI_SECOND';
  361. * mess 'concalc -> AFT 'time1b'=> FNL1b='; list resum FNL1b;
  362. Fext1b = Fext1b + FNL1b;
  363. FINSI;
  364.  
  365. * Variation des forces non-lineaires des structures variables BRASERO
  366. * rem : calcul explicite pour l instant
  367. * on pourrait avoir directement les matrices ..? todo + tard
  368. SI (EXIS TAB1 'STRUCT_VARIABLE');
  369. * si w=0, avec structure variable amortissante, on n'a pas de
  370. * "raideur" -> system irresoluble si second membre non nul
  371. si (neg w 0.);
  372. TVARI1b = BRASE403 TAB1 DEP1 w1b (mots 'RESI');
  373. FVARI1b = TVARI1b . 'ADDI_SECOND';
  374. Fext1b = Fext1b + FVARI1b;
  375. finsi;
  376. FINSI;
  377.  
  378. * on calcule la variation des efforts en maintenant les u^imp a t
  379. * DFext1 = ((dtime / dt1b) * ((Fext1b - Fext1) ENLE 'FLX'))
  380. * + (EXCO Fext1 'FLX' 'NOID' 'FLX' 'NATURE' 'DISCRET');
  381. * on calcule la variation des efforts et des u^imp
  382. si (ega dtime 0.);
  383. DFext1 = 0.*(Fext1b - Fext1);
  384. sino;
  385. DFext1 = (dtime / dt1b) * (Fext1b - Fext1);
  386. finsi;
  387.  
  388. * Variation des Forces internes :
  389. * on les suppose ne dependre que de u et pas de t --> DFint1 = 0
  390. DFint1 = VIDE 'CHPOINT'/'DISCRET';
  391.  
  392. * rappel : en frequentiel, dK^dyn/dt = dw/dt * [C + 2*w*M]
  393. * et donc dF = dtime * [dK^dyn/dt]*(u enle LX)
  394. SI (NON IFREQ);
  395. FLdyn1 = faux;
  396. SI (EXIS TAB1 moamor);
  397. dKdyn1 = TAB1 . moamor;
  398. FLdyn1 = vrai;
  399. SINON;
  400. dKdyn1 = VIDE 'RIGIDITE';
  401. FINSI;
  402. SI (EXIS TAB1 momass);
  403. dKdyn1 = dKdyn1 et ((2. * w) * TAB1 . momass);
  404. FLdyn1 = vrai;
  405. FINSI;
  406. SI ((neg dtime 0.) et (neg dw1b 0.) et FLdyn1);
  407. DFdyn1 = dtime * dw1b * (dKdyn1 * DEP1);
  408. DFint1 = DFint1 + DFdyn1;
  409. FINSI;
  410. FINSI;
  411.  
  412.  
  413. * calcul et stockage effectif de dResidu/dt*dt
  414. DRes1 = DFext1 - DFint1;
  415. WTAB . 'DERIVEE_RESIDU' = DRes1;
  416.  
  417. FINSI;
  418.  
  419.  
  420.  
  421.  
  422. FINPROC ;
  423.  
  424.  
  425.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales