You are here
apply periodic boundary condition
Tue, 2015-06-23 23:32 - Sarah A. Eassa
hi
could any one please tell me how to apply periodic boundary conditions for off axis loading on unidirectional fiber composite.
I am trying to apply the PBC in this paper "Yunfa Zhang ,Zihui Xia,Micromechanical Analysis of Interphase Damage for Fiber Reinforced Composite
Laminates", but the model does not seam right.
Could any one please tell me the steps that I shall take
Forums:
re: PBCs
If you're using Abaqus it's actually quite simple. Use the equations constraints (Seatrch for *EQUATION keyword in the manual) for nodes that are opposite facewise along with a common reference point. Now what does that mean, let's say your geometry is cubic. Think about just the front and rear face. Imagine you meshed it such that each face 9 nodes. 4 corner nodes, 4 mid edge nodes and one center node. Now you would make a reference point and then join the displacement degrees of freedom of the reference point to each of the nodes on the front face and the rear face using equation constraints so it looks like u1(rear)-u1(front)+u1(reference)=0, and similarly u2(read )-u2(front)+u2(reference point)= 0 and so on and so forth.
I've found the best way to do it is using python scripting and sweep meshing so that the nodes on each opposite face look the same. I make sets for each opposite pair of nodes and then script the equations. If you're interested I might make post that python script.
Ashwin
Hello Ashwin please share the
Hello Ashwin please share the python script for PBC
thanks in advance
Ravi
thank you
Thank you so much for your reply.
I have already written the constraint equation in the input file .inp. for the inner nodes of opposite faces .
But in the paper I mentioned, do i have to put sigma23=0 as a boundary condition? Can I put e23=0 instead?
Because I can not put a load equals zero as a boundary condition. I have to use a user defined fileld, which I am not familiar with.
and how can I put e23 to be zero?
sorry if my questions are not making sence, but I am new in this field :)
sweep meshing
hello ashwin
could you describe how to sweep mesh if say there are 3 layers of materials making a cube.
plz help
Applying periodic boundary condition
Attached paper supplies detailed procedure and necessary scipt for applying PBC.
http://imechanica.org/files/pbc.pdf
--
Anshul
Thanks alot Anshulf for your
Thanks alot Anshulf for your help.
I have a quastion though regarding this paper.
If I am using a 3D, do i have to deal with nodes of inner faces separated from nodes of inner edges seperated from nodes at the corners,
as in the book of finote element of Barbero?
Thank you once again :)
Function in Abaqus Script
I have recently uploaded a function to apply periodic boundary condtions written in Abaqus Script (including examples) on my website, you can find it here:
http://www.overvelde.com/#downloads
The same page also contains a short introduction to Abaqus script if you are not familiar with it.
You can also direclty download the files through the following link:
http://www.overvelde.com/wp-content/uploads/2015/02/Periodic-Boundary-Condtions-Script-Abaqus.zip
Any comments to improve the function are appreciated!
thank you for sharing
thank you so much for sharing your script. I appreciate it :)
the limits
Dear Johannes
why the code doesn't work for more elements ? what if I apply a quadratic structured mesh ?
The code will run, but
The code will run, but probably slowly. What you can do is define a set with only those nodes that need to be included in the periodic boundary conditions. The code will loop over all nodes in the set to determine if it has another node in that set that is located periodically, so the less nodes the faster.
PBCs
Haha, I see Mr. Overvelde beat me, not only in answering your question more throughly but also doing a much better job than me with providing a script. Interestingly enough Johannes, you beat me with the Python scripting document as well. I made a similar document in my undergrad which I would distribute among friends but just before I could publish it, I saw your website and your document. Now I just redirect people to your document!
Nonetheless, here is my script as well, this works for 3D implementation (though for 2D all you have to do is remove some lines), here is the matlab code which writes the python code.
Conditions:- 1) Your model is cubic with the dimension set by the variable 'ulimit'.
2) You have meshed your geometry using sweep meshing or structured mesh.
3) You have named your part Part-1 (the default name)
Change the value of ulimit and run the below code in matlab.
Copy and paste the following code in Abaqus' message area.
%ulimit=size of the cubic geometry, for example if it's 10
ulimit=10;
%PERIODIC BC Implementation
fprintf('Nodes=mdb.models[''Model-1''].rootAssembly.instances[''Part-1-1''].nodes\n');
fprintf('k=0\n')
fprintf('ctr=0\n')
fprintf('\nfor i in range(len(Nodes)):\n');
fprintf(' Cx = Nodes[i].coordinates[0]\n');
fprintf(' Cz = Nodes[i].coordinates[2]\n');
fprintf(' if Cx==0 and Cz>=0.01 and Cz<=%d-0.01:\n',ulimit);
fprintf(' k=k+1\n')
fprintf(' ctr=ctr+1\n')
fprintf(' Nodetemp=mdb.models[''Model-1''].rootAssembly.instances[''Part-1-1''].nodes.getByBoundingBox(Cx+%d-0.01,Nodes[i].coordinates[1]-0.01,Nodes[i].coordinates[2]-0.01,Cx+%d+0.01,Nodes[i].coordinates[1]+0.01,Nodes[i].coordinates[2]+0.01)\n',ulimit,ulimit);
fprintf(' mdb.models[''Model-1''].rootAssembly.Set(name=''simple%%d'' %%k,nodes=Nodetemp)\n')
fprintf(' highlight(mdb.models[''Model-1''].rootAssembly.sets[''simple%%d'' %%k])\n')
fprintf(' k=k+1\n')
fprintf(' mdb.models[''Model-1''].rootAssembly.Set(name=''simple%%d'' %%k,nodes=Nodes[i:i+1])\n')
fprintf(' mdb.models[''Model-1''].Equation(name=''Constraint-%%d'' %%ctr, terms=((-1.0, ''simple%%d'' %%(k-1), 1), (1.0, ''simple%%d'' %%k, 1), (1.0, ''ReferencePt2'', 1)))\n');
fprintf(' ctr=ctr+1\n')
fprintf(' mdb.models[''Model-1''].Equation(name=''Constraint-%%d'' %%ctr, terms=((-1.0, ''simple%%d'' %%(k-1), 2), (1.0, ''simple%%d'' %%k, 2), (1.0, ''ReferencePt2'', 2)))\n');
fprintf(' ctr=ctr+1\n');
fprintf(' mdb.models[''Model-1''].Equation(name=''Constraint-%%d'' %%ctr, terms=((-1.0, ''simple%%d'' %%(k-1), 3), (1.0, ''simple%%d'' %%k, 3), (1.0, ''ReferencePt2'', 3)))\n');
fprintf(' highlight(Nodes[i])\n')
fprintf('\nfor i in range(len(Nodes)):\n');
fprintf(' Cx = Nodes[i].coordinates[0]\n');
fprintf(' Cz = Nodes[i].coordinates[2]\n');
fprintf(' if Cx==0 and Cz>=0.01 and Cz<=%d-0.01:\n',ulimit);
fprintf(' Nodetemp=mdb.models[''Model-1''].rootAssembly.instances[''Part-1-1''].nodes.getByBoundingBox(Cx+%d-0.01,Nodes[i].coordinates[1],Nodes[i].coordinates[2],Cx+%d+0.01,Nodes[i].coordinates[1],Nodes[i].coordinates[2])\n',ulimit,ulimit);
fprintf(' mdb.models[''Model-1''].rootAssembly.Set(name=''junkset'',nodes=Nodetemp)\n')
fprintf(' unhighlight(Nodes[i])\n')
fprintf(' unhighlight(mdb.models[''Model-1''].rootAssembly.sets[''junkset''])\n')
fprintf('\nfor i in range(len(Nodes)):\n');
fprintf(' Cy = Nodes[i].coordinates[1]\n');
fprintf(' Cx = Nodes[i].coordinates[0]\n');
fprintf(' Cz = Nodes[i].coordinates[2]\n');
fprintf(' if Cy==0 and Cz>=0.01 and Cz<=%d-0.01 and Cx>=0.01 and Cx<=%d-0.01:\n',ulimit,ulimit);
fprintf(' k=k+1\n')
fprintf(' ctr=ctr+1\n')
fprintf(' Nodetemp=mdb.models[''Model-1''].rootAssembly.instances[''Part-1-1''].nodes.getByBoundingBox(Nodes[i].coordinates[0]-0.01,Cy+%d-0.01,Nodes[i].coordinates[2]-0.01,Nodes[i].coordinates[0]+0.01,Cy+%d+0.01,Nodes[i].coordinates[2]+0.01)\n',ulimit,ulimit);
fprintf(' mdb.models[''Model-1''].rootAssembly.Set(name=''simple%%d'' %%k,nodes=Nodetemp)\n')
fprintf(' highlight(mdb.models[''Model-1''].rootAssembly.sets[''simple%%d'' %%k])\n')
fprintf(' k=k+1\n')
fprintf(' mdb.models[''Model-1''].rootAssembly.Set(name=''simple%%d'' %%k,nodes=Nodes[i:i+1])\n')
fprintf(' highlight(Nodes[i])\n')
fprintf(' mdb.models[''Model-1''].Equation(name=''Constraint-%%d'' %%ctr, terms=((-1.0, ''simple%%d'' %%(k-1), 1), (1.0, ''simple%%d'' %%k, 1), (1.0, ''ReferencePt3'', 1)))\n');
fprintf(' ctr=ctr+1\n')
fprintf(' mdb.models[''Model-1''].Equation(name=''Constraint-%%d'' %%ctr, terms=((-1.0, ''simple%%d'' %%(k-1), 2), (1.0, ''simple%%d'' %%k, 2), (1.0, ''ReferencePt3'', 2)))\n');
fprintf(' ctr=ctr+1\n')
fprintf(' mdb.models[''Model-1''].Equation(name=''Constraint-%%d'' %%ctr, terms=((-1.0, ''simple%%d'' %%(k-1), 3), (1.0, ''simple%%d'' %%k, 3), (1.0, ''ReferencePt3'', 3)))\n');
fprintf('\nfor i in range(len(Nodes)):\n');
fprintf(' Cy = Nodes[i].coordinates[1]\n');
fprintf(' Cz = Nodes[i].coordinates[2]\n');
fprintf(' if Cy==0 and Cz>=0.01 and Cz<=%d-0.01:\n',ulimit);
fprintf(' Nodetemp=mdb.models[''Model-1''].rootAssembly.instances[''Part-1-1''].nodes.getByBoundingBox(Nodes[i].coordinates[0],Cy+%d-0.01,Nodes[i].coordinates[2],Nodes[i].coordinates[0],Cy+%d+0.01,Nodes[i].coordinates[2])\n',ulimit,ulimit);
fprintf(' mdb.models[''Model-1''].rootAssembly.Set(name=''junkset'',nodes=Nodetemp)\n')
fprintf(' unhighlight(Nodes[i])\n')
fprintf(' unhighlight(mdb.models[''Model-1''].rootAssembly.sets[''junkset''])\n')
fprintf('\nfor i in range(len(Nodes)):\n');
fprintf(' Cz = Nodes[i].coordinates[2]\n');
fprintf(' if Cz<=0.001:\n');
fprintf(' k=k+1\n')
fprintf(' ctr=ctr+1\n')
fprintf(' Nodetemp=mdb.models[''Model-1''].rootAssembly.instances[''Part-1-1''].nodes.getByBoundingBox(Nodes[i].coordinates[0]-0.01,Nodes[i].coordinates[1]-0.01,Cz+%d-0.01,Nodes[i].coordinates[0]+0.01,Nodes[i].coordinates[1]+0.01,Cz+%d+0.01)\n',ulimit,ulimit);
fprintf(' mdb.models[''Model-1''].rootAssembly.Set(name=''simple%%d'' %%k,nodes=Nodetemp)\n')
fprintf(' highlight(mdb.models[''Model-1''].rootAssembly.sets[''simple%%d'' %%k])\n')
fprintf(' k=k+1\n')
fprintf(' mdb.models[''Model-1''].rootAssembly.Set(name=''simple%%d'' %%k,nodes=Nodes[i:i+1])\n')
fprintf(' highlight(Nodes[i])\n')
fprintf(' mdb.models[''Model-1''].Equation(name=''Constraint-%%d'' %%ctr, terms=((-1.0, ''simple%%d'' %%(k-1), 1), (1.0, ''simple%%d'' %%k, 1), (1.0, ''ReferencePt1'', 1)))\n');
fprintf(' ctr=ctr+1\n')
fprintf(' mdb.models[''Model-1''].Equation(name=''Constraint-%%d'' %%ctr, terms=((-1.0, ''simple%%d'' %%(k-1), 2), (1.0, ''simple%%d'' %%k, 2), (1.0, ''ReferencePt1'', 2)))\n');
fprintf(' ctr=ctr+1\n')
fprintf(' mdb.models[''Model-1''].Equation(name=''Constraint-%%d'' %%ctr, terms=((-1.0, ''simple%%d'' %%(k-1), 3), (1.0, ''simple%%d'' %%k, 3), (1.0, ''ReferencePt1'', 3)))\n');
fprintf('\nfor i in range(len(Nodes)):\n');
fprintf(' Cz = Nodes[i].coordinates[2]\n');
fprintf(' if Cz<=0.001:\n');
fprintf(' Nodetemp=mdb.models[''Model-1''].rootAssembly.instances[''Part-1-1''].nodes.getByBoundingBox(Nodes[i].coordinates[0],Nodes[i].coordinates[1],Cz+%d-0.01,Nodes[i].coordinates[0],Nodes[i].coordinates[1],Cz+%d+0.01)\n',ulimit,ulimit);
fprintf(' mdb.models[''Model-1''].rootAssembly.Set(name=''junkset'',nodes=Nodetemp)\n')
fprintf(' unhighlight(Nodes[i])\n')
fprintf(' unhighlight(mdb.models[''Model-1''].rootAssembly.sets[''junkset''])\n\n')
These are all single line fprintf commands, so after you copy paste it, make sure to put the fprintf commands back to the original single line format or you'll get print errors.
Following is the python code, WARNING: In Python, indentations are very important. Make sure after you copy and paste it, the code format is conserved. This following code works for a cube of size 10 by 10 by 10 units. Change the value in the below code for different cubic sizes.
Nodes=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].nodes
k=0
ctr=0
for i in range(len(Nodes)):
Cx = Nodes[i].coordinates[0]
Cz = Nodes[i].coordinates[2]
if Cx==0 and Cz>=0.01 and Cz<=10-0.01:
k=k+1
ctr=ctr+1
Nodetemp=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].nodes.getByBoundingBox(Cx+10-0.01,Nodes[i].coordinates[1]-0.01,Nodes[i].coordinates[2]-0.01,Cx+10+0.01,Nodes[i].coordinates[1]+0.01,Nodes[i].coordinates[2]+0.01)
mdb.models['Model-1'].rootAssembly.Set(name='simple%d' %k,nodes=Nodetemp)
highlight(mdb.models['Model-1'].rootAssembly.sets['simple%d' %k])
k=k+1
mdb.models['Model-1'].rootAssembly.Set(name='simple%d' %k,nodes=Nodes[i:i+1])
mdb.models['Model-1'].Equation(name='Constraint-%d' %ctr, terms=((-1.0, 'simple%d' %(k-1), 1), (1.0, 'simple%d' %k, 1), (1.0, 'ReferencePt2', 1)))
ctr=ctr+1
mdb.models['Model-1'].Equation(name='Constraint-%d' %ctr, terms=((-1.0, 'simple%d' %(k-1), 2), (1.0, 'simple%d' %k, 2), (1.0, 'ReferencePt2', 2)))
ctr=ctr+1
mdb.models['Model-1'].Equation(name='Constraint-%d' %ctr, terms=((-1.0, 'simple%d' %(k-1), 3), (1.0, 'simple%d' %k, 3), (1.0, 'ReferencePt2', 3)))
highlight(Nodes[i])
for i in range(len(Nodes)):
Cx = Nodes[i].coordinates[0]
Cz = Nodes[i].coordinates[2]
if Cx==0 and Cz>=0.01 and Cz<=10-0.01:
Nodetemp=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].nodes.getByBoundingBox(Cx+10-0.01,Nodes[i].coordinates[1],Nodes[i].coordinates[2],Cx+10+0.01,Nodes[i].coordinates[1],Nodes[i].coordinates[2])
mdb.models['Model-1'].rootAssembly.Set(name='junkset',nodes=Nodetemp)
unhighlight(Nodes[i])
unhighlight(mdb.models['Model-1'].rootAssembly.sets['junkset'])
for i in range(len(Nodes)):
Cy = Nodes[i].coordinates[1]
Cx = Nodes[i].coordinates[0]
Cz = Nodes[i].coordinates[2]
if Cy==0 and Cz>=0.01 and Cz<=10-0.01 and Cx>=0.01 and Cx<=10-0.01:
k=k+1
ctr=ctr+1
Nodetemp=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].nodes.getByBoundingBox(Nodes[i].coordinates[0]-0.01,Cy+10-0.01,Nodes[i].coordinates[2]-0.01,Nodes[i].coordinates[0]+0.01,Cy+10+0.01,Nodes[i].coordinates[2]+0.01)
mdb.models['Model-1'].rootAssembly.Set(name='simple%d' %k,nodes=Nodetemp)
highlight(mdb.models['Model-1'].rootAssembly.sets['simple%d' %k])
k=k+1
mdb.models['Model-1'].rootAssembly.Set(name='simple%d' %k,nodes=Nodes[i:i+1])
highlight(Nodes[i])
mdb.models['Model-1'].Equation(name='Constraint-%d' %ctr, terms=((-1.0, 'simple%d' %(k-1), 1), (1.0, 'simple%d' %k, 1), (1.0, 'ReferencePt3', 1)))
ctr=ctr+1
mdb.models['Model-1'].Equation(name='Constraint-%d' %ctr, terms=((-1.0, 'simple%d' %(k-1), 2), (1.0, 'simple%d' %k, 2), (1.0, 'ReferencePt3', 2)))
ctr=ctr+1
mdb.models['Model-1'].Equation(name='Constraint-%d' %ctr, terms=((-1.0, 'simple%d' %(k-1), 3), (1.0, 'simple%d' %k, 3), (1.0, 'ReferencePt3', 3)))
for i in range(len(Nodes)):
Cy = Nodes[i].coordinates[1]
Cz = Nodes[i].coordinates[2]
if Cy==0 and Cz>=0.01 and Cz<=10-0.01:
Nodetemp=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].nodes.getByBoundingBox(Nodes[i].coordinates[0],Cy+10-0.01,Nodes[i].coordinates[2],Nodes[i].coordinates[0],Cy+10+0.01,Nodes[i].coordinates[2])
mdb.models['Model-1'].rootAssembly.Set(name='junkset',nodes=Nodetemp)
unhighlight(Nodes[i])
unhighlight(mdb.models['Model-1'].rootAssembly.sets['junkset'])
for i in range(len(Nodes)):
Cz = Nodes[i].coordinates[2]
if Cz<=0.001:
k=k+1
ctr=ctr+1
Nodetemp=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].nodes.getByBoundingBox(Nodes[i].coordinates[0]-0.01,Nodes[i].coordinates[1]-0.01,Cz+10-0.01,Nodes[i].coordinates[0]+0.01,Nodes[i].coordinates[1]+0.01,Cz+10+0.01)
mdb.models['Model-1'].rootAssembly.Set(name='simple%d' %k,nodes=Nodetemp)
highlight(mdb.models['Model-1'].rootAssembly.sets['simple%d' %k])
k=k+1
mdb.models['Model-1'].rootAssembly.Set(name='simple%d' %k,nodes=Nodes[i:i+1])
highlight(Nodes[i])
mdb.models['Model-1'].Equation(name='Constraint-%d' %ctr, terms=((-1.0, 'simple%d' %(k-1), 1), (1.0, 'simple%d' %k, 1), (1.0, 'ReferencePt1', 1)))
ctr=ctr+1
mdb.models['Model-1'].Equation(name='Constraint-%d' %ctr, terms=((-1.0, 'simple%d' %(k-1), 2), (1.0, 'simple%d' %k, 2), (1.0, 'ReferencePt1', 2)))
ctr=ctr+1
mdb.models['Model-1'].Equation(name='Constraint-%d' %ctr, terms=((-1.0, 'simple%d' %(k-1), 3), (1.0, 'simple%d' %k, 3), (1.0, 'ReferencePt1', 3)))
for i in range(len(Nodes)):
Cz = Nodes[i].coordinates[2]
if Cz<=0.001:
Nodetemp=mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].nodes.getByBoundingBox(Nodes[i].coordinates[0],Nodes[i].coordinates[1],Cz+10-0.01,Nodes[i].coordinates[0],Nodes[i].coordinates[1],Cz+10+0.01)
mdb.models['Model-1'].rootAssembly.Set(name='junkset',nodes=Nodetemp)
unhighlight(Nodes[i])
unhighlight(mdb.models['Model-1'].rootAssembly.sets['junkset'])
Edit
Another condition is to have have 3 Reference points. Make sets out of the Reference points as well, so you may use the following code in Matlab. The creating Set's portion of the code may or may not work for you depending on which array position your referemce point is. If it doesn't work, do it manually, call the Z axis reference point ReferencePt1 , the X axis ReferencePt2 and the Y axis ReferencePt3
fprintf('mdb.models[''Model-1''].rootAssembly.ReferencePoint(point=(%f, %f, %f))\n',ulimit/2,ulimit/2,2*ulimit);
fprintf('mdb.models[''Model-1''].rootAssembly.ReferencePoint(point=(%f, %f, %f))\n',ulimit*2,ulimit/2,ulimit/2);
fprintf('mdb.models[''Model-1''].rootAssembly.ReferencePoint(point=(%f, %f, %f))\n',ulimit/2,ulimit*2,ulimit/2);
fprintf('mdb.models[''Model-1''].rootAssembly.Set(name=''ReferencePt1'', referencePoints=(mdb.models[''Model-1''].rootAssembly.referencePoints[4], ))\n');
fprintf('mdb.models[''Model-1''].rootAssembly.Set(name=''ReferencePt2'', referencePoints=(mdb.models[''Model-1''].rootAssembly.referencePoints[5], ))\n');
fprintf('mdb.models[''Model-1''].rootAssembly.Set(name=''ReferencePt3'', referencePoints=(mdb.models[''Model-1''].rootAssembly.referencePoints[6], ))\n');
Thank you so much for your
Thank you so much for your detailed script.
I have another question if you don't mind.
I would like to apply e2, e3 and e23 on a representative volume element of a fiber composite.
I use a reference point to apply the strains, fo example reference point 1 for the e1, and reference point 2 for e2.
The question is that, how can I put e23 on a reference point 3, is it rotation around 1? will it be UR1 on reference point 3?
Thanks in advance
uncoupled thermal-stress analysis using periodic boundary condi
hi
could any one please tel me if I shall apply periodic boundary condition in the two sequential analysis for uncoupled thermal-stress analysis or not?
Thanks
hi
hi
I have a question.
If the mesh is the same for two surfaces, can I use inner node set for each of the faces to apply the constraint equations instead of applying them to each pair of nodes?
It will really save effort and time :)
Thanks
Thermal residual stress
hi
Could any one please tell me how to apply periodic boundary conditions for thermal analysis of fiber composite.
I am getting the stresses from the difference of the thermal expansion coefficients for both the fiber and the matrix.
What would be the dispalcement of the reference point??????
Thanks