Difference between revisions of "Contrib:KeesWouters/bc/pythonlist"
Keeswouters (Talk | contribs) m (→The cylindrical boundary conditions) |
Keeswouters (Talk | contribs) m (→The cylindrical boundary conditions) |
||
Line 51: | Line 51: | ||
LIAISON_DDL=(_F(NOEUD=('Ni','Ni','Ni'),DDL=('DX','DY','DZ'),COEF_MULT=(m2x,m2y,m2z),COEF_IMPO=0) | LIAISON_DDL=(_F(NOEUD=('Ni','Ni','Ni'),DDL=('DX','DY','DZ'),COEF_MULT=(m2x,m2y,m2z),COEF_IMPO=0) | ||
where 'Ni' is a node number located on the cylinder area, eg, 'N123'. This command has to be implemented for every node on the cylinder, so we use Python to make this a bit easier and faster. | where 'Ni' is a node number located on the cylinder area, eg, 'N123'. This command has to be implemented for every node on the cylinder, so we use Python to make this a bit easier and faster. | ||
− | |||
Now remember that the part in the _F(...) operator can be placed in a Python list, dictionary, or tuple, eg: | Now remember that the part in the _F(...) operator can be placed in a Python list, dictionary, or tuple, eg: | ||
− | |||
− | |||
− | |||
_F(NOEUD=('N808','N808','N808'),DDL=('DX','DY','DZ'),COEF_MULT=(0.0, -0.36, -0.26),COEF_IMPO=0.00) | _F(NOEUD=('N808','N808','N808'),DDL=('DX','DY','DZ'),COEF_MULT=(0.0, -0.36, -0.26),COEF_IMPO=0.00) | ||
− | |||
is the same as, or equivalent to: | is the same as, or equivalent to: | ||
− | |||
− | |||
− | |||
{'NOEUD': ['N808','N808','N808'],'COEF_MULT':[0.0,-0.36,-0.26], 'COEF_IMPO':0.0,'DDL':['DX','DY','DZ']} | {'NOEUD': ['N808','N808','N808'],'COEF_MULT':[0.0,-0.36,-0.26], 'COEF_IMPO':0.0,'DDL':['DX','DY','DZ']} | ||
− | + | So in the Python script we iterate over all nodes on the cylinder area and append | |
− | + | bc.append({'NOEUD': ['N808','N808','N808'],'COEF_MULT':[0.0,-0.36,-0.26], 'COEF_IMPO':0.0,'DDL':['DX','DY','DZ']}) | |
+ | this sequence for all selected nodes. | ||
===''Python code for list for free rotation''=== | ===''Python code for list for free rotation''=== |
Revision as of 15:08, 27 January 2013
Geometry and mesh of the block with cylindrical hole
Applied on a simple block with a cylindrical hole the use of LIAISON_DLL is shown to simulate a cylindrical coordinate system
Contrib:KeesWouters/bc/cylinder
Geometry and mesh of the block with cylindrical hole
- key words
o LIAISON_DDL o simulated cylinder coordinates and boundary conditions
Contrib:KeesWouters/bc/pythonlist
Rotation axis defined by Python list
- key words
o LIAISON_DDL o simulated cylinder coordinates and boundary conditions o Python list for applying boundary conditions
Contents
- 1 Geometry and mesh of the block with cylindrical hole
- 2 The cylindrical boundary conditions
- 3 The boundary conditions applied to the block
- 4 Updated comm file for Code-Aster v11
- 5 Input files for the FE Analysis
- 6 References
- 7 Material properties of the block
- 8 The boundary conditions
- 9 The results of the calculation
- 10 Input files for the FE Analysis and references
Geometry and mesh of the block with cylindrical hole
This construction shows the use of LIAISON_DDL to simulate boundary condition on a cylindrical hole. It just shows the use of it.
The boundary conditions ore build around a cylindrical hole in a block. The nodes attached to the cylinder are restricted to rotate around its central axis. The definition of this axis is done by two points P1 and P2 that need to be provided.
The geometry consists of a block with a cylindrical hole near the bottom side. The overall dimensions of the block are [Lx * Ly * Lz] = [2.0 * 3.5 * 20.0 ]. The hole is placed on the x-plane at position [yc, zc] = [2.0, 3.0]. The radius of the hole is R=0.45.
Some images of the construction:
A number of groups has been defined (P for plane, L for line segments):
- Ptop (not used)
- Pbot (not used)
- Pcyl (could be used to define nodes on this surface)
- Lcyl (not used in final version) and
- Ltop, used to apply boundary conditions in x and y direction.
The cylindrical boundary conditions
Below an image to illustrate how to restrict these degrees of freedom:
In vector notation or Numpy - Python notation :
m1 = p2-p1 ## rotation axis m2 = Node-p1 ## m2* in the picture m3 = cross(m2,m1) ## Node may not be on rotation axis--> m2=a*m1 -->m3=0??? m2 = cross(m1,m3) ## normal vector from line(P1,P2) to Node lm2 = sqrt(dot(m2,m2)) ## length of m2 m2 /=lm2 ## normalise m2 to unit vector
The two points P1 and P2 define the rotation axis (and, but not neccessarely, the cylinder). The nodes on the cylinder area may free move in the direction of the vector m1 and m3. In the direction of vector m2, normal to the rotation axis, the displacement is restricted. This restriction may be described by setting the inner product dot of the dispacement vector dX = [dx, dy, dz] and the normal vector m2 = [m2x, m2y, m2z] to 0:
<dX, m2> = <[dx, dy, dz], [m2x, m2y, m2z]> = 0.
This is exactly the restricted that can be described by the C-Aster command LIAISON_DDL:
LIAISON_DDL=(_F(NOEUD=('Ni','Ni','Ni'),DDL=('DX','DY','DZ'),COEF_MULT=(alpha1,alpha2,alpha3),COEF_IMPO=beta) or LIAISON_DDL=(_F(NOEUD=('Ni','Ni','Ni'),DDL=('DX','DY','DZ'),COEF_MULT=(m2x,m2y,m2z),COEF_IMPO=0)
where 'Ni' is a node number located on the cylinder area, eg, 'N123'. This command has to be implemented for every node on the cylinder, so we use Python to make this a bit easier and faster.
Now remember that the part in the _F(...) operator can be placed in a Python list, dictionary, or tuple, eg:
_F(NOEUD=('N808','N808','N808'),DDL=('DX','DY','DZ'),COEF_MULT=(0.0, -0.36, -0.26),COEF_IMPO=0.00)
is the same as, or equivalent to:
{'NOEUD': ['N808','N808','N808'],'COEF_MULT':[0.0,-0.36,-0.26], 'COEF_IMPO':0.0,'DDL':['DX','DY','DZ']}
So in the Python script we iterate over all nodes on the cylinder area and append
bc.append({'NOEUD': ['N808','N808','N808'],'COEF_MULT':[0.0,-0.36,-0.26], 'COEF_IMPO':0.0,'DDL':['DX','DY','DZ']})
this sequence for all selected nodes.
Python code for list for free rotation
Since we are going to mesh information in the Python script language, we need to define:
DEBUT(PAR_LOT='NON')
We will also use the some modules defined in the ...Utilitai.partition folder:
from Utilitai.partition import *
need_MAIL_PY_help=False ## obvious other choise is True if need_MAIL_PY_help: from Utilitai.partition import * help(MAIL_PY)
PyMesh = MAIL_PY() ## convert CodeAster mesh to Python PyMesh.FromAster(mesh); nonu = PyMesh.dime_maillage[0] ## number of nodes elnu = PyMesh.dime_maillage[2] ## number of elements test=CmeshCA.dime_maillage[5] ## [min,max] dimension: [0,5] NodeCoord = PyMesh.cn ## xyz coordinates of nodes (nonu*3 matrix) ElemConnect = PyMesh.co ## node numbers of elements (elnu*x matrix) NodeList = list(PyMesh.correspondance_noeuds) ElemList = list(PyMesh.correspondance_mailles) ElemType = list(PyMesh.tm) NodeGroup = PyMesh.gno ## list of GROUP_NO groups (see help(MAIL_PY) for object methods reference) ElemGroup = PyMesh.gma ## list of GROUP_MA groups, eg groupstr
This part of the script extract the nodes and coordinates of the cylinder area that need to have the rotational degrees of freedom. Note that this group of nodes has already been defined in Salome. The group is called 'Ncyl'
NodeCount = 0 ddl_condition=[] ## initialise ddl_condition list for ii in xrange(len(NodeGroup[groupstr])): NodeNumber = NodeGroup[groupstr][ii] #print 'NodeGroup: ',ii,NodeNumber,NodeCount NodeCount+=1 NNxyz = NodeCoord[NodeNumber] nx,ny,nz = cylinder_normal(P1,P2,NNxyz) NodeCount+=1
NNp1 = NodeNumber+1 # node numbers are increased by 1: Salome <--> Code Aster ddl_condition.append({'NOEUD': ['N%d'%(NNp1),'N%d'%(NNp1),'N%d'%(NNp1)], 'DDL':['DX','DY','DZ'], 'COEF_MULT': [nx,ny,nz],'COEF_IMPO':0.00}) if info>1: #print {'NOEUD': ['N%d'%(NNp1),'N%d'%(NNp1),'N%d'%(NNp1)], 'DDL':['DX','DY','DZ'], 'COEF_MULT': [0.0,ny-y1,nz-z1],'COEF_IMPO':0.00} print ddl_condition[ii]
if info>0: # check format of bc_list print ic,countntop,': ',ddl_condition[0] print '::::::::::::::::::::::::::' print ic,countntop,': ',ddl_condition[-1]
The AFFE_CHAR_MECA part can be compared to this (small part of the extraction) in the standard way:
ddl_condition=AFFE_CHAR_MECA(MODELE=Cmod, LIAISON_DDL = _F(NOEUD=('N162 ','N162 '),DDL=('DY','DZ'),COEF_MULT=( 0.312402110475 , -10.3238902922 ),COEF_IMPO=0.00), _F(NOEUD=('N12 ','N12 '),DDL=('DY','DZ'),COEF_MULT=( 0.0 , -9.55 ),COEF_IMPO=0.00), _F(NOEUD=('N695 ','N695 '),DDL=('DY','DZ'),COEF_MULT=(-0.443605677987 , -9.92440897899 ),COEF_IMPO=0.00), ... _F(NOEUD=('N765 ','N765 '),DDL=('DY','DZ'),COEF_MULT=(-0.162737538758 , -10.4195431962 ),COEF_IMPO=0.00), _F(NOEUD=('N551 ','N551 '),DDL=('DY','DZ'),COEF_MULT=( 0.107926195474 , -9.56313396066 ),COEF_IMPO=0.00),),);
Now the list can be used in the AFFE_CHAR_MECA command just as the _F(...) operator:
rotcyl=AFFE_CHAR_MECA(MODELE=Cmod,LIAISON_DDL = ddl_condition);
The boundary conditions applied to the block
We apply a displacement on a line of the top plane in y direction DY=0.5 and keep it fixed in x direction: DX=0:
- The boundary condition on the line segment Ltop:
- bcforce=AFFE_CHAR_MECA(MODELE=Cmod,DDL_IMPO=(_F(GROUP_MA='Ltop',DY=0.5,DX=0.0000),),);
Then we apply the cylindrical boundary conditions to the node on the cylinder and let is freely rotate around the cylinder axis:
- The cylindrical boundary condition on the nodes of the cylindrical hole are defined by the LIAISON_DDL keyword in stead of DDL_IMPO:
- ... LIAISON_DDL=(_F(NOEUD=('Ni','Ni','Ni'),DDL=('DX','DY','DZ'),COEF_MULT=(alpha1,alpha2,alpha3),COEF_IMPO=beta),
lright, it took me several days to find out just this. But once you know, it is easy.
Adding the boundary conditions to the calculation
Finally, the list is applied to the MECA_STATIQUE command.
result=MECA_STATIQUE(MODELE=Cmod, CHAM_MATER=Amat, EXCIT=(_F(CHARGE=bcforce,), _F(CHARGE=rotcyl),),);
Updated comm file for Code-Aster v11
updated january 2013
Code-Aster evolves, so a lot of commands have changed since v9-10 of Code Aster. The input files for v11 reflect these changes as well as some node selection.
Since we have to import a few Python files to define the rotation, in DEBUT we need to advice Code-Aster to talk to Python in an proper way:
DEBUT(PAR_LOT='NON');
import sys sys.path.append("/cae_sg500/caexample/caelinux/csys1/python2a") ## define your own path here for import RotationAxisBC from CA_geometry import RotationAxisBC import numpy
The selection of nodes on the cylinder is now included in the command file:
Cmesh=LIRE_MAILLAGE(UNITE=20,FORMAT='MED',);
Cmesh=DEFI_GROUP(reuse =Cmesh, MAILLAGE=Cmesh, CREA_GROUP_NO=_F(NOM='Ncyl',GROUP_MA='Pcyl',),);
The geometry group Pcyl [GROUP_MA='Pcyl'] need to be defined in the mesh file. This is easy in eg Salome. The selection of the node group GROUP_NO='Ncyl', previously done in Salome, is now carried out by the DEFI_GROUP command (which I find easier).
Now returning to the definition of the free rotation of a cylinder around its axis: two points P1 and P2 define the centre line of the rotation axis (which may coincide with the centre axis of a cylinder). Ideally, these two point would be based on the selection of nodes, but so far I have figured out how to do that relatively easy.
Anyway, the points P1 and P1 are numpy arrays and define the rotation axis. The Python function RotationAxisBC(...) has five arguments:
- mesh: Cmesh
- group of nodes GROUP_NO 'Ncyl' that are restriction to rotate around the given axis
- two points: numpy arrays, P1 and P2 defining the rotation axis and
- info: [0|1|2] for printing controle
P1 = numpy.array([0.0, 2.0, 3.0]) P2 = numpy.array([2.0, 2.0, 3.0]) info = 0 ddl_condition = RotationAxisBC(Cmesh,'Ncyl',P1,P2,info) rotcyl = AFFE_CHAR_MECA(MODELE=Cmod,LIAISON_DDL=ddl_condition);
The result is a Python list ddl_condition that contains the familiar parts of the LIAISON_DDL for all the nodes contained in the selection 'Ncyl'.
As always, the hard work is being done by the Python function RotationAxisBC(...). More on this later ..., see for now the download files.
Note that the CALC_ELEM procedure has been made obsolete.
Input files for the FE Analysis
Download the files here:
- >>>Media:kw_blockrot_pythonlist.zip for C-Aster version 10-
- >>>Media:kw_rotation_cylinder_ca11.zip for C-Aster version 11
Input files:
- Python script for defining the geometry and mesh (blockrot4.py), load by File --> load script (cntrl T in the object browser), refresh (F5) after running
- Save the mesh file by right clicking in the object browser by right clicking on the mesh name Mbcyl; select Export to MED and accept or change the default values (v11 saves mesh file automaticly in working directory)
- ASTK file (cylrot.astk, you need to edit the path to your requirements ...)
- command file (blockb.comm)
- blockrot6.py basicly construct the geometry and selects nodes on the cylinder. blockrot6.py directly selects the nodes on the cylinder area by filters. See also [Salome] for more reference (and thank you JMB for the question asked there).
References
CodeAster documentation:
[Document u4.44.01 - AFFE_CHAR_MECA, page 20/109]
[Document u1.03.02 - Methodes Python d'access aux objects Aster - french]
CodeAster forum:
[Using python with Aster]
[Applying nodal forces node by node]
[Accessing Mesh Data Structure]
==
=
==
Geometry and mesh of the block with cylindrical hole
Applied on a simple block with a cylindrical hole the use of LIAISON_DLL is shown to simulate a cylindrical coordinate system
Contrib:KeesWouters/bc/cylinder
Geometry and mesh of the block with cylindrical hole
- key words
o LIAISON_DDL o simulated cylinder coordinates and boundary conditions
Contrib:KeesWouters/bc/pythonlist
Rotation axis defined by Python list
- key words
o LIAISON_DDL o simulated cylinder coordinates and boundary conditions o Python list for applying boundary conditions
Material properties of the block
The material property of the block is set to steel.
The boundary conditions
- On the top line segment Ltop a non zero displacement in y direction is prescribed. The displacement is z direction is fixed.
- On the nodes connected to the cylindrical hole a tangential displacement is allowed. The radial component is fixed. The displacement in axial or x direction is free. Due to this restriction the displacement of the geometry in z direction is defined.
The boundary conditions in detail
- The first boundary condition on the line segment Ltop is easy to define:
- bcforce=AFFE_CHAR_MECA(MODELE=Cmod,DDL_IMPO=(_F(GROUP_MA='Ltop',DY=0.5,DX=0.0000),),);
- The cylindrical boundary condition on the nodes of the cylindrical hole are defined by the LIAISON_DDL keyword in stead of DDL_IMPO:
- ... LIAISON_DDL=(_F(NOEUD=('Ni','Ni'),DDL=('DY','DZ'),COEF_MULT=(alpha1,alpha2),COEF_IMPO=beta),
Using this keyword and selecting all the nodes for Ni can simulate a tangential boundary condition. Details follow now.
- ... LIAISON_DDL=(_F(NOEUD=('Ni','Ni'),DDL=('DY','DZ'),COEF_MULT=(alpha1,alpha2),COEF_IMPO=beta),
With LIAISON_DDL it is possible to couple displacements of various nodes to each other. The general form is:
- sum(alpha(i)*du(j)_node(k)) = beta, sum over all alpha(i)
- alpha(i) is a given value, du(j)_node(k) is the displacement compoment du(j) (dx, dy, dz, and possibly drx, dry, drz) of node k
- To restrict the displacement of a node in a given direction, we use the fact that the inner product of two perpendicular vectors are zero. Hence defining a vector tau perpendicular to an allowed displacement is suitable for this purpose.
- from the prevous it follows that tau is a vector in the restricted direction.
- For this case: we define a vector tau in the radial direction of the cylinder and the allowable displacement is perpendicular to this, ie tangential to the cylinder surface. Now make this concrete with the LIAISON_DDL keyword.
- The LIAISON_DDL keyword with NOEUD defines the following restriction for the slightly more general case:
- LIAISON_DDL=(_F(NOEUD=('Ni','Nj'),DDL=('DY','DZ'),COEF_MULT=(alpha1,alpha2),COEF_IMPO=beta),
- component 'DY' of node 'Ni' --> dy_ni
- component 'DZ' of node 'Nj' --> dz_nj
- alpha1*dy_ni + alpha2*dz_nj = beta
So if we choose alpha1 and alpha2 the components of the vector from the centre of the cylindrical hole to the node itself,
say tau = [tau1,tau2] = [(ynode-yc), (znode-zc)], - [alpha1,alpha2] = [tau1,tau2] and putting beta to zero we have:
- <[tau1,tau2],[DY,DZ]>=0.
- This means we effectively have: tau1*DY=-tau2*DZ or a tangential restriction of the movement around the cylinder centre.
- LIAISON_DDL=(_F(NOEUD=('Ni','Nj'),DDL=('DY','DZ'),COEF_MULT=(alpha1,alpha2),COEF_IMPO=beta),
Remarks
- the selection of the cylindrical coordinate system is severly restricted here. It should be relatively easy to implement a general cylindrical coordinate system by selecting two points that define the centre axis and a radius. But I really need to activate Numpy first before I try. Although it should work on Salome5.1.4 ...
- the output by the Python script needs to be improved. That I will call work in progress.
- In stead of looping through all the nodes, preselecting the nodes on the cylinder by a filter (already done) is possible. But at the moment I cannot access the selected nodes (again that means I donot know how to, see also Python script at the end).
Selection of the nodes (Salome)
As usual, the geometry and mesh is defined by a Python script. This makes it relatively easy to select the nodes on the cylinder face. The following snippet of the code selects the nodes and writes the keyword and corresponding LIAISON_DDL part to the Python output window:
... y1 = 2 z1 = 3 ... Rcylinder = 0.45 Lcylinder = 2.00 ... countall=0 countcyl=0 node_oncyl = [] Rcyl2 = Rcylinder**2 eps = 1e-4 print "\ncopy and paste next part into command file", print "\n******************************************\n" print "cylco2=AFFE_CHAR_MECA(MODELE=Cmod,\n LIAISON_DDL = (", for allnode in Mbcyl.GetNodesId(): countall+=1 q1, q2, q3 = Mbcyl.GetNodeXYZ(allnode) # on cylinder yes/no: dy2 = (q2-y1)**2 dz2 = (q3-z1)**2 if abs(dy2+dz2-Rcyl2)<eps: countcyl+=1 cnode = to_string(allnode) node_oncyl.append(allnode) #ss = [] ss = "\n _F(NOEUD=('N"+cnode+"','N"+cnode+"'),DDL=('DY','DZ'),COEF_MULT=(" print ss, print (q2-y1), print ",", print (q3-z1), print "),COEF_IMPO=0.00),", pass pass print "),);" print "\nend of copy and past", print "\n********************\n\n\n"
the main part of the code print output to the Python console in a suitable formatted way.
Copy and paste the indicated output from the Python console into the C-A command file:
All should run fine. Well, it did here after more than a few tries. It looks like:
cylco2=AFFE_CHAR_MECA(MODELE=Cmod, LIAISON_DDL = _F(NOEUD=('N162 ','N162 '),DDL=('DY','DZ'),COEF_MULT=( 0.312402110475 , -10.3238902922 ),COEF_IMPO=0.00), _F(NOEUD=('N12 ','N12 '),DDL=('DY','DZ'),COEF_MULT=( 0.0 , -9.55 ),COEF_IMPO=0.00), _F(NOEUD=('N695 ','N695 '),DDL=('DY','DZ'),COEF_MULT=(-0.443605677987 , -9.92440897899 ),COEF_IMPO=0.00), ... _F(NOEUD=('N765 ','N765 '),DDL=('DY','DZ'),COEF_MULT=(-0.162737538758 , -10.4195431962 ),COEF_IMPO=0.00), _F(NOEUD=('N551 ','N551 '),DDL=('DY','DZ'),COEF_MULT=( 0.107926195474 , -9.56313396066 ),COEF_IMPO=0.00),),);
The results of the calculation
The following picture show the rotation of the block around the cylindrical hole:
- The left part of the picture displays the von Mises stresses in the construction: they are zero to the working precision as we expect. The rotation around the cylinder is controlled by the displacement at the top edge and can be performed by a rigid body movement.
- The right part shows the displacement in y direction: the construction nicely rotates around the cylindrical hole.
A bit of fooling around ...
Selecting the same nodes around the cylindrical hole but choosing the rotating point with an offset of 10 in z direction yields the following displacement:
*![]()
The code to generate this during the Salome session:
... zoffset = 10.0 ## change to 10.0 for second picture ... for allnode in Mbcyl.GetNodesId(): countall+=1 q1, q2, q3 = Mbcyl.GetNodeXYZ(allnode) # on cylinder yes/no: dy2 = (q2-y1)**2 dz2 = (q3-z1)**2 if abs(dy2+dz2-Rcyl2)<eps: countcyl+=1 cnode = to_string(allnode) node_oncyl.append(allnode) #ss = [] ss = "\n _F(NOEUD=('N"+cnode+"','N"+cnode+"'),DDL=('DY','DZ'),COEF_MULT=(" print ss, print (q2-y1), print ",", print (q3-z1-zoffset), print "),COEF_IMPO=0.00),", pass pass
Input files for the FE Analysis and references
Input files:
- Python script for defining the geometry and mesh (blockrot4.py), load by File --> load script (cntrl T in the object browser), refresh (F5) after running
- Save the mesh file by right clicking in the object browser by right clicking on the mesh name Mbcyl; select Export to MED and accept or change the default values
- ASTK file (cylrot.astk, you need to edit the path to your requirements ...)
- command file (blockb.comm)
Download the files here:
- blockrot4.py and blockrot6.py basicly perform the same operation: selecting nodes on the cylinder. Though the way it is implemented is quite different: blockrot4 depends on the measures of the geometry whereas blockrot6 directly selects the nodes on the cylinder area by filters. See also [Salome] for more reference (and thank you JMB for the question asked there).
- blockrot7.py writes the liaison_ddl information to a file. Remember to change the path at line 109 to a suitable value for your system.
Reference:
[AFFE_CHAR_MECA, page 20/109]