Commit 59c5b73d authored by matecellab's avatar matecellab
Browse files

changed line numbering in documentation

parent 3a62d1f3
% !TeX root = BasilLab.tex
e% !TeX root = BasilLab.tex
\begin{savequote}
I can't be as confident about computer science as I can about biology. Biology easily has 500\,years of exciting problems to work on. It's at that level.
\qauthor{Donald Knuth, Computer Literacy Bookshops Interview (1993)}
......@@ -69,7 +69,7 @@ The Poisson's ratio $\nu$ can be calculated from $\mu_0$ and $K_0$ with (see sec
\nu = \frac{3\frac{K_0}{\mu_0}-2}{6\frac{K_0}{\mu_0}+2}
\end{equation}
BasilLab provides the function \bas{D2nu} to calculate the Poisson's ratio $\nu$ from $\mu_0$ and $D$, and the function \bas{nu2D} to calculate $D$ from $\mu_0$ and a Poisson's ratio $\nu$:
\begin{lstlisting}[language=basillab, numbers=none]
\begin{lstlisting}[language=basillab, firstnumber=last]
>> nu = bas_D2nu(0.0156, 60)
nu =
......@@ -122,7 +122,7 @@ param.indOff = 0.000000;
Next to the Abaqus python script \texttt{basil<basilid>\_job.py\index{basil<basilid>\_job.py@\texttt{basil<basilid>\_job.py}}} that builds the model, \bas{buildFEM} also writes a file \texttt{basil<basilid>\_modelparam.m\index{basil<basilid>\_modelparam.m@\texttt{basil<basilid>\_modelparam.m}}} in the current working directory. This file contains all the optional parameters that \bas{buildFEM} used to write the Abaqus python script. The default values are substituted for parameters that were not explicitly defined in the input parameter structure. In addition, the first line of \texttt{basil<basilid>\_modelparam.m\index{basil<basilid>\_modelparam.m@\texttt{basil<basilid>\_modelparam.m}}} contains the the date and time that \bas{buildFEM} was called to write \texttt{basil<basilid>\_job.py\index{basil<basilid>\_job.py@\texttt{basil<basilid>\_job.py}}} (\texttt{buildon}, see listing \ref{modelparam} for an example).
With this file, you can always and easily check the parameters that were used to build a particular model. As long as you keep this file, you can always ask \bas{buildFEM} to build the model again from scratch and be sure that you use the same options as you used before. Likewise, you can be sure to use the same options for (certain) different models. E.g.\ to build the model for \texttt{<basilid> = 140611} with the same parameters that were used to build a model for \texttt{<basilid> = 140613}, run:
\begin{lstlisting}[language=basillab, numbers=none]
\begin{lstlisting}[language=basillab, firstnumber=last]
>> run basil140613_modelparam.m, bas_buildFEM(140611,param)
\end{lstlisting}
......@@ -139,7 +139,7 @@ Furthermore, \bas{rebuildFEM} also accepts the parameter structure as (optional
\end{lstlisting}
This function \bas{rebuildFEM} can be very useful when you wish to change a single or few parameters in all your models, while keeping the (other) original parameters for each individual \texttt{<basilid>} the same as before. E.g.\ when, you call
\begin{lstlisting}[language=basillab, numbers=none]
\begin{lstlisting}[language=basillab, firstnumber=last]
>> param.seed = 0.5; bas_rebuildFEM([140611:140613], param)
\end{lstlisting}
all parameters for \texttt{<basilid> = 140611}, \texttt{<basilid> = 140612} and \texttt{<basilid> = 140613} will be the same as before, except for the mesh seed size.\\
......@@ -197,7 +197,7 @@ The first line contains the 3D coordinates of the (100 for the skin, or 50 for t
\subsection{Initialise the model}
At his point, \bas{buildFEM} starts writing to the main \texttt{basil<basilid>\_job.py\index{basil<basilid>\_job.py@\texttt{basil<basilid>\_job.py}}} file. This file starts with a comment header, imports necessary (Abaqus specific) python modules, and defines a global meshing option:
\begin{lstlisting}[language=apython, numbers=none]
\begin{lstlisting}[language=apython, firstnumber=last]
"""===============================
Script to make an abaqus job of a rat leg with tibia,
indented with cylindrical indenter with spherical head
......@@ -208,24 +208,27 @@ from abaqusConstants import *
from math import *
session.defaultMesherOptions.setValues(allowMapped=TRUE)
\end{lstlisting}
Next, we define our model, and delete the 'default' model \texttt{'Model-1'}:
\begin{lstlisting}[language=apython, numbers=none]
\begin{lstlisting}[language=apython, firstnumber=last]
# Create a model
Mdb()
myModel = mdb.Model(name='Basil140611')
del mdb.models['Model-1']
\end{lstlisting}
\subsection{Model geometry}
We are now ready to import the Abaqus \texttt{part}-module and define our first \texttt{part}, the rat's leg:
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# Create first part, the leg
import part
myPart = myModel.Part(name='leg', dimensionality=THREE_D, type=DEFORMABLE_BODY)
# import contours
execfile('basil140611_meshContours.py')
\end{lstlisting}
This last line \texttt{5} imports the previously written \texttt{basil<basilid>\_meshContours.py\index{basil<basilid>\_meshContours.py@\texttt{basil<basilid>\_meshContours.py}}} and that should result in an equal amount of Abaqus `wires' for the skin and tibia contours. If you run the main \texttt{basil<basilid>\-\_job\-.py\index{basil<basilid>\_job.py@\texttt{basil<basilid>\_job.py}}}-script up until this point, Abaqus should display something like figure \eqref{abaqus_meshContours}.
\begin{figure}[h!]
......@@ -234,7 +237,7 @@ This last line \texttt{5} imports the previously written \texttt{basil<basilid>\
\end{figure}
Next, we convert the wires to two surfaces: one for the outer skin wires, and one for the tibia wires. This surface fitting through the wires is called 'lofting' in Abaqus.
\begin{lstlisting}[language=apython, , escapechar=!]
\begin{lstlisting}[language=apython, firstnumber=last, escapechar=!]
# loft skin
e = myPart.edges
# get ordered list of egdes (because: not the same as wire-numbers, and prone to change when the part changes)
......@@ -246,6 +249,7 @@ for i in range(0, 19):
ledges.append( t[0] );
# create loft of edge list
myPart.ShellLoft(loftsections=((ledges[0], ), (ledges[1], ), (ledges[2], ), (ledges[3], ), (ledges[4], ), (ledges[5], ), (ledges[6], ), (ledges[7], ), (ledges[8], ), (ledges[9], ), (ledges[10], ), (ledges[11], ), (ledges[12], ), (ledges[13], ), (ledges[14], ), (ledges[15], ), (ledges[16], ), (ledges[17], ), (ledges[18], ), ), startCondition=NONE,endCondition=NONE)
\end{lstlisting}
In order to `loft' a number of wires, we need to know the internal \textsl{edge} numbering that Abaqus assigned to our wires. This numbering is not trivial (even appears to be arbitrary) and the internal numbering of the edges changes when anything happens with any wire: i.e.\ when the skin wires/edges are lofted, the internal numbering of the tibia wires/edges will have changed.
......@@ -260,7 +264,7 @@ When the outer skin contours have been lofted, we use the same approach to find
\end{figure}
The model now needs to be 'closed': we need to surface loft the two open sides of the model to create a closed volume. Similar to the lofting of the skin and tibia surfaces, we need to query the edge numbers (lines \texttt{4} and \texttt{9}, again: book-keeping by \bas{buildFEM}) of the sides that we wish to loft (line \texttt{14}):
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=46]
# close lofts
e = myPart.edges
# find edge on skin
......@@ -275,6 +279,7 @@ t = eb[0]
bedges.append( t[0] )
# close this side
myPart.ShellLoft(loftsections=((bedges[0], ), (sedges[0], )), startCondition=NONE,endCondition=NONE)
\end{lstlisting}
We do this twice, once for each side of the model. If you run the main \texttt{basil<basilid>\-\_job\-.py\index{basil<basilid>\_job.py@\texttt{basil<basilid>\_job.py}}}-script up until this point, Abaqus should display something like figure \eqref{abaqus_closedLofts}.
\begin{figure}[h!]
......@@ -283,17 +288,18 @@ We do this twice, once for each side of the model. If you run the main \texttt{
\end{figure}
Finally, we can tell Abaqus to turn the volume enclosed by the 4 `faces' (tibia surface, skin surface, and the two plain end-surfaces) into a true `volume'-part:
\begin{lstlisting}[language=apython, numbers=none]
\begin{lstlisting}[language=apython, firstnumber=74]
# create volume
f = myPart.faces
myPart.AddCells(faceList = f[0:3])
\end{lstlisting}
And with that we have defined the geometry of the 3D deformable part, the leg.
\subsection{Sets, material and section}
It is now useful to assign skin and tibia sets. These sets will be based on the geometry and are independent of meshing choices later on. E.g.\ all nodes that will later be placed on the \texttt{bone}-set, will be part of the \texttt{bone} node-set. For the skin, we also define a surface that may be used for contact modelling (with the alginate cast) later on. For the bone, we also define a `reference point' (lines \texttt{10--12}), always at coordinates \texttt{(0, 0, 0)}. The bone set will be transformed to a `rigid surface' later on and boundary conditions will be applied to the bone's reference point.
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# add skin set
faces = f.getSequenceFromMask(mask=('[#8 ]', ), )
myPart.Set(name='skin', faces=faces)
......@@ -306,19 +312,21 @@ myPart.Set(name='bone', faces=faces)
myPart.ReferencePoint(point=(0, 0, 0))
lastRP = myPart.referencePoints.values()
myPart.Set(referencePoints=lastRP, name='boneRP')
\end{lstlisting}
The definition of the material parameters for the soft tissue is straightforward:
\begin{lstlisting}[language=apython, numbers=none]
\begin{lstlisting}[language=apython, firstnumber=last]
# define material
import material
myModel.Material(name='softTissue')
myModel.materials['softTissue'].Density(table=((1e-9, ), ))
myModel.materials['softTissue'].Hyperelastic(testData=OFF, type=OGDEN, volumetricResponse=VOLUMETRIC_DATA, table=((0.003600, 5.000000, 57.471264), ))
\end{lstlisting}
Finally, we define a (homogeneous solid) `section' for the soft tissue (leg) with the material properties that we have just defined (line \texttt{4}). The volume that represents the soft tissue is assigned to a new set (lines \texttt{5--7}), and this set is assigned to the `section' (lines \texttt{8--9}):
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# section and material assignment
import section
import regionToolset
......@@ -328,12 +336,13 @@ c = p.cells
p.Set(cells=c, name='softTissue')
region = p.sets['softTissue']
p.SectionAssignment(region=region, sectionName='tissueSection', offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
\end{lstlisting}
\subsection{The indenter}
The indenter is modelled in a new \texttt{part} as a analytical rigid surface: it will not be meshed and it requires a reference point to apply boundary conditions. Because the geometry of the indenter is known (given `shape' with `size' \texttt{param.ri}) we can draw it in the Abaqus `sketch' module:
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# build indenter
import sketch
myModel.Part(name='indenter', dimensionality=THREE_D, type=ANALYTIC_RIGID_SURFACE)
......@@ -355,7 +364,7 @@ The sketch itself (lines \texttt{5--7}, figure \ref{indentersketch}) is axisymme
\end{figure}
Next, we assign a surface to the indenter \texttt{part} for contact modelling later on:
\begin{lstlisting}[language=apython, numbers=none]
\begin{lstlisting}[language=apython, firstnumber=last]
# indenter to surface
p = myModel.parts['indenter']
s = p.faces
......@@ -363,11 +372,12 @@ side2Faces = s.getSequenceFromMask(mask=('[#3 ]', ), )
p.Surface(side1Faces=side2Faces, name='indenter')
\end{lstlisting}
and we assign the indenter's reference point to a set for easy future reference:
\begin{lstlisting}[language=apython, numbers=none]
\begin{lstlisting}[language=apython, firstnumber=last]
p = myModel.parts['indenter']
r = p.referencePoints
refPoints=(r[2], )
p.Set(referencePoints=refPoints, name='indenterRP')
\end{lstlisting}
\noindent Note that the indenter is not yet at it's correct position: it is oriented in the \texttt{(0, 1, 0)}-direction with the tip at \texttt{(0, 0, 0)} (figure \ref{abaqusindenter}). Positioning af the indenter \texttt{part} relative to the leg part will be done later when we `instance' the model.
......@@ -377,7 +387,7 @@ p.Set(referencePoints=refPoints, name='indenterRP')
The `assist-disc' has no role at all in the actual simulation and is only used during the building of the model. In the assembly of the model, the edges of the (repositioned) assist-disc are projected onto the skin surface. The contour of that projection will be used to define the area on the skin that is modelled with contact boundary conditions.
The size of the assist-disc scales with the size of the indenter (\texttt{param.ri}, since \bas{buildFEM} does the book-keeping). The assist-disc gets an inner radius of 4\,\texttt{param.ri} and an outer radius of 4.5\,\texttt{param.ri}. It is represented in the 2D sketch as a line at a `height' of \texttt{y = 1\,param.ri}:
\begin{lstlisting}[language=apython, numbers=none]
\begin{lstlisting}[language=apython, firstnumber=last]
# assist disc to help partitioning of skin for contact
myModel.Part(name='assistDisc', dimensionality=THREE_D, type=ANALYTIC_RIGID_SURFACE)
s1 = myModel.ConstrainedSketch(name='discSketch', sheetSize=10.0)
......@@ -386,6 +396,7 @@ s1.ConstructionLine(point1=(0.0, -1.5), point2=(0.0, 0.0))
p = myModel.parts['assistDisc']
p.AnalyticRigidSurfRevolve(sketch=s1)
p.ReferencePoint(point=(0, 0, 0))
\end{lstlisting}
......@@ -400,7 +411,7 @@ p.ReferencePoint(point=(0, 0, 0))
\subsection{Assemble the parts}
We are now going to assemble model: collect the three parts (leg, indenter and assist-disc) and position them relative to each other. A \texttt{part} that is imported into the \texttt{assembly} is called an \texttt{instance}, in Abaqus:
\begin{lstlisting}[language=apython, numbers=none]
\begin{lstlisting}[language=apython, firstnumber=136]
# assemble
import assembly
a1 = myModel.rootAssembly
......@@ -409,7 +420,7 @@ p = myModel.parts['leg']
a1.Instance(name='leg-1', part=p, dependent=OFF)
\end{lstlisting}
The leg is positioned (instanced) `as is'. The indenter needs to be rotated and translated to the correct position:
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
p = myModel.parts['indenter']
a1.Instance(name='indenter-1', part=p, dependent=OFF)
# position indenter
......@@ -426,7 +437,7 @@ a1.translate(instanceList=('indenter-1', ), vector=(-4.811105, -5.031821, 0.5520
The indenter is first rotated from \texttt{(0, 1, 0)} to \texttt{(0, 0, 1)} (line \texttt{6}). Next, \bas{buildFEM} does the book-keeping to rotate the indenter to the correct orientation from \texttt{(0, 0, 1)} (line \texttt{8}); and to translate it to the correct position for indentation (line \texttt{11}). Note that \bas{buildFEM} uses the rotation matrix stored in \texttt{Indenter\_tipNR.txt\index{Indenter\_tipNR.txt@\texttt{Indenter\_tipNR.txt}}} and the \gibbon-function \gib{rot2VecAngle} to calculate the \texttt{AxisDirection} and \texttt{angle} for the \texttt{rotate}-command for Abaqus (lines \texttt{6} and \texttt{8}).\\
\noindent The orientation of the assist-disc depends on the local curvature of the leg close to the indenter. The function \bas{buildFEM} looks for the point in the skin contour definition that is closest to the indenter tip start position. It then finds the closest point on the same contour, and the closest point on an adjacent contour. With these three points, \bas{buildFEM} estimates the local normal vector of the leg underneath the indenter. The assist-disc is therefore positioned 'parallel to the skin':
\begin{lstlisting}[language=apython, numbers=none]
\begin{lstlisting}[language=apython, firstnumber=last]
# instance and position assist disc
a1 = myModel.rootAssembly
p = myModel.parts['assistDisc']
......@@ -438,6 +449,7 @@ a1.rotate(instanceList=('assistDisc-1', ), axisPoint=(0.0, 0.0, 0.0), axisDirect
# translate wire to correct initial position
a1 = myModel.rootAssembly
a1.translate(instanceList=('assistDisc-1', ), vector=(-4.811105, -5.031821, 0.552060))
\end{lstlisting}
The translation vector (last line) is the same as that for the indenter. However, because we sketched the assist-disc at \texttt{y = 1\,param.ri}, the centre of the disc will not end up at the tip of the indenter, but one indenter radius `higher' instead.\\
......@@ -453,7 +465,7 @@ The translation vector (last line) is the same as that for the indenter. However
\subsection{Define area on skin for indenter contact}
The next step is to `partition' the skin in order to define the contact area for the indenter. As mentioned above, this is the (only) purpose of the assist-disc:
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# partition skin for indenter contact region (surface)
a = myModel.rootAssembly
f1 = a.instances['leg-1'].faces
......@@ -466,6 +478,7 @@ a.PartitionFaceByProjectingEdges(faces=pickedFaces, edges=pickedEdges, extendEdg
s1 = a.instances['leg-1'].faces
side1Faces1 = s1.getSequenceFromMask(mask=('[#1 ]', ), )
a.Surface(side1Faces=side1Faces1, name='skinContactIn')
\end{lstlisting}
We select the inner radius of the assist-disc (lines \texttt{2--4}) and project this ring onto the skin (lines \texttt{6--8}, figure \ref{abaqus_skinPartitioned}). Finally, the newly defined area gets a surface definition for contact (lines \texttt{10--12}).
\begin{figure}[h!]
......@@ -476,7 +489,7 @@ We select the inner radius of the assist-disc (lines \texttt{2--4}) and project
\subsection{Meshing the leg}
After partitioning of the skin surface, we can have Abaqus mesh the leg \texttt{instance}\footnote{In Abaqus you can mesh on the \texttt{part} \textsl{before} it is instanced in the model assembly, or on the \texttt{instance} \textsl{after} the part has been imported into the model assembly. We use the latter method because the assist-disc for partitioning the skin surface is not available in the \texttt{part} but only in the assembly (\texttt{instance}).}. First, we need to `seed' the leg: seeds will be placed at the edges in the model with a density determined by \texttt{param.seed}. These seeds will serve as the corner points of the elements that form a part of the edge: for linear elements nodes will coincide with the seeds; for quadratic elements there will also be nodes between the seeds.
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# mesh leg
import mesh
a = myModel.rootAssembly
......@@ -491,6 +504,7 @@ elemType2 = mesh.ElemType(elemCode=C3D15, elemLibrary=STANDARD)
elemType3 = mesh.ElemType(elemCode=C3D10M, elemLibrary=STANDARD)
a.setElementType(regions=pickedRegions, elemTypes=(elemType1, elemType2, elemType3))
a.generateMesh(regions=partInstances)
\end{lstlisting}
After seeding the leg (lines \texttt{3--5}, figure \ref{abaqus_legSeeded}), we tell Abaqus which meshing technique to use (lines \texttt{6--8}) and which element types can be used for meshing (lines \texttt{10--13}). The mesh is finally generated in line \texttt{14} (figure \ref{abaqus_legMeshed}).\\
......@@ -506,7 +520,7 @@ After seeding the leg (lines \texttt{3--5}, figure \ref{abaqus_legSeeded}), we t
\subsection{Contact modelling for indenter and skin}
Contact between the indenter and the skin has to modelled and this is done in Abaqus' \texttt{interaction} module. First we need to define the type of `interaction' (contact) itself, and then we define the surfaces for which this contact interaction needs to be active:
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# indenter skin contact interaction
import interaction
myModel.ContactProperty('contact')
......@@ -516,6 +530,7 @@ a = myModel.rootAssembly
region1=a.instances['indenter-1'].surfaces['indenter']
region2=a.surfaces['skinContactIn']
myModel.SurfaceToSurfaceContactStd(name='indentInteraction', createStepName='Initial', master=region1, slave=region2, sliding=FINITE, thickness=ON, interactionProperty='contact', adjustMethod=NONE, initialClearance=OMIT, datumAxis=None, contactTracking=ONE_CONFIG, clearanceRegion=None)
\end{lstlisting}
Tangential contact between the indenter and skin is modelled as frictionless (line \texttt{4}, the skin can freely slide over the indenter surface), and the normal contact behaviour is `hard' (line \texttt{5}, skin nodes are not allowed to penetrate the indenter surface). Furthermore, nodes of the skin are allowed to separate from the indenter (line \texttt{5}) when contact between the indenter and skin is 'released' (i.e.\ the skin nodes will not `stick' to the indenter).
......@@ -526,12 +541,13 @@ We finally use our earlier defined surface sets (lines \texttt{6--8}) to assign
Also in the \texttt{interaction} module, we need to convert the nodes on the tibia into a rigid surface. The relative positions of these tibia surface nodes will not change, modelling the (hollow) tibia bone as completely rigid (undeformable). Boundary conditions for these coupled nodes will be applied to the tibia's reference point.
We defined the tibia surface set and the reference point (set) earlier, so this is easy:
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# bone to rigid surface
a = myModel.rootAssembly
region4=a.instances['leg-1'].sets['bone']
region5=a.instances['leg-1'].sets['boneRP']
myModel.RigidBody(name='rigidBone', refPointRegion=region5, tieRegion=region4)
\end{lstlisting}
\subsection{The cast}
......@@ -541,7 +557,7 @@ When \texttt{param.cast = 1}, we also model the alginate cast around the leg. Th
\noindent The first task for \bas{buildFEM} when the cast is included in the model is to calculate the contours for the cast. For every (regularised) contour in the skin set, \bas{buildFEM} calculates the centre, and with this centre each of the 100 points in the skin contour are moved 0.1\,mm in the radial direction. Thus, on average there will be a 0.1\,mm gap between leg and cast at the start of the simulation.
These cast contours are written to a file \texttt{basil<basilid>\_castContours.py\index{basil<basilid>\_castContours.py@\texttt{basil<basilid>\_castContours.py}}} that is imported and lofted similar to the leg contours in a new \texttt{part}:
\begin{lstlisting}[language=apython,escapechar=!]
\begin{lstlisting}[language=apython, firstnumber=last, escapechar=!]
# build cast
myPart = myModel.Part(name='cast', dimensionality=THREE_D, type=DEFORMABLE_BODY)
# import contours
......@@ -561,7 +577,7 @@ myPart.ShellLoft(loftsections=((ledges[0], ), (ledges[1], ), (ledges[2], ), (led
Note that we do not close the cast to form a volume as we did with the leg (figure \ref{abaqus_castLofted}): the cast only limits radial displacement of the soft tissue, not the longitudinal displacement. Note also that the part is initially modelled as `deformable' (line \texttt{2}): we will convert the cast to a rigid body later on.\\
\noindent Next, we make a hole in the cast for the indenter. With the local leg orientation that was calculated earlier for partitioning of the skin, \bas{buildFEM} now calculates 3D coordinates for a `partitioning wire' that is projected onto the cast (similar as before with the assist-disc, figure \ref{abaqus_castWire}). The part of the cast inside the projected edge can then be removed:
\begin{lstlisting}[language=apython, escapechar=!]
\begin{lstlisting}[language=apython, firstnumber=last, escapechar=!]
# write partitioning wire
pts = [(-2.724500, -7.424954, -0.529141), (-3.159280, -7.252343, -1.342870), !\textnormal{(20 triplets of 3D coordinates)}! ]
myPart.WireSpline(points=pts, mergeType=IMPRINT, meshable=ON, smoothClosedSpline=ON)
......@@ -585,6 +601,7 @@ t = fa[0]
dface.append( t[0] );
# cut hole
myPart.RemoveFaces(faceList = dface, deleteCells=False)
\end{lstlisting}
Abaqus expects a `list' for \texttt{PartitionFaceByProjectingEdges}, so even if there is only one edge, we put it into a list (lines \texttt{10--13}). When the cast is partitioned with the wire (line \texttt{15}), we find the surface inside the contour and put it into a list (lines \texttt{17--21}). When we delete this partitioned surface (line \texttt{23}), there is a hole in the cast for the indenter.\\
......@@ -599,7 +616,7 @@ Abaqus expects a `list' for \texttt{PartitionFaceByProjectingEdges}, so even if
\noindent Contrary to the assist-disc, we cannot leave the partitioning wire `as is'. If we delete wire after partitioning, also the partition itself is deleted and the hole that we created will be closed again. If we leave it hovering in the model similar to the assist-disc, the simulation will not run. This is because the partitioning wire is modelled in the \texttt{part}, not as an individual \texttt{instance} as the assist-disc.
Therefore we connect the wire to the cast with a loft, creating a `ridge' in the proccess (figure \ref{abaqus_castRidge}):
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# create ridge
p = myModel.parts['cast']
e = p.edges
......@@ -608,11 +625,12 @@ ed = e.getClosest(coordinates=((-4.377070, -4.263596, 1.000000),))
t = ed[0]
# loft
p.ShellLoft(loftsections=((e[0], ), (t[0], )), startCondition=NONE, endCondition=NONE)
\end{lstlisting}
Note that the indenter can pass freely through the cast (ridge), as long as we do not define any interaction between the indenter and the cast (ridge). \\
\noindent With the geometry defined (figure \ref{abaqus_castLofted}), we can assign a surface and some sets to use later on:
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# surface and sets
f = myPart.faces
faces = f.getSequenceFromMask(mask=('[#3 ]', ), )
......@@ -624,11 +642,12 @@ myPart.Surface(side1Faces=side1Faces, side2Faces=side2Faces, name='cast')
myPart.ReferencePoint(point=(3, 3, 3))
lastRP = myPart.referencePoints.values()
myPart.Set(referencePoints=lastRP, name='castRP')
\end{lstlisting}
The (inside) surface of the cast needs to be defined for modelling contact with the skin later on; and because we are going to convert the cast to a rigid body, we need to define a reference point (lines \texttt{9-11}). The actual position of this reference point is irrelevant and \bas{buildFEM} puts it at \texttt{(3, 3, 3)} (line \texttt{9}).\\
\noindent The cast needs to be meshed before it can be converted to rigid surface because of the geometry (similar to the rigid surface of the tibia). So we assign a section to the cast, instance it in the model assembly, and mesh it with linear shell elements:
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# cast section assignment
myModel.SurfaceSection(name='cast', useDensity=OFF)
p = myModel.parts['cast']
......@@ -652,13 +671,18 @@ partInstances =(a.instances['cast-1'], )
a.seedPartInstance(regions=partInstances, size=0.500000, deviationFactor=0.1, minSizeFactor=0.1)
partInstances =(a.instances['cast-1'], )
a.generateMesh(regions=partInstances)
\end{lstlisting}
After the section assignment (lines \texttt{2--6}), the cast can be instanced `as is' (lines \texttt{8--10}).
\begin{figure}[t!]\center
\includegraphics[width=\linewidth]{abaqus_fullAssembly.png}\\
\caption{The full model assembly (including cast) in Abaqus (rotated for clarity). The yellow lines are the symmetry axis of the indenter and assist-disc (construction lines in the sketches), RP's are the reference points for the indenter, bone, and cast. The assist-disc will be dismissed when the simulation is submitted for analysis. \label{abaqus_fullAssembly}}
\end{figure}
The size for the seed of the cast mesh (line \texttt{21}) depends on the element order of the leg mesh. For linear elements in the leg, the seed size is equal to that of the leg: \texttt{param.seed}. For quadratic elements, \bas{buildFEM} uses a seed size that is half the seed size of the leg: \texttt{param.seed/2}. This ensures that the node density is similar between the (either linear or quadratic) leg mesh and (linear) cast mesh, which will give Abaqus less trouble with the contact interaction between leg and cast.\\
\noindent Once meshed, we can tie the nodes of the cast into a rigid surface, and define the contact interaction between leg and cast:
\begin{lstlisting}[language=apython]
\begin{lstlisting}[language=apython, firstnumber=last]
# cast to rigid surface
a = myModel.rootAssembly
region4=a.instances['cast-1'].sets['cast']
......@@ -671,21 +695,20 @@ a = myModel.rootAssembly
region1=a.instances['cast-1'].surfaces['cast']
region2=a.instances['leg-1'].surfaces['skin']
myModel.SurfaceToSurfaceContactStd(name='castInteraction', createStepName='Initial', master=region1, slave=region2, sliding=FINITE, thickness=ON, interactionProperty='contact', adjustMethod=NONE, initialClearance=OMIT, datumAxis=None, clearanceRegion=None)
\end{lstlisting}
The cast is converted into a rigid surface similar as the tibia surface (lines \texttt{2--6}). Contact between cast and leg (line \texttt{12}) is modelled with the same contact interaction as that between the indenter and leg.\\
\noindent If you run the main \texttt{basil<basilid>\-\_job\-.py\index{basil<basilid>\_job.py@\texttt{basil<basilid>\_job.py}}}-script up until this point, Abaqus should display something like figure \eqref{abaqus_fullAssembly}.
\begin{figure}[t!]\center
\includegraphics[width=\linewidth]{abaqus_fullAssembly.png}\\
\caption{The full model assembly (including cast) in Abaqus (rotated for clarity). The yellow lines are the symmetry axis of the indenter and assist-disc (construction lines in the sketches), RP's are the reference points for the indenter, bone, and cast. The assist-disc will be dismissed when the simulation is submitted for analysis. \label{abaqus_fullAssembly}}
\end{figure}
\subsection{Loading steps (essential boundary conditions)}
Abaqus creates a special initial step at the beginning of the model's step sequence and names it `Initial'. Abaqus creates only one initial step for your model, and it cannot be renamed, edited, replaced, copied, or deleted. This initial step allows you to define boundary conditions, predefined fields, and interactions that are applicable at the very beginning of the analysis. If a boundary condition or interaction is applied throughout the analysis, it is usually convenient to apply such conditions in the initial step. This 'initial' step is not an analysis step.
There are no prescribed pressures or forces in the model. We prescribe the movement of the tibia (rotation and translation) and the movement of the indenter (translation). When the cast is modelled we also prescribe that the cast is not allowed to move in any way (no rotations, no displacement, in Abaqus: `encastre'). The soft tissue of the leg is free to move about, except of course when movement happens to be limited by contact.
These essential boundary conditions are prescribed in 4 different steps. The first three steps \textsl{rotate} the tibia over three axis and have a duration of 1\,s each. In the fourth step we translate both the indenter and the tibia, simultaneously. This last step has a duration of \texttt{param.timePeriod} (100\,s by default):
\begin{lstlisting}[language=apython]
The (non-zero) translations and rotations are prescribed in four different analysis steps. The first three analysis steps rotate the tibia over three axis and have a duration of 1\,s each. In the fourth step we translate both the indenter and the tibia, simultaneously. This last step has a duration of \texttt{param.timePeriod} (100\,s by default). All boundary conditions in this model are applied to rigid bodies, so we assign the boundary conditions to the reference points:
\begin{lstlisting}[language=apython, firstnumber=last]
# create loading steps, BNC's
import step
Iregion=myModel.rootAssembly.instances['indenter-1'].sets['indenterRP']
......@@ -716,5 +739,13 @@ myModel.steps['step-1'].setValues(matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC
myModel.steps['step-2'].setValues(matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC, solutionTechnique=FULL_NEWTON)
myModel.steps['step-3'].setValues(matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC, solutionTechnique=FULL_NEWTON)
myModel.steps['step-4'].setValues(matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC, solutionTechnique=FULL_NEWTON)
\end{lstlisting}
The first step is used to rotate the tibia over the $z$-axis (line \texttt{6}). Translation of the tibia is prohibited (line \texttt{8})
\ No newline at end of file
In the initial step, we fix the rotations of the indenter (line \texttt{6}) and any movement of the cast (if applicable, lines \texttt{8--9}).
In the first analysis step (lines \texttt{10--14}), we first prohibit translations of the indenter (line \texttt{12}) and tibia (line \texttt{13}). These zero values for translations will be propagated throughout the following steps, until we will change them in the fourth (last) analysis step. In this first step, we only rotate the tibia over the $z$-axis (line \texttt{14}). Note that the tibia rotations are prescribed as rotation \textsl{velocities}. So when the duration of the step is changed, these velocities will also have to change.
In the second analysis step, rotation velocity over the $z$-axis is (re)set to zero, and we rotate the tibia over the $y$-axis (line \texttt{17}). Similarly, in the third analysis step the rotation velocity over the $y$-axis is (re)set to zero, and we rotate the tibia over the $x$-axis (line \texttt{20}).
The fourth analysis step is used to translate the tibia (which now has the `deformed' orientation) and to actually indent the leg. First we (re)set the tibia rotation over the $x$-axis to zero (line \texttt{23}), and next we prescribe the tibia (line \texttt{24}) and indenter (line \texttt{25}). Note that \texttt{param}-parameters for \bas{buildFEM} for the step duration and -incrementation only affect this last (indentation) step.
......@@ -131,7 +131,7 @@ morekeywords = {WireSpline,
maxInc,
nlgeom,
DisplacementBC,
u1, u2, u3, ur1, ur2, ur3, vr1, vr2, vr3
u1, u2, u3, ur1, ur2, ur3, vr1, vr2, vr3,
amplitude,
fixed,
distributionType,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment