×

Notice

The forum is in read only mode.

Script to calculate a deformation with a rotating force vector.

  • Patufet
  • Topic Author
  • Offline
  • New Member
  • New Member
More
15 years 3 months ago #3967 by Patufet
Hello,

I have been doing some calculations to analyse the deformation of a mechanical part with Salome-Meca and Code Aster. I am interested in the deformation of several points of the structure depending on the angle of the applied force. For instance I run several identical scripts by just changing the variable ANGLE:

[code:1]ANGLE = 30 ; # angle in degrees
AMPLITUDE=5 ; # amplitude in N
PI=3.141592654;
COMP_X = cos(ANGLE*PI/180)*AMPLITUDE ;
COMP_Y = sin(ANGLE*PI/180)*AMPLITUDE ;

...

CHAR=AFFE_CHAR_MECA(MODELE=MODE,
DDL_IMPO=_F(GROUP_MA='Hold',
DX=0,
DY=0,
DZ=0,),
FORCE_FACE=_F(GROUP_MA='Press',
FX=COMP_X,FY=COMP_Y,),);[/code:1]


This works well, but it is not practical to measure many points. Is there a way to automate that process by introducing some kind of iteration on the script?

Thank you for your help.
More
15 years 3 months ago #3968 by kwou
Hoi

you can use a Python script. Eg. something along the line:

[code:1]
DEBUT();

import Numeric
# initialise counters
maxcount = 8;
counter = range(0,maxcount); # --> [0,1,2....7]

# initialise parameters
CHAR=[None]*maxcount
result=[None]*maxcount

AMPLITUDE=5 ; # amplitude in N
PI=Numeric.pi;


for count in counter:
ANGLE = count*360.0/maxcount; # make sure one number is float ...
COMP_X = cos(ANGLE*PI/180.0)*AMPLITUDE ;
COMP_Y = sin(ANGLE*PI/180.0)*AMPLITUDE ;

#....

CHAR[count]=AFFE_CHAR_MECA(MODELE=MODE,
DDL_IMPO=_F(GROUP_MA='Hold',
DX=0,
DY=0,
DZ=0,),
FORCE_FACE=_F(GROUP_MA='Press',
FX=COMP_X,FY=COMP_Y,),);
#.....
result[counter]=MECA_STATIQUE(MODELE=model,
CHAM_MATER=Amat,
EXCIT=(_F(CHARGE=CHAR[count],),),);


result[count] = CALC_ELEM(....);

IMPR_RESU(....,RESULTATE=result[count],....)
#....
# end of for loop (end of indent)

FIN()
[/code:1]

Please let me know if this works.
You can also have a look at the example zzzz121a.comm in the /opt/aster<version>/STA<version>/astest directory.
I used the Python loop for the Homard mesh refinement here: www.caelinux.org/wiki/index.php/Contrib:...s/Homard/lshape/loop

Although this can probably done a lot easier by using a list for the time and use a multiplification function for the load in MECA_STATIQUE or equivalent command. I am sure others will add that.

kind regards - kees<br /><br />Post edited by: Kees Wouters, at: 2010/03/21 00:41

Interest: structural mechanics, solar energy (picture at 'my location' shows too little pv panels)

--
kind regards - kees
  • Patufet
  • Topic Author
  • Offline
  • New Member
  • New Member
More
15 years 3 months ago #3969 by Patufet
Hello Kees,

thank you very much for your reply.
I have seen as you suggest some examples with a list for the time and a multiplication function for the load, but I did not know how to apply it to a rotating vector.
Anyway your script is exactly what I was looking for!

Thanks again.

Patufet.
More
15 years 3 months ago #3984 by kwou
Hoi Patufet,

I tried the other way also. I need to improve my french so I dived into the manuels and translated to and fro french/english/dutch.
I must say that I like the way to use VALE_PARA={list} and VALE_FONC={list} in DEFI_FONCTION very much. Much easier with Python to generate it.

[code:1]
DEBUT(PAR_LOT='NON');
import Numeric
....
pie = Numeric.pi;
Cload = pie*0.85*2;
Floadx=AFFE_CHAR_MECA(MODELE=modmod,
FORCE_ARETE=(_F(GROUP_MA=('Lload',),FX=2500.0/Cload,),),);
Floadz=AFFE_CHAR_MECA(MODELE=modmod,
FORCE_ARETE=(_F(GROUP_MA=('Lload',),FZ=2500.0/Cload,),),);


# define time function from [0,...1&gt;
# devided into tsteps
tsteps = 8; ## number of intervals
tend = float(tsteps-1)/(tsteps)

print 'tsteps: ',tsteps
print 'tend: ',tend
print 'pie: ',pie

time=DEFI_LIST_REEL(DEBUT=0.0,
INTERVALLE=_F(JUSQU_A=tend,NOMBRE=tsteps-1,),
INFO=2,TITRE='time',);

lenTime = len(time.Valeurs());

fcs = DEFI_LIST_REEL(INFO=2,TITRE='cos ',VALE = [(Numeric.cos(2.0*pie*ii/lenTime)) for ii in range(0,lenTime)], )
fsn = DEFI_LIST_REEL(INFO=2,TITRE='sin ',VALE = [(Numeric.sin(2.0*pie*ii/lenTime)) for ii in range(0,lenTime)], )
ftm = DEFI_LIST_REEL(INFO=2,TITRE='time',VALE = [(float(ii)/lenTime) for ii in range(0,lenTime)], )

# time and ftm are now effectively the same

rampx2 = DEFI_FONCTION(NOM_PARA='INST',
VALE_PARA=ftm,
VALE_FONC=fcs,
PROL_DROITE = 'CONSTANT',
PROL_GAUCHE = 'CONSTANT',
INFO=2, # nice to check the values
TITRE='rampx2',);

rampz2 = DEFI_FONCTION(NOM_PARA='INST',
VALE_PARA=ftm,
VALE_FONC=fsn,
PROL_DROITE = 'CONSTANT',
PROL_GAUCHE = 'CONSTANT',
INFO=2, # nice to check the values
TITRE='rampz2',);

resu=MECA_STATIQUE(MODELE=modmod,
CHAM_MATER=material,
CARA_ELEM=shellch,
LIST_INST=time,
INFO=2,
EXCIT=(_F(CHARGE=clamped,),
_F(CHARGE=Floadx,FONC_MULT=rampx2,),
_F(CHARGE=Floadz,FONC_MULT=rampz2,),),
OPTION='SIEF_ELGA_DEPL',);
....
[/code:1]

This way is probably much faster to solve since -I guess- the system matrix is only solved once for as many load steps as needed, defined in LIST_INST.

Kind regards - kees<br /><br />Post edited by: Kees Wouters, at: 2010/03/19 12:14

Interest: structural mechanics, solar energy (picture at 'my location' shows too little pv panels)

--
kind regards - kees
  • Patufet
  • Topic Author
  • Offline
  • New Member
  • New Member
More
15 years 3 months ago #3992 by Patufet
Hello Kees,

thank you very much for your comments.

I have tested the first approach that works well.

As:
[code:1] IMPR_RESU(FORMAT='MED',
UNITE=80,
RESU=_F(MAILLAGE=MAIL,
RESULTAT=RESU[count],
NOM_CHAM= ('SIGM_NOEU_DEPL','EQUI_NOEU_SIGM','DEPL',),),);[/code:1]

is in the loop, CA gives me a warning at each iteration after the first one:

-&gt; Le maillage MAIL est déjà présent dans le fichier med fort.80.


I am now testing your second approach.
For instance I just get an error at

[code:1]time=DEFI_LIST_REEL(DEBUT=0.0,
INTERVALLE=_F(JUSQU_A=tend,NOMBRE=tsteps-1,),
INFO=2,TITRE='time',);[/code:1]

tsteps: 8
tend: 0.875
JDC.py : ERREUR A L'INTERPRETATION DANS ACCAS - INTERRUPTION
&gt;&gt; JDC.py : DEBUT RAPPORT
CR phase d'initialisation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Erreur dans listr8.Valeurs en PAR_LOT='OUI' !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
fin CR phase d'initialisation


(I am using Salome-Meca 1.4.1)

Best Regards.
More
15 years 3 months ago #3997 by kwou
Hoi Patufet

[code:1]
IMPR_RESU(FORMAT='MED',
UNITE=80,
RESU=_F(MAILLAGE=MAIL,
RESULTAT=RESU[count],
NOM_CHAM= ('SIGM_NOEU_DEPL','EQUI_NOEU_SIGM','DEPL',),),);
[/code:1]
is in the loop, CA gives me a warning at each iteration after the first one:
-&gt; Le maillage MAIL est déjà présent dans le fichier med fort.80.


That warning is true of course. The IMPR_RESU command adds new results to the existing file MAIL/fort.80. You see it also in the post pro Salome window at result fields. I think it is nothing to worry about. The command is not re-entrant, so you cannot use IMPR_RESU(...reuse=MAIL....) or alike (as far as I can tell).

I am now testing your second approach.
For instance I just get an error at

[code:1]
time=DEFI_LIST_REEL(DEBUT=0.0,
INTERVALLE=_F(JUSQU_A=tend,NOMBRE=tsteps-1,),
INFO=2,TITRE='time',);

tsteps: 8
tend: 0.875

JDC.py : ERREUR A L'INTERPRETATION DANS ACCAS - INTERRUPTION
&gt;&gt; JDC.py : DEBUT RAPPORT
CR phase d'initialisation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Erreur dans listr8.Valeurs en PAR_LOT='OUI' !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
fin CR phase d'initialisation
[/code:1]


Oeps, sorry about that. I forgot that you need to change the start command:
DEBUT(PAR_LOT='NON');
Although I do not completely understand the French explanation, in case you use listr8 (and some of its 'python' methods) you need this. And I think Eficas is not 100% effective anymore when you use this kind of commands. I added it in the previous post.

(I am using Salome-Meca 1.4.1)

I use Salome 5.1.3 and Code Aster 10.1.x stand alone. A new version of SalomeMeca will be a available soon.

Kind regards - kees<br /><br />Post edited by: Kees Wouters, at: 2010/05/22 15:22

Interest: structural mechanics, solar energy (picture at 'my location' shows too little pv panels)

--
kind regards - kees
Moderators: catux
Time to create page: 0.169 seconds
Powered by Kunena Forum