​

En este trabajo se propone un algoritmo mortar para el estudio del contacto mecánico en problemas de elasticidad tridimensional. La superficie de proyección utilizada para la integración de las ecuaciones se define a través de una base local cartesiana en cada elemento de contacto, de esta forma se facilita la implementación del algoritmo y la linealización de las ecuaciones. Se proponen ejemplos donde se demuestra que el algoritmo satisface los test de la parcela. Finalmente se utiliza el algoritmo en un caso de interés industrial, el contacto de una válvula de motor de combustión interna con su asiento.

​

==Abstract==

​

In this paper we propose a mortar algorithm for the study of contact mechanics in three dimensional elasticity problems. The projection surface used for integrating the equations is selected through a local cartesiana base defined in each contact element. In this way, some difficulties in the algorithm implementation as well as in the linearization of the equations are avoided. The proposed examples show that the algorithm satisfies the patch tests. Finally, we use the algorithm in an industrial application, the contact of an internal combustion engine valve with its seat.

​

==Palabras clave==

​

Contacto ; Lagrangiano aumentado ; Algoritmo mortar ; Elementos finitos

​

==Keywords==

​

Contact ; Augmented Lagrangian ; Mortar ; Finite elements

​

==1. Introducción==

​

La simulación de componentes mecánicos y procesos que involucran el contacto entre cuerpos tiene gran aplicación en diferentes áreas de la ingeniería, pudiendo citarse el diseño de engranajes [[#bib0005|[1]]] , procesos de embutición [[#bib0010|[2]]] , fatiga por contacto [[#bib0015|[3]]] , entre otras. Para la descripción del desplazamiento relativo entre 2 cuerpos en contacto mediante el Método de los Elementos Finitos (MEF), una técnica ampliamente difundida es la del nodo-segmento, en la cual a un nodo de un cuerpo, denominado ''esclavo'' , se asocia una zona de un segmento o superficie del otro cuerpo denominado ''master''[[#bib0020|[4]]] . El principal inconveniente de la estrategia nodo-segmento es que no garantiza que una superficie transmita una presión de contacto uniforme a la otra superficie; en otras palabras, esta estrategia no supera los denominados tests de la parcela de contacto [[#bib0025|[5]]] . Las aproximaciones del tipo nodo-segmento con doble pasada verifican los tests de la parcela, pero pueden bloquearse debido a las sobre restricciones introducidas en la formulación; adicionalmente, con superficies de contacto no suaves, generan saltos de la presión cuando los nodos esclavos se deslizan entre segmentos adyacentes [[#bib0030|[6]]]  y no cumplen con la condición de Babuska-Brezzi [[#bib0035|[7]]]  causando un mal condicionamiento de las matrices y pobre tasa de convergencia, típicas manifestaciones de la sobre-restricción. A pesar de ello, estos métodos son a menudo elegidos en problemas bidimensionales y en ciertas configuraciones de mallas tridimensionales con elementos de bajo orden, para analizar problemas de contacto entre cuerpos flexibles por su capacidad de verificar el test de la parcela. Sin embargo, no pueden superar los tests de la parcela para elementos de orden superior  [[#bib0030|[6]]]  and [[#bib0040|[8]]] .      

​

Otra aproximación es la de los métodos segmento-segmento que relacionan el segmento de un cuerpo (segmento ''master'' ) con algún segmento del otro cuerpo (segmento ''esclavo'' ). La mayoría de los métodos segmento-segmento utilizan algún tipo de superficie de contacto intermedia, o proyección de superficies e integran el trabajo virtual de contacto usando una cuadratura numérica con alguna interpolación en la presión. Estos métodos fueron aplicados inicialmente a interfases de contacto planas para ejemplos 2D y luego ampliados a mallas 3D  [[#bib0045|[9]]] , [[#bib0050|[10]]]  and [[#bib0055|[11]]] . En el trabajo de Park y Felippa [[#bib0060|[12]]]  se propone una formulación con una superficie de contacto intermedia entre los cuerpos de contacto determinada por un conjunto de ecuaciones auxiliares tal que se satisfaga el equilibrio en dicha superficie. La ventaja de este método es que evita la integración de los campos de multiplicadores de Lagrange, pero tiene como inconveniente la complejidad para ser extendido a problemas 3D. Chen y Hisada [[#bib0065|[13]]]  desarrollaron una modificación a un esquema del tipo nodo segmento el cual es capaz de superar los test de la parcela de contacto. En este caso, la contribución del contacto al trabajo virtual se integra en las superficies de contacto a través de la presión de contacto nodal. Como desventaja, la matriz de rigidez tangente no es simétrica. Kim et al. desarrollaron una nueva estrategia para problemas de contacto bidimensionales dentro del rango de las deformaciones infinitesimales [[#bib0070|[14]]] . Este algoritmo utiliza elementos finitos de variables nodales. La idea es transformar los problemas de contacto nodo-superficie en problemas de contacto nodo-nodo. Para el caso de realizar una extensión a problemas tridimensionales, es necesario proponer un algoritmo robusto de búsqueda de segmentos como así también desarrollar elementos poligonales de variables nodales para que pueda ser tratado en cualquier superficie arbitraria. La formulación es capaz de superar los test de la parcela, sin embargo, no se han reportado hasta el momento extensiones a geometrías tridimensionales y en el rango de grandes deformaciones. Más recientemente, Zavarise y Lorenzis [[#bib0075|[15]]] , propusieron un algoritmo nodo segmento con una regularización del tipo penalidad capaz de superar los test de la parcela. Sin embargo, al igual que en el caso anterior, no se reportaron casos de problemas tridimensionales.      

​

Las diferentes aproximaciones propuestas de métodos segmento-segmento difieren generalmente en la forma de interpolar la presión de contacto y en la manera en que se efectúa la proyección entre los elementos. Los primeros trabajos desarrollados utilizaron el método de penalidad como regularización del problema variacional y con aplicaciones solo a problemas bidimensionales [[#bib0080|[16]]] , [[#bib0085|[17]]]  and [[#bib0025|[5]]] . El método superficie-superficie tipo mortar fue originalmente propuesto como un método de descomposición de dominios y utilizado para resolver problemas de elementos finitos con discretizaciones no conformes [[#bib0045|[9]]]  and [[#bib0090|[18]]] . Específicamente, el primer trabajo que ha utilizado el método mortar fue el propuesto por Bernardi et al. [[#bib0095|[19]]] , donde se muestra la estabilidad del método en relación a las condiciones de Babuska-Brezzi. En el trabajo de Puso y Laursen [[#bib0100|[20]]] , se presenta una versión del método mortar aplicado al contacto mecánico y utilizado en problemas 3D con grandes deformaciones. La proyección de las superficies de los elementos en contacto se realiza sobre una superficie plana; en cambio, en la propuesta de Puso [[#bib0105|[21]]] , la proyección se lleva a cabo sobre una malla intermedia suave definida por un parche con funciones de forma de Hermite que garantiza continuidad ''C''<sup>1</sup>  solo cuando la malla es regular. El método mortar ha continuado su desarrollo extendiéndolo a aplicaciones con fricción o en problemas de auto-contacto  [[#bib0110|[22]]] , [[#bib0115|[23]]]  and [[#bib0120|[24]]] .      

​

En este trabajo, siguiendo algunas de las ideas presentadas por Puso y Laursen [[#bib0100|[20]]] , se propone un algoritmo de contacto tipo mortar con un esquema de penalidad aplicable a problemas tridimensionales en donde, a diferencia de la propuesta de Puso y Laursen [[#bib0100|[20]]] , la superficie topológica de proyección de los elementos se caracteriza mediante una base local cartesiana definida para cada elemento no-mortar en contacto. Esta característica simplifica el tratamiento algebraico de las ecuaciones y la implementación computacional. El algoritmo fue implementado en el código de elementos finitos Oofelie [[#bib0125|[25]]] . Los ejemplos numéricos propuestos muestran que el algoritmo supera los tests de la parcela con resultados de tensiones muy regulares para mallas arbitrarias. En todos los casos las simulaciones fueron llevadas a cabo asumiendo un modelo mecánico elástico lineal, sin consideración de efectos dinámicos ni procesos de fricción entre los cuerpos en contacto.      

​

==2. Descripción del problema==

​

El movimiento de los cuerpos <math display="inline">B</math><sup>''α''</sup> , con ''α''  = 1, 2 se expresa por el mapeo <math display="inline">{\chi }^{\alpha }:{\Omega }^{\alpha }\times [0,T]\longrightarrow {\mathbb{R}}^3</math> , y la posición actual de las partículas materiales se calcula por '''x'''<sup>''α''</sup>  = ''χ''<sup>''α''</sup> (''X''<sup>''α''</sup> , ''t'' ). A diferencia de los métodos nodo-segmento, en el contacto tipo mortar la porción de frontera en contacto ''Γ''<sub>''c''</sub>  está dada por una intersección de las superficies <math display="inline">{\Gamma }_c^{\alpha }</math> , esto es <math display="inline">{\Gamma }_c={\Gamma }_c^1\cap {\Gamma }_c^2</math> , ver [[#fig0005|figura 1]] . La convención adoptada en este trabajo para la definición de las superficies mortar y no mortar es la siguiente: <math display="inline">{\Gamma }_c^1</math>  denota la superficie no-mortar o esclava, en tanto que <math display="inline">{\Gamma }_c^2</math>  es la superficie mortar o ''master'' . La estrategia utilizada para discretizar el dominio es por medio del MEF, de esta forma, la malla de la superficie de contacto puede ser parametrizada como

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{x}^{\alpha }=\sum_{A=1}^{n^{\alpha }}N_A^{\alpha }({\boldsymbol{\xi }}^{\alpha })\boldsymbol{x}_A^{\alpha },</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)

|}

​

donde: <math display="inline">\boldsymbol{x}^{\alpha }\in {\Gamma }_c^{\alpha }\rightarrow {\mathbb{R}}^3</math> , <math display="inline">\boldsymbol{x}_A^{\alpha }\in {\Gamma }_c^{\alpha }\rightarrow {\mathbb{R}}^3</math> , ''n''<sup>''α''</sup>  es el número de nodos de la malla en <math display="inline">{\Gamma }_c^{\alpha }</math>  y <math display="inline">N_A^{\alpha }:{\Gamma }_c^{\alpha }\rightarrow \mathbb{R}</math>  son las clásicas funciones de forma, bilineales en caso de utilizar mallas con hexaedros, o lineales para el caso de mallas formadas por tetraedros.

​

<span id='fig0005'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_843322741-1-s2.0-S0213131512000156-gr1.jpg|center|488px|Notación utilizada en el método mortar para los cuerpos en contacto.]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 1.

​

Notación utilizada en el método mortar para los cuerpos en contacto.

​

</span>

|}

​

El vector de coordenadas nodales en <math display="inline">{\mathbb{R}}^3</math>  para las superficies de contacto se define de la siguiente manera

​

<span id='eq0010'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\Phi }=\left[\begin{array}{c}

\boldsymbol{x}_1^1\\

\boldsymbol{x}_2^1\\

\vdots \\

\boldsymbol{x}_n^1\\

\boldsymbol{x}_1^2\\

\boldsymbol{x}_2^2\\

\vdots \\

\boldsymbol{x}_n^2\\

​

\end{array}\right].</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)

|}

​

Notar que '''Φ'''  se ordena de forma tal que en las primeras componentes se encuentran las coordenadas de los nodos referidos a la frontera <math display="inline">{\Gamma }_c^1</math>  (no mortar) y luego las correspondientes a <math display="inline">{\Gamma }_c^2</math>  (mortar).      

​

==3. Trabajo virtual asociado al contacto==

​

Cuando los cuerpos <math display="inline">B</math><sup>1</sup>  y <math display="inline">B</math><sup>2</sup>  interactúan mecánicamente sobre la superficie de contacto, el trabajo virtual asociado al contacto queda definido por

​

<span id='eq0015'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\delta {\Pi }^c=\sum_{\alpha =1}^2{\int }_{{\Gamma }_c^{\alpha }}\boldsymbol{t}^{\alpha }\cdot \delta \boldsymbol{x}^{\alpha }\quad d{\Gamma }_c^{\alpha }=</math><math>\delta \boldsymbol{\Phi }\cdot \boldsymbol{F},</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)

|}

​

donde '''''t'''''<sup>''α''</sup>  denota el vector tracción de Cauchy. Asumiendo que existe conservación de la cantidad de movimiento lineal en la interfase de contacto, <math display="inline">\boldsymbol{t}^1\quad d{\Gamma }_c^1=</math><math>-\boldsymbol{t}^2\quad d{\Gamma }_c^2</math> ; y que '''''x'''''<sup>1</sup>  − '''''x'''''<sup>2</sup>  ≈ 0, la ecuación [[#eq0015|(3)]]  puede ser escrita como

​

<span id='eq0020'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\delta {\Pi }^c={\int }_{{\Gamma }_c^1}\boldsymbol{t}^1\cdot (\delta \boldsymbol{x}^1-</math><math>\delta \boldsymbol{x}^2)\quad d{\Gamma }_c^1=\delta \boldsymbol{\Phi }\cdot F.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)

|}

​

Utilizando las siguientes parametrizaciones para las superficies de contacto y el vector tracción

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{x}^1=\sum_{B=1}^{n^1}N_B^1({\boldsymbol{\xi }}^1)\boldsymbol{x}_B^1,</math>

|-

|<math>\boldsymbol{x}^2=\sum_{C=1}^{n^2}N_C^2({\boldsymbol{\xi }}^2)\boldsymbol{x}_C^2,</math>

|-

|<math>\boldsymbol{t}=\sum_{A=1}^{n^1}N_A^1({\boldsymbol{\xi }}^1)\boldsymbol{t}_A,</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)

|}

​

junto con la ecuación [[#eq0020|(4)]] , el principio de los trabajos virtuales asociado al contacto se aproxima de la siguiente manera

​

<span id='eq0030'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\delta {\Pi }^c=\sum_{A=1}^{n^1}\sum_{B=1}^{n^1}\boldsymbol{t}_A\cdot {\int }_{{\Gamma }_c^1}N_A^1({\boldsymbol{\xi }}^1)N_B^1({\boldsymbol{\xi }}^1)\quad d{\Gamma }_c^1\quad \delta \boldsymbol{x}_B^1-</math><math></math>

|-

|<math>\sum_{A=1}^{n^1}\sum_{C=1}^{n^2}\boldsymbol{t}_A\cdot {\int }_{{\Gamma }_c^1}N_A^1({\boldsymbol{\xi }}^1)N_C^2({\boldsymbol{\xi }}^2)\quad d{\Gamma }_c^1\quad \delta \boldsymbol{x}_C^2.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)

|}

​

Luego, la ecuación [[#eq0030|(6)]]  puede ser escrita en forma más compacta

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\delta {\Pi }^c=\sum_{A=1}^{n^1}\boldsymbol{t}_A\cdot \delta \boldsymbol{g}_A,</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)

|}

​

donde

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\delta \boldsymbol{g}_A=\sum_{B=1}^{n^1}n_{AB}^1\quad \delta \boldsymbol{x}_B^1-</math><math>\sum_{C=1}^{n^2}n_{AC}^2\quad \delta \boldsymbol{x}_C^2,</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)

|}

​

denota la primera variación del huelgo o ''gap''  del nodo esclavo ''A'' , y

​

<span id='eq0045'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>n_{AB}^1={\int }_{{\Gamma }_c^1}N_A^1({\boldsymbol{\xi }}^1)N_B^1({\boldsymbol{\xi }}^1)\quad d{\Gamma }_c^1,</math>

|-

|<math>n_{AC}^2={\int }_{{\Gamma }_c^1}N_A^1({\boldsymbol{\xi }}^1)N_C^2({\boldsymbol{\xi }}^2)\quad d{\Gamma }_c^1,</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)

|}

​

son denominados factores de peso.      

​

Desarrollando las sumatorias y teniendo en cuenta la notación definida para el vector de coordenadas nodales '''Φ'''  dada por la ecuación [[#eq0010|(2)]] , la variación de la energía potencial elástica debida al contacto en el caso de que los elementos de la malla sean hexaedros lineales tiene la siguiente expresión matricial

​

<span id='eq0050'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\delta {\Pi }^c=\left[\begin{array}{c}

\delta \boldsymbol{x}_1^1\\

\delta \boldsymbol{x}_2^1\\

\delta \boldsymbol{x}_3^1\\

\delta \boldsymbol{x}_4^1\\

\delta \boldsymbol{x}_1^2\\

\delta \boldsymbol{x}_2^2\\

\delta \boldsymbol{x}_3^2\\

\delta \boldsymbol{x}_4^2\\

​

\end{array}\right]\cdot \left[\begin{array}{c}

\sum_{A=1}^{n^1}\boldsymbol{t}_A\quad n_{A1}^1\\

\sum_{A=1}^{n^1}\boldsymbol{t}_A\quad n_{A2}^1\\

\sum_{A=1}^{n^1}\boldsymbol{t}_A\quad n_{A3}^1\\

\sum_{A=1}^{n^1}\boldsymbol{t}_A\quad n_{A4}^1\\

\sum_{A=1}^{n^1}-\boldsymbol{t}_A\quad n_{A1}^2\\

\sum_{A=1}^{n^1}-\boldsymbol{t}_A\quad n_{A2}^2\\

\sum_{A=1}^{n^1}-\boldsymbol{t}_A\quad n_{A3}^2\\

\sum_{A=1}^{n^1}-\boldsymbol{t}_A\quad n_{A4}^2\\

​

\end{array}\right],</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)

|}

​

pudiendo extenderse de manera similar a otras topologías de elementos. De la ecuación [[#eq0050|(10)]]  se desprende que el vector de fuerzas internas del contacto para el caso de mallas formadas por elementos hexaedros es:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{F}=\left[\begin{array}{c}

\boldsymbol{F}_1^1\\

\boldsymbol{F}_2^1\\

\boldsymbol{F}_3^1\\

\boldsymbol{F}_4^1\\

-\boldsymbol{F}_1^2\\

-\boldsymbol{F}_2^2\\

-\boldsymbol{F}_3^2\\

-\boldsymbol{F}_4^2\\

​

\end{array}\right]</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)

|}

​

con

​

<span id='eq0060'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{F}_B^{\alpha }=\sum_{A=1}^{n^1}\boldsymbol{t}_A\quad n_{AB}^{\alpha }.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)

|}

​

Los desarrollos que se presentaron hasta aquí no tienen en cuenta efectos dinámicos, tales como el impacto. Para su consideración es necesario realizar los cálculos detallados a cada paso de tiempo. Se requiere utilizar un algoritmo de integración temporal apropiado, resultando aconsejable el empleo de un algoritmo de integración que garantice la conservación o disipación de energía [[#bib0130|[26]]] , [[#bib0135|[27]]] , [[#bib0140|[28]]]  and [[#bib0145|[29]]] . La formulación de contacto mostrada puede adaptarse fácilmente a un esquema de este tipo.      

​

==4. Regularización de las restricciones normales==

​

El requerimiento físico de impenetrabilidad e interacción de compresión entre los cuerpos es impuesto por medio de una restricción de contacto normal. El vector tracción de Cauchy '''''t'''''  puede ser divido en una componente normal '''''t'''''<sub>''N''</sub>  y en otra tangencial '''''t'''''<sub>''T''</sub> . Luego, con las expresiones discretas del trabajo virtual, la restricción del contacto normal puede establecerse con las condiciones de Karush-Kuhn-Tucker (KKT) como

​

<span id='eq0065'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>t_{N_A}\geq 0,\quad \quad g_A\leq 0,\quad \quad t_{N_A}g_A=</math><math>0,</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)

|}

​

expresadas aquí en forma discreta para cada nodo ''A''  no-mortar. Las ecuaciones  [[#eq0065|(13)]]  representan un proceso de contacto sin fricción, ya que no se ha tenido en cuenta la componente tangencial '''''t'''''<sub>''T''</sub> . La extensión del algoritmo mortar que aquí se presenta a problemas con fricción mediante un ''return mapping''  no resulta dificultosa  [[#bib0115|[23]]] , ya que solo se debe incluir la matriz de fricción que surge de considerar a '''''t'''''<sub>''T''</sub>  como función del desplazamiento tangencial. De esta forma, no hay necesidad de incluir variables adicionales.      

​

El huelgo ''g''<sub>''A''</sub>  de un nodo ''A''  no-mortar es: ''g''<sub>''A''</sub>  = '''''g'''''<sub>''A''</sub>  · '''''ν'''''<sub>''A''</sub> , donde '''''ν'''''<sub>''A''</sub>  es el vector normal promedio en un nodo no-mortar ''A''  tal como se muestra en la  [[#fig0010|figura 2]] . El promediado de las normales concurrentes a un nodo es una operación relativamente sencilla de implementar y aporta resultados satisfactorios como lo demuestran a continuación los ejemplos numéricos propuestos. Para asegurar el cumplimiento de las condiciones de KKT dadas en la ecuación [[#eq0065|(13)]] , se ha seleccionado una regularización por medio del método de penalidad. De esta forma, la tracción normal en el nodo ''A''  tiene la siguiente expresión

​

<span id='eq0070'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>t_{N_A}={\kappa }_A\quad \boldsymbol{g}_A\cdot {\boldsymbol{\nu }}_A,\quad \quad (sin\quad suma\quad en\quad \mbox{ A})</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)

|}

​

donde

​

<span id='eq0075'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{g}_A=\sum_{B=1}^{n^1}n_{AB}^1\quad \boldsymbol{x}_B^1-</math><math>\sum_{C=1}^{n^2}n_{AC}^2\quad \boldsymbol{x}_C^2,</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)

|}

​

y ''κ''<sub>''A''</sub>  es un factor definido como

​

<span id='eq0080'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\kappa }_A=\frac{\varepsilon }{{\sum }_Dn_{AD}^1},</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16)

|}

​

siendo ''ɛ''  un coeficiente de penalidad propuesto por el usuario. En la ecuación  [[#eq0080|(16)]] , el coeficiente de penalidad ''ɛ''  se debe ajustar a valores típicos de rigidez de los cuerpos en contacto (un valor apropiado está en el orden de 10''Eh'' , siendo ''E''  el módulo de elasticidad y ''h''  el tamaño de malla). Esta característica es fundamental para lograr una buena tasa de convergencia. Luego, con las expresiones de las ecuaciones ( [[#eq0060|12]] ,[[#eq0070|14]] ), el vector de fuerzas internas resulta en

​

<span id='eq0085'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{F}_B^{\alpha }=\sum_{A=1}^{n^1}\quad {\kappa }_A\quad n_{AB}^{\alpha }\left({\boldsymbol{\nu }}_A\otimes {\boldsymbol{\nu }}_A\right)\boldsymbol{g}_A.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 17)

|}

​

<span id='fig0010'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_843322741-1-s2.0-S0213131512000156-gr2.jpg|center|281px|Vector normal promedio νA definido por los elementos: (J−1), (J) y (J+1) en el ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 2.

​

Vector normal promedio '''''ν'''''<sub>''A''</sub>  definido por los elementos: (''J''  − 1), (''J'' ) y (''J''  + 1) en el nodo ''A'' .                  

​

</span>

|}

​

==5. Cómputo de los factores de peso. Proyecciones==

​

Para llevar a cabo las integraciones en la interface ''Γ''<sub>''c''</sub> , es necesario definir una superficie topológica sobre la cual realizar las aproximaciones. En este trabajo se aproxima la cara de un elemento ''k''  no mortar, ver  [[#fig0005|figura 1]] , con un elemento <math display="inline">\tilde{k}</math>  que pertenece a un plano '''''p'''''  definido por un vector normal '''''e'''''<sub>3</sub> , ver [[#fig0015|figura 3]] -a. Dicho vector normal se obtiene por el producto vectorial de las diagonales de la cara del elemento ''k  '' . La suma de todas las superficies <math display="inline">\tilde{k}</math>  representa la superficie de integración ''Γ''<sub>''c''</sub>  de la ecuación [[#eq0015|(3)]] . Una alternativa para la proyección de los nodos no-mortar y mortar en el plano '''''p'''''  es a través de una base global definida en <math display="inline">{\mathbb{R}}^3</math> . De esta manera, es posible utilizar el algoritmo de intersección lineal paramétrico Cyrus-Beck descrito en la referencia [[#bib0150|[30]]] . En este trabajo se propone otra alternativa que facilita el desarrollo de las ecuaciones y que consiste en la definición de una base local para cada elemento no mortar que define un plano o superficie de proyección '''''p''''' , tal como lo representa la [[#fig0015|figura 3]] -a. De esta forma, es posible utilizar un algoritmo de intersección de polígonos basado en reglas de avance y en la definición paramétrica de rectas en <math display="inline">{\mathbb{R}}^2</math> , ver por ejemplo [[#bib0160|[31]]] . Tomando una base <math display="inline">V=\lbrace {\overset{\textasciicaron}{\boldsymbol{e}}}_1,{\overset{\textasciicaron}{\boldsymbol{e}}}_2,{\overset{\textasciicaron}{\boldsymbol{e}}}_3\rbrace </math> , donde los versores <math display="inline">{\overset{\textasciicaron}{\boldsymbol{e}}}_1</math>  y <math display="inline">{\overset{\textasciicaron}{\boldsymbol{e}}}_3</math>  son calculados por medio de las diagonales para cada elemento no-mortar y el vector <math display="inline">{\overset{\textasciicaron}{\boldsymbol{e}}}_2</math> , como <math display="inline">{\overset{\textasciicaron}{\boldsymbol{e}}}_1\times {\overset{\textasciicaron}{\boldsymbol{e}}}_3</math> , ver [[#fig0015|figura 3]] -a, las proyecciones de las coordenadas de los nodos <math display="inline">\boldsymbol{x}_A^{\alpha }</math>  de un elemento mortar o no-mortar en '''''p'''''  pueden calcularse como

​

<span id='eq0090'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{y}_A^{\alpha }=\left[\begin{array}{c}

{\overset{\textasciicaron}{\boldsymbol{e}}}_1\cdot (\boldsymbol{x}_A^{\alpha }-\boldsymbol{x}_0^1)\\

{\overset{\textasciicaron}{\boldsymbol{e}}}_2\cdot (\boldsymbol{x}_A^{\alpha }-\boldsymbol{x}_0^1)\\

​

\end{array}\right],</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)

|}

​

donde <math display="inline">\boldsymbol{x}_0^1</math>  es un punto que pertenece a un elemento no-mortar ''k  '' , por ejemplo, uno de sus vértices. Así, las coordenadas <math display="inline">\boldsymbol{y}_A^{\alpha }</math>  de los nodos no-mortar y mortar quedan definidas con solo 2 componentes.

​

<span id='fig0015'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_843322741-1-s2.0-S0213131512000156-gr3.jpg|center|491px|(a) Elemento no-mortar k, elemento mortar l, plano de proyección p, base ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 3.

​

(a) Elemento no-mortar ''k'' , elemento mortar ''l'' , plano de proyección '''''p''''' , base <math display="inline">V=\lbrace {\overset{\textasciicaron}{\boldsymbol{e}}}_1,{\overset{\textasciicaron}{\boldsymbol{e}}}_2,{\overset{\textasciicaron}{\boldsymbol{e}}}_3\rbrace </math> ; (b) proyección de los elementos ''k''  y ''l''  en el plano '''''p''''' , (c); polígono <math display="inline">P</math>  formado con el algoritmo de intersección; (d) división de <math display="inline">P</math>  en triángulos.                  

​

</span>

|}

​

===5.1. Algoritmo para el cálculo de los factores de peso===

​

En esta sección se describe el algoritmo utilizado para el cómputo de los factores de peso <math display="inline">n_{AB}^{\alpha }</math> .

* '''Para cada elemento no-mortar ''k''''' ,

* ''Formar el plano'''''''p'''''''con base  ''<math display="inline">V</math>''.''

* ''Proyectar las coordenadas de los nodos no-mortar  ''<math display="inline">\boldsymbol{x}_A^1\rightarrow {\mathbb{R}}^3</math>''en el plano'''''''p''''''', obteniendo los puntos:  ''<math display="inline">\boldsymbol{y}_A^1\rightarrow {\mathbb{R}}^2\in \tilde{k}</math>'', referido a  ''<math display="inline">V</math>'', ver figura''[[#fig0015|''3'']]''-b.''

* '''''Para cada elemento mortar''''''''''l''''' ,

* ''Proyectar las coordenadas de los nodos mortar  ''<math display="inline">\boldsymbol{x}_A^2\rightarrow {\mathbb{R}}^3</math>''en el plano'''''''p''''''', obteniendo los puntos:  ''<math display="inline">\boldsymbol{y}_A^2\rightarrow {\mathbb{R}}^2\in \tilde{l}</math>'', referido a  ''<math display="inline">V</math>''.''

* ''Utilizar el algoritmo de intersección para formar un polígono  ''<math display="inline">P=(\tilde{k}\cap \tilde{l})\subset \boldsymbol{p}</math>'', ver figura''[[#fig0015|''3'']]''-c'' .                                                  

* ''Ubicar el centro geométrico  ''<math display="inline">\boldsymbol{y}_P^C</math>''del polígono  ''<math display="inline">P</math> .                                                  

* ''Formar n''<sup>''P''</sup>''triángulos de vértices  ''<math display="inline">\boldsymbol{y}_I^P\in P</math>''y  ''<math display="inline">\boldsymbol{y}_P^C</math>'', ver figura''[[#fig0015|''3'']]''-d.''

* ''Cada triángulo P admite la siguiente parametrización:''

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{y}^P=\sum_{I=1}^3N_I^P(\boldsymbol{\xi })\quad \boldsymbol{y}_I^P,</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)

|}

​

''donde N''<sub>''I''</sub>''son las funciones de forma clásicas para un triángulo.''

* ''Utilizar los n''<sub>''g''</sub>''puntos de Gauss, conocidos para un triángulo lineal, y ubicarlos en cada triángulo P.''

* ''Para establecer una relación entre los puntos de Gauss de un triángulo P de un elemento no-mortar k o mortar l, es necesario calcular las coordenadas locales paramétricas  ''<math display="inline">{\boldsymbol{\xi }}_g^{\alpha }</math>''como''

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\sum_{I=1}^3N_I^P(\boldsymbol{\xi })\quad \boldsymbol{y}_I^P=</math><math>\sum_{A=1}^{3\overset{\textasciiacute}{o}4}N_A^{\alpha }({\boldsymbol{\xi }}_g^{\alpha })\quad \boldsymbol{y}_A^{\alpha }.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)

|}

​

''El sistema de ecuaciones algebraicas es resuelto por medio de un esquema de Newton-Raphson si los elementos k o l son cuadrangulares, o en forma directa si son triangulares. (En este último caso el sistema es lineal).''

* ''Para cada triángulo P, calcular los factores de peso definidos de la siguiente manera:''

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>n_{AB}^{1,P}=2\quad A^P\quad \sum_{g=1}^{n_g}N_A^1({\boldsymbol{\xi }}_g^1)N_B^1({\boldsymbol{\xi }}_g^1),</math>

|-

|<math>n_{AC}^{2,P}=2\quad A^P\quad \sum_{g=1}^{n_g}N_A^1({\boldsymbol{\xi }}_g^1)N_C^2({\boldsymbol{\xi }}_g^2),</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)

|}

​

''donde A''<sup>''P''</sup>''es el área de cada triángulo P.''

* ''Finalmente sumar todas las contribuciones  ''<math display="inline">n_{AB}^{\alpha ,P}</math>''de cada triángulo'' ,

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>n_{AB}^1=n_{AB}^1+\sum_{p=1}^{n_p}n_{AB}^{1,P},</math>

|-

|<math>n_{AC}^1=n_{AC}^1+\sum_{p=1}^{n_p}n_{AC}^{1,P}.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 22)

|}

​

==6. Matriz tangente==

​

La matriz tangente se obtiene a partir de la linealización del vector de fuerzas internas presentado en la ecuación [[#eq0085|(17)]] ,

​

<span id='eq0115'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta \boldsymbol{F}_B^{\alpha }=t_{N_A}{\boldsymbol{\nu }}_A\Delta n_{AB}^{\alpha }+</math><math>{\kappa }_An_{AB}^{\alpha }\left({\boldsymbol{\nu }}_A\otimes {\boldsymbol{\nu }}_A\right)\Delta \boldsymbol{g}_A</math>

|-

|<math>+{\kappa }_An_{AB}^{\alpha }\left(t_{N_A}\boldsymbol{I}+\right. </math><math>\left. {\boldsymbol{\nu }}_A\otimes \boldsymbol{g}_A\right)\Delta {\boldsymbol{\nu }}_A.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 23)

|}

​

La linealización del huelgo '''''g'''''<sub>''A''</sub>  de la ecuación [[#eq0115|(23)]]  se obtiene utilizando las ecuaciones ([[#eq0045|9]] ,[[#eq0075|15]] ) y resulta en

​

<span id='eq0120'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta \boldsymbol{g}_A=\sum_{B=1}^{n^1}\left(n_{AB}^1\Delta \boldsymbol{x}_B^1+\right. </math><math>\left. \Delta n_{AB}^1\boldsymbol{x}_B^1\right)</math>

|-

|<math>-\sum_{C=1}^{n^2}\left(n_{AC}^2\Delta \boldsymbol{x}_C^2+\right. </math><math>\left. \Delta n_{AC}^2\boldsymbol{x}_C^2\right).</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 24)

|}

​

Luego, para la derivación de <math display="inline">n_{AB}^{\alpha }</math> , se requiere derivar los factores de peso de cada triángulo ''P  ''  (<math display="inline">n_{AB}^{\alpha ,P}</math> ), y ensamblar cada contribución en <math display="inline">n_{AB}^{\alpha }</math> . Utilizando la ecuación [[#eq0045|(9)]] , la linealización de los factores de peso <math display="inline">n_{AB}^{\alpha ,P}</math>  puede ser escrita de la siguiente manera

​

<span id='eq0125'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta n_{AB}^{\alpha ,P}=\Delta {\int }_{{\Gamma }_c^{\alpha }}N_A^1({\boldsymbol{\xi }}^1)N_B^{\alpha }({\boldsymbol{\xi }}^{\alpha })\quad d{\Gamma }_c^1</math>

|-

|<math>=\sum_{g=1}^{ng}\Delta N_A^1({\boldsymbol{\xi }}^1)N_B^{\alpha }({\boldsymbol{\xi }}^{\alpha })+</math><math></math>

|-

|<math>N_A^1({\boldsymbol{\xi }}^1)\Delta N_B^{\alpha }({\boldsymbol{\xi }}^{\alpha })w_gA^P+</math><math>n_{AB}^{\alpha ,p}\Delta A^P/A^P.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 25)

|}

​

donde <math display="inline">w_g</math>  son los factores de peso. Para completar el desarrollo de la ecuación [[#eq0125|(25)]] , es necesario definir las linealizaciones de las funciones de forma <math display="inline">N_A^{\alpha }</math>  y la correspondiente al área de cada triángulo ''A''<sup>''P''</sup>  del polígono <math display="inline">P</math> .      

​

<span id='sec0040'></span>

===6.1. Funciones de forma===

​

La coordenada local de los puntos de Gauss de un triángulo ''P''  y un elemento no mortar o mortar son establecidas por la siguiente relación

​

<span id='eq0130'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\sum_{I=1}^3N_I^P(\boldsymbol{\xi })\quad \boldsymbol{y}_I^P=</math><math>\sum_{A=1}^{3\overset{\textasciiacute}{o}4}N_A^{\alpha }({\boldsymbol{\xi }}^{\alpha })\quad \boldsymbol{y}_A^{\alpha }.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 26)

|}

​

Tomando la derivada parcial de las funciones de forma en la dirección <math display="inline">{\xi }_{\beta }^{\alpha }</math> , con ''β''  = 1, 2,

​

<span id='eq0135'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta N_A^{\alpha }=\frac{\partial N_A^{\alpha }({\boldsymbol{\xi }}^{\alpha })}{\partial {\xi }_{\beta }^{\alpha }}\Delta {\xi }_{\beta }^{\alpha }=</math><math>N_{A,\beta }^{\alpha }\Delta {\boldsymbol{\xi }}_{\beta }^{\alpha },</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 27)

|}

​

y junto con la ecuación [[#eq0130|(26)]] , la variación de las coordenadas '''''ξ'''''<sup>''α''</sup>  para la superficie no-mortar o mortar en la dirección ''β''  se obtiene como

​

<span id='eq0140'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta {\xi }_{\beta }^{\alpha }={{\boldsymbol{\tau }}_{\beta }^{\alpha }}^T\cdot \left(\sum_{I=1}^3N_I^P(\boldsymbol{\xi })\Delta \boldsymbol{y}_I^P-\right. </math><math>\left. \sum_{A=1}^{3\overset{\textasciiacute}{o}4}N_A^{\alpha }({\boldsymbol{\xi }}^{\alpha })\Delta \boldsymbol{y}_A^P\right),</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 28)

|}

​

donde

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{{\boldsymbol{\tau }}_{\beta }^{\alpha }}^T={\left[\frac{\partial N_A^{\alpha }({\boldsymbol{\xi }}^{\alpha })}{\partial {\xi }_{\beta }^{\alpha }}\boldsymbol{y}_A^{\alpha }\right]}^{-1}.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 29)

|}

​

Notar que para la obtención de los coeficientes <math display="inline">{{\boldsymbol{\tau }}_{\beta }^{\alpha }}^T</math> , simplemente hay que invertir un sistema de 2 × 2. Con la elección de una base local como aquí se propone, se evita una definición de una base covariante que complica el desarrollo de las ecuaciones. Por lo tanto, con las ecuaciones ([[#eq0135|27]] ,[[#eq0140|28]] ) la variación de la función de forma evaluada en los puntos de Gauss <math display="inline">{\boldsymbol{\xi }}_g^{\alpha }</math>  puede ser escrita como

​

<span id='eq0150'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta N_A^{\alpha }({\boldsymbol{\xi }}_g^{\alpha })=</math><math>N_A^{\alpha }({\boldsymbol{\xi }}_g^{\alpha }){{\boldsymbol{\tau }}_{\beta }^{\alpha }}^T\cdot \sum_{I=1}^3N_I^P(\boldsymbol{\xi })\Delta \boldsymbol{y}_I^P-</math><math></math>

|-

|<math>\sum_{A=1}^{3\overset{\textasciiacute}{o}4}N_A^{\alpha }({\boldsymbol{\xi }}_g^{\alpha })\Delta \boldsymbol{y}_A^P.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 30)

|}

​

La ecuación [[#eq0150|(30)]]  depende de las variaciones de las coordenadas de los nodos de cada triángulo ''P  '' , (<math display="inline">\boldsymbol{y}_I^P\in P</math> ), y también de las variaciones de las coordenadas nodales proyectadas de los elementos no-mortar y mortar, (<math display="inline">\boldsymbol{y}_A^P\in P</math> ). Sin embargo, para derivar la matriz tangente es necesario transformar las ecuaciones para que sean expresadas solo en términos de las coordenadas de los nodos de los elementos no-mortar y mortar, <math display="inline">\boldsymbol{x}_A^{\alpha }\quad \in {\mathbb{R}}^3</math> , y ordenados según el vector '''''Φ'''''  descrito en la ecuación [[#eq0010|(2)]] .      

​

Se define <math display="inline">\overset{\mbox{ˆ}}{\boldsymbol{\Phi }}</math>  al vector de coordenadas nodales de 2 componentes de los elementos <math display="inline">\tilde{k}</math>  y <math display="inline">\tilde{l}</math>  referido a la base <math display="inline">V</math> ,

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\overset{\mbox{ˆ}}{\boldsymbol{\Phi }}=\left[\begin{array}{c}

\boldsymbol{y}_1^1\\

\boldsymbol{y}_2^1\\

\boldsymbol{y}_3^1\\

\boldsymbol{y}_4^1\\

\boldsymbol{y}_1^2\\

\boldsymbol{y}_2^2\\

\boldsymbol{y}_3^2\\

\boldsymbol{y}_4^2\\

​

\end{array}\right],</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 31)

|}

​

en forma similar se procede para otras topologías de elementos.      

​

La matriz que relaciona las variaciones de las coordenadas de los nodos de cada triángulo del polígono <math display="inline">P</math>  con el vector de coordenadas nodales <math display="inline">\overset{\mbox{ˆ}}{\Phi }</math>  es denominada <math display="inline">\boldsymbol{D}_{IK}^P</math> . Lo expresado anteriormente resulta en la siguiente expresión genérica,

​

<span id='eq0160'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta \boldsymbol{y}_I^P=\boldsymbol{D}_{IK}^P\Delta {\overset{\mbox{ˆ}}{\boldsymbol{\Phi }}}_{KL}=</math><math>\boldsymbol{D}_{IK}^{P1}\Delta \boldsymbol{y}_K^1+</math><math>\boldsymbol{D}_{IK}^{P2}\Delta \boldsymbol{y}_K^2,</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 32)

|}

​

donde la obtención de <math display="inline">\boldsymbol{D}_{IK}^P</math>  se detalla en el apéndice. Los vectores <math display="inline">\boldsymbol{y}_K^1</math>  y <math display="inline">\boldsymbol{y}_L^2</math>  denotan las coordenadas nodales de los elementos no-mortar y mortar, respectivamente, <math display="inline">\boldsymbol{D}_{IK}^{P1}</math>  y <math display="inline">\boldsymbol{D}_{IK}^{P2}</math>  son particiones de la matriz <math display="inline">\boldsymbol{D}_{IK}^P</math>  correspondientes a los lados no-mortar y mortar, respectivamente. Con la ecuación [[#eq0160|(32)]] , la linealización de las funciones de forma de la ecuación [[#eq0135|(27)]]  referida a las coordenadas de los nodos de los elementos no-mortar y mortar puede escribirse como

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta N_A^{\alpha }\left({\boldsymbol{\xi }}_g^{\alpha }\right)=</math><math>\boldsymbol{P}_A^{\alpha ,P}\left({\boldsymbol{\xi }}_g^{\alpha }\right)\cdot \Delta {\overset{\mbox{ˆ}}{\boldsymbol{\Phi }}}_{KL},</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 33)

|}

​

donde

​

<span id='eq0170'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{P}_A^{1,P}({\boldsymbol{\xi }}_g^1)=N_{A,\beta }^1\quad \left[(N_I^P({\boldsymbol{\xi }}^1)\boldsymbol{D}_{IK}^{P1}-\right. </math><math>\left. N_K^1({\boldsymbol{\xi }}_g^1)){{\boldsymbol{\tau }}_{\beta }^1}^T\mid N_I^P({\boldsymbol{\xi }}^1)\boldsymbol{D}_{IL}^{P2}{{\boldsymbol{\tau }}_{\beta }^1}^T\right],</math>

|-

|<math>\boldsymbol{P}_C^{2,P}({\boldsymbol{\xi }}_g^2)=N_{C,\beta }^2\quad \left[N_I^P({\boldsymbol{\xi }}^1)\boldsymbol{D}_{IL}^{P1}{{\boldsymbol{\tau }}_{\beta }^2}^T\mid (N_I^P({\boldsymbol{\xi }}^1)\boldsymbol{D}_{IK}^{P2}-\right. </math><math>\left. N_K^2({\boldsymbol{\xi }}_g^2)){{\boldsymbol{\tau }}_{\beta }^2}^T\right].</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 34)

|}

​

===6.2. Área de los triángulos===

​

La variación del área de cada triángulo ''P  ''  que pertenece al polígono <math display="inline">P</math> , es también requerida para completar el desarrollo de la ecuación [[#eq0125|(25)]] . El área de un triángulo ''P''  puede ser calculada con el producto vectorial de sus lados, es decir,

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>A^P=\frac{1}{2}\left\|\left[\left(\boldsymbol{y}_3^P-\right. \right. \right. </math><math>\left. \left. \left. \boldsymbol{y}_1^P\right)\times \left(-\right. \right. \right. </math><math>\left. \left. \left. \boldsymbol{y}_2^P+\boldsymbol{y}_1^P\right)\right]\right\|.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 35)

|}

​

Luego, teniendo en cuenta que cada triángulo pertenece al plano de proyección '''''p''''' , definido por la base <math display="inline">V</math> , la linealización de ''A''<sup>''P''</sup>  puede ser calculada de la siguiente manera

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta A^P=\frac{1}{2}\left[-Y_{32,2}^P\mid Y_{32,1}^P\mid -\right. </math><math>\left. Y_{13,2}^P\mid Y_{13,1}^P\mid -Y_{21,2}^P\mid Y_{21,1}^P\right]\cdot \left[\begin{array}{c}

\Delta \boldsymbol{y}_1^P\\

\Delta \boldsymbol{y}_2^P\\

\Delta \boldsymbol{y}_3^P\\

​

\end{array}\right],</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 36)

|}

​

donde <math display="inline">Y_{32,1}^P</math>  es la componente en la dirección <math display="inline">{\overset{\textasciicaron}{\boldsymbol{e}}}_1</math>  del vector <math display="inline">y_3^P-y_2^P</math>  que pertenece a la base <math display="inline">V</math> , <math display="inline">Y_{32,2}^P</math>  es la segunda componente del mismo vector pero en la dirección <math display="inline">{\overset{\textasciicaron}{\boldsymbol{e}}}_2</math> . En forma análoga se obtienen: <math display="inline">Y_{13,1}^P</math> , <math display="inline">Y_{13,2}^P</math> , <math display="inline">Y_{21,1}^P</math>  y <math display="inline">Y_{21,2}^P</math> .      

​

Teniendo en cuenta la ecuación [[#eq0160|(32)]]  que relaciona los coordenadas de los nodos del triángulo ''P  ''  con las coordenadas del elemento no-mortar y mortar (<math display="inline">\Delta \boldsymbol{y}_I^P=\boldsymbol{D}_{IK}^P\Delta \overset{\mbox{ˆ}}{\boldsymbol{\Phi }}</math> ), la linealización del área resulta en la siguiente ecuación

​

<span id='eq0185'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\frac{\Delta A^P}{A^P}=\boldsymbol{J}^P\boldsymbol{D}_{KL}^P\Delta {\overset{\mbox{ˆ}}{\boldsymbol{\Phi }}}_{KL},</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 37)

|}

​

con <math display="inline">\boldsymbol{J}^P=1/(2A^P)\left[-\right. </math><math>\left. Y_{32,2}^P\quad \vert \quad Y_{32,1}^P\quad \vert \quad -\right. </math><math>\left. Y_{13,2}^P\quad \vert \quad Y_{13,2}^P\quad \vert \quad -\right. </math><math>\left. Y_{21,2}^P\quad \vert \quad Y_{21,1}^P\right]</math> .      

​

===6.3. Factores de peso===

​

La linealización de los factores de peso evaluados en <math display="inline">{\boldsymbol{\xi }}_g^{\alpha }</math>  puede ser calculada con las ecuaciones ([[#eq0125|25]] ,[[#eq0170|34]] ,[[#eq0185|37]] ) como sigue

​

<span id='eq0190'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta n_{AB}^{\alpha ,P}=\boldsymbol{P}_{AB}^{\alpha ,P}\cdot \Delta {\overset{\mbox{ˆ}}{\boldsymbol{\Phi }}}_{KL},</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 38)

|}

​

con

​

<span id='eq0195'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{P}_{AB}^{\alpha ,P}=\sum_{g=1}^{n_g}N_B^{\alpha }({\boldsymbol{\xi }}_g^{\alpha })\quad \boldsymbol{P}_A^{1,P}({\boldsymbol{\xi }}_g^1)+</math><math></math>

|-

|<math>N_A^{\alpha }({\boldsymbol{\xi }}_g^{\alpha })\boldsymbol{P}_B^{\alpha ,p}({\boldsymbol{\xi }}_g^{\alpha })w_gA^P+</math><math>n_{AB}^{\alpha ,P}\boldsymbol{J}^P\boldsymbol{D}_{KL}^P.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 39)

|}

​

===6.4. Huelgo===

​

La variación de la función huelgo '''''g'''''<sub>''A''</sub>  se linealiza en cada triángulo y luego se suman todas las contribuciones en forma completamente análoga a lo descrito para los factores de peso. Por medio de la ecuación [[#eq0120|(24)]]  y las ecuaciones ([[#eq0190|38]] ,[[#eq0195|39]] ), la linealización de '''''g'''''<sub>''A''</sub>  para un triángulo ''P''  puede ser expresada como

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta \boldsymbol{g}_A^P=\left[\boldsymbol{N}_P+{\overline{\boldsymbol{N}}}_P\right]\Delta {\overset{\mbox{ˆ}}{\boldsymbol{\Phi }}}_{KL},</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 40)

|}

​

con

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{N}_P=\left[n_{A1}^{1,P}\boldsymbol{I},\ldots n_{A4}^{1,P}\boldsymbol{I},n_{A1}^{2,P}\boldsymbol{I}\ldots n_{A4}^{2,P}\boldsymbol{I}\right],</math>

|-

|<math>{\overline{\boldsymbol{N}}}_P=\sum_{D=1}^4\left(\boldsymbol{x}_D^1\otimes \boldsymbol{P}_{AD}^{1,P}-\right. </math><math>\left. \boldsymbol{x}_D^2\otimes \boldsymbol{P}_{AD}^{2,P}\right).</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 41)

|}

​

===6.5. Vector normal promedio===

​

La linealización del vector normal promedio para cada nodo de un vértice ''A''  de un determinado grupo de elementos no mortar, al que se denominará '''''ν'''''<sub>''A''</sub>  es calculado de la siguiente manera

​

<span id='eq0210'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\boldsymbol{\nu }}_A=\frac{{\sum }_j^{n_A}\boldsymbol{n}_j}{\parallel \boldsymbol{n}_j\parallel },</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 42)

|}

​

(ver [[#fig0010|fig. 2]] ), donde

​

<span id='eq0215'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>n_j=\frac{\boldsymbol{v}_j\times \boldsymbol{v}_j+1}{\left\|\boldsymbol{v}_j\times \boldsymbol{v}_j+1\right\|}\mbox{.}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 43)

|}

​

Aquí '''''n'''''<sub>''j''</sub>  es el vector normal del elemento no-mortar ''j  '' , y los vectores <math display="inline">\boldsymbol{v}_j</math>  y <math display="inline">\boldsymbol{v}_{j+1}</math>  denotan la diferencia de las coordenadas nodales referidas al nodo ''A'' , esto es,

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\begin{array}{c}

\boldsymbol{v}_{j+1}=\boldsymbol{x}_C^1-\boldsymbol{x}_A^2\\

\boldsymbol{v}_{j+1}=\boldsymbol{x}_B^1-\boldsymbol{x}_A^2

\end{array}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 44)

|}

​

Véase en la [[#fig0010|figura 2]]  para una mejor interpretación.      

​

La linealización de la ecuación [[#eq0210|(42)]]  resulta en

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta {\boldsymbol{\nu }}_A=\frac{{\sum }_j^{n_A}\Delta \boldsymbol{n}_j(t)\parallel \boldsymbol{n}_j\parallel -{\sum }_j^{n_A}\boldsymbol{n}_j\Delta \parallel \boldsymbol{n}_j\parallel }{\parallel \boldsymbol{n}_j{\parallel }^2}.</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 45)

|}

​

Luego, la linealización del vector normal '''''n'''''<sub>''j''</sub>  se obtiene a través de la ecuación [[#eq0215|(43)]]  como

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\Delta n_j=\frac{\left[I-n_j\otimes n_j\right]}{\left\|v_j\times v_{j+1}\right\|}\overset{\mbox{ˆ}}{\boldsymbol{W}}\Delta {\boldsymbol{\Phi }}_A=</math><math>\boldsymbol{S}\overset{\mbox{ˆ}}{\boldsymbol{W}}\Delta {\boldsymbol{\Phi }}_A</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 46)

|}

​

con '''Φ'''<sub>''A''</sub>  el vector de coordenadas '''''x'''''<sup>1</sup>  de los primeros nodos «vecinos» del nodo ''A''  y las matrices,

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{S}=\frac{\left[\boldsymbol{I}-\boldsymbol{n}_j\otimes \boldsymbol{n}_j\right]}{\parallel \boldsymbol{v}_j\times \boldsymbol{v}_{j+1}\parallel },</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 47)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\overset{\mbox{ˆ}}{\boldsymbol{W}}=\left[{\tilde{\boldsymbol{V}}}_{j+1}\quad \vert \quad -\right. </math><math>\left. {\tilde{\boldsymbol{V}}}_j\quad \vert \quad -\right. </math><math>\left. {\tilde{\boldsymbol{V}}}_{j+1}\quad \vert \quad {\tilde{\boldsymbol{V}}}_j\right],</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 48)

|}

​

donde <math display="inline">{\tilde{\boldsymbol{V}}}_j</math>  corresponde a la matriz antisimétrica asociada al vector <math display="inline">\boldsymbol{v}_j</math> . Luego puede demostrarse que <math display="inline">\Delta \parallel \boldsymbol{n}_j\parallel =</math><math>\boldsymbol{n}_j/\parallel \boldsymbol{n}_j\parallel \Delta \boldsymbol{n}_j=</math><math>{\overset{\textasciicaron}{\boldsymbol{n}}}_jS\overset{\mbox{ˆ}}{W}</math> , con lo cual la variación del vector normal promedio es expresado de la siguiente manera