Differences

This shows you the differences between two versions of the page.

--- commit:futur:db [2018/07/12 12:15] – [How to use the new feature] boemer
+++ commit:futur:db [2018/07/17 12:03] (current) – boemer
@@ Line 1: / Line 1: @@
-====== Commit 2018-07-12 ======
+====== Commit 2018-07-17 ======
-In this commit, the traction element was created to set the traction force to zero at the nodes that are in contact within a specified contact interaction.  This feature is required in the coupling procedure with Metalub (see commit 3186 of Metafor and commit 1456 of Metalub for more details).
+===== Description =====
+In this commit, the out-of-plane thickness evolution in the generalized plane strain state is coherently accounted for in the computation.  Previously, the possible increase or decrease of the out-of-plane thickness via the ''depend'' command was not taken into account in the contact area computation.  This led to an incoherently changing gap in the following scenario, where ''AIC_ONCEPERSTEP'' is activated.
-===== How to use the new feature =====
-A ''ContactTraction2DElement'' or ''ContactTraction3DElement'' is created like a ''Traction2DElement'' or ''Traction3DElement''.  The only additional property that has to be set is ''NBR_CONTACT_INTERACTION'', which is the number of the contact interaction.  For instance,
-  prp = ElementProperties(ContactTraction2DElement)
+In the previous figure, the square is represented in 3D although the computation is performed in the generalized plane strain state.  The lateral movement of the left and right sides of the square is blocked.  A pressure $p$, identical in both cases, pushes the square against the contact tool (frictionless).  When the out-of-plane thickness of the square increases, as shown on the right, the resultat force created by the pressure increases due to the increase of the surface on which this pressure is applied.  Moreover, if the increase of the out-of-plane thickness is not taken into account in the contact interaction, i.e. in the computation of the area in contact, the increase of the resultat pressure force will increase the gap.  Hence,  $\text{gap}_{1} < \text{gap}_{2}$  previously.  In the new version,  $\text{gap}_{1} = \text{gap}_{2}$ .  In fact, the gap should remain constant, if the contact pressure remains constant.
-  prp.put(PRESSURE, -p['pMax'])
-  prp.put(NBR_CONTACT_INTERACTION, 3)
-Hence, if nodes in this ''ContactTraction2DElement'' are slave nodes in contact in the contact interaction 3, the traction force applied to these nodes is set to zero.
-===== File structure =====
+===== Test case =====
-Instead of implementing the new feature directly in the ''TractionElement'' of the ''mtElement'' package, a new template class ''ContactTractionElement'' was created in the ''mtContact'' (elements) package.  It had to be moved to the ''mtContact'' package (instead of the ''mtElement'' package) in order to have access to the contact features via ''include "ContactInteraction.h"''.  Moreover, the ''ContactTractionElement'' derives of ''TractionElement'' to minimize duplicate code.  In a similar way to the ''TractionElement'', 2D and 3D elements were created: ''ContactTraction2DElement'' and ''ContactTraction3DElement''.  This required to export all methods of ''Traction2DElement'' and ''Traction3DElement'' via ''MTELEMENTS_API'' (+ ''ELEMENT_CPP_MACRO2'' extension), which only had to be added in ''Traction3DElement'' since it was already done for ''Traction2DElement'', certainly to derive it in ''mtXFEM''.
+The test case in the following figure is identical to what was decribed previously.  A square is pushed against a contact tool (frictionless) with a pressure  $p$ , its lateral movement is blocked and its out-of-plane thickness increases (as a function of time).
-Some minor modifications:
-  * The new property ''NBR_CONTACT_INTERACTION'' was added in ''mtElements.h/cpp'' as an integer property.
-  * It was necessary to add ''ContactTraction2DElement'' and ''ContactTraction3DElement'' to the functions ''is2D'' and ''is3D'' in ''LoadingInteraction.cpp'' in order to be able to add them to the ''LoadInteraction'' in the Python data set.
-  * Untabify in ''ContactElement.h/cpp''.
-  * ''getSlaveNode()'', previously private method in ''ContactElement.h/cpp'', becomes protected to use it in ''ContactTractionElement.inl''.
-  * The method ''isInContactPrev()'' was added in ''ContactElement.h/cpp'' to get the contact status at the end of the previous time iteration.  Using this contact status has the advantage that the contact traction activation/deactivation does not happen during the same time step, i.e. during the mechanical iteration.  It might be possible to use the current contact status but this could introduce some undesired oscillations.
+The gap can be computed in various ways, numerically, analytically and semi-analytically/numerically.  This is what is done in three different extractors to test the good implementation of the previous feature:
+  * ''numericalGap'': First, the numerical computation simply consists in extracting the gaps in the Metafor computation and choosing the one at  $P_{4}$ .
+  * ''analyticalGap'': Second, the analytical computation provides the value of the gap via the following formula:  $\text{gap} = p/k_{n}$ , where  $k_{n}$  is the normal penalty coefficient.
+  * ''semiNumAnalyGap'': Thirdly, the semi-analytical/numerical computation consists in dividing the nodal contact force along  $\mathbf{y}$ ,  $f_{y}$ , at  $P_{4}$  by the nodal area  $l_{x}l_{z}/n_{x}$  (where  $l_{z} = \epsilon_{z}$ , since the initial out-of-plane thickness is 1, and where  $n_{x}$  is the number of elements along  $\mathbf{x}$ ) and the normal penalty  $k_{n}$ :  $\text{gap} = f_{y}/( k_{n}l_{z}l_{x}/(2n_{x}) )$
+In all three cases, the gap should be the same.
-===== Implementation =====
-Concerning the implementation, the ''ContactTractionElement'' is very similar to the ''TractionElement'' except from parts of the following virtual functions:
-  * ''checkProps'': This function checks whether the new property ''NBR_CONTACT_INTERACTION'' was set during the initialization of the elements.  If this property was set, it checks whether the specified number corresponds to an existing contact interaction.  Once this interaction is found, pointers to the slave nodes corresponding to the nodes of the ''ContactionTractionElement'' are saved in the new ''ces'' private vector of contact elements in the ''ContactTractionElement''.
-  * ''fillExternalForces'': This function uses the ''fillExternalForces'' function of ''TractionElement'' first and then it sets all contact forces corresponding to nodes in contact equal to zero.  As mentionned earlier, the contact state is determined based on the results of the previous time integration for reasons of possible oscillation minimization.
-  * ''noAnalyticalStiffnessAvailable'':  An analytical stiffness matrix was defined for the ''ContactTractionElement''.  Hence, this function can return ''false''.
-  * ''fillAnalyticalStiffness'': Similar to the ''fillExternalForces'' function, the ''fillAnalyticalStiffness'' of ''TractionElement'' is used first and then the rows corresponding to the traction forces that were set equal to zero are also set equal to zero.  The structure of this matrix is illustrated in a comment of the ''fillAnalyticalStiffness'' function to understand the details.
-===== Test cases =====
-Two test cases were created in the folder ''mtContact/tests'' to test the ''ContactTractionElement'' in 2D and 3D.  Hence, they have the same names than the elements that are tested: ''contactTraction2DElement.py'' and ''contactTraction3DElement''.  The 3D test is the natural extension of the 2D case.  Hence, only some details about the 2D case are explained here.
-As shown in the previous figure, two squares are crushed by a compression pressure and a contact tool that moves downwards by  $2 \; dy$ .  In the case on the left, the traditional ''Traction2DElement'' is used, while the new ''ContactTraction2DElement'' is used on the right.  The ''ContactTraction2DElement'' works correctly, if the total force on node  $P_{14}$  along  $\mathbf{y}$  is equal to the contact force along  $\mathbf{y}$  at the end of the computation, i.e. that the traction force has been set to zero for the node in contact.  To check this equality, the difference between both forces (total and contact) is computation in the extractor ''diffTotalAndContactForceInCase2''.  The results of the cases with the ''Traction2DElement'' and the ''ContactTraction2DElement'' are also compared by computing the differences of the total force along  $\mathbf{y}$  at the nodes  $P_{4}$  and  $P_{14}$  in the extractor ''diffTotalForcesBetween2Cases''.  This difference should be different from 0, unlike the previous one.