B.3 Propiedades del material. Modelo elasto-plástico cíclico de Hujeux (ECP)

Empleado por los elementos sólidos 2D y 3D (NPAR(1) = 2 ó 3).

Para este modelo se tiene: MODEL=NPAR(15)=2 $\longrightarrow$ modelo elasto-plástico de Hujeux.

El modelo 2 es un modelo no lineal.

El programa considera NCON=36 y JDETAT=35.

Se requieren 6 líneas de datos

lectura en EFNL2 y EFNL3
	READ (LECG,1010) (PROP(j,n), j=1,NCON)
1010	FORMAT(8F10.0)
	READ (LECG,1070) (INIPEL(j), j=1,4)
1070	FORMAT(4I5)

nota	columnas	variable	definición
ligne 1:
1	1-10	PROP(1,n)	(XKI) módulo de compresión isotrópico bajo la presión de referencia PREF(=PROP(36)) ($K_i$)
	11-20	PROP(2,n)	(XGI) módulo de corte bajo la presión de referencia PREF(=PROP(36)) ($G_i$)
	21-30	PROP(3,n)	(XN) exponente de la ley elástica ($n$)
	31-40	PROP(4,n)	(PHI) ángulo de fricción en plasticidad perfecta ($\phi$) (en grados)
	41-50	PROP(5,n)	(PSI) ángulo de dilatancia (en grados) ($\psi$)
	51-60	PROP(6,n)	(BETA) coeficiente de la ley de estado crítico, que caracteriza el endurecimiento en densidad. ($\beta$)
	61-70	PROP(7,n)	(PCI) presión crítica inicial ($p_{ci}$)
	71-80	PROP(8,n)	(A) coeficiente de la ley de endurecimiento de los mecanismos desviadores ($a_m$)
ligne 2:
2	1-10	PROP(9,n)	(B) coeficiente que caracteriza la organización de los granos empleado en la expresiones de la función de fluencia desviadora ($b$)
	11-20	PROP(10,n)	(ACYC) coeficiente $a$ cíclico ($a_{cyc}$)
	21-30	PROP(11,n)	(ALFA) coeficiente para el dominio de comportamiento (dilatancia) ($\alpha$)
	31-40	PROP(12,n)	(RAYELA) radio del dominio elástico ($r_{ela}$) $0 < r_{ela} < r_{hys} < r_{mbl} \leq 1$
	41-50	PROP(13,n)	(RAYHYS) radio del dominio histerético ($r_{hys}$ ) $0 < r_{hys} \leq 1$
	51-60	PROP(14,n)	(RAYMBL) radio del dominio movilizado ($r_{mbl}$) $0 < r_{mbl} \leq 1$
	61-70	PROP(15,n)	(C) coeficiente de endurecimiento del mecanismo isótropo ($c$)
	71-80	PROP(16,n)	(D) coeficiente del mecanismo isótropo que relaciona la recta de estado crítico y la recta de consolidación isótropa ($d$)
ligne 3:
	1-10	PROP(17,n)	(CCYC) coeficiente $c$ en cíclico ($c_{cyc}$)
	11-20	PROP(18,n)	(DLTELA) radio del dominio elástico del mecanismo isótropo ($r_4$)
	21-30	PROP(19,n)	(XKIMIN) coeficiente de la ley de endurecimiento isótropo valor usual = 0
	31-40	PROP(20,n)	(XM) coeficiente del dominio de comportamiento ($m$)
	41-50	PROP(21,n)	(FACINC) factor para el tamaõ de los sub-incrementos valor usual = 0.2
	51-60	PROP(22,n)	(XKAUX) módulo de compresión isótropa para el operador auxiliar $\texttt{XKAUX} \geq \texttt{XKI}$
3	61-70	PROP(23,n)	(XGAUX) módulo de corte para el operador auxiliar $\texttt{XGAUX} \geq \texttt{XGI}$
	71-80	PROP(24,n)	(IECOUL)indicador del tipo de regla de flujo valores posibles :$\pm 1$,$\pm 2$,$\pm3$ valor usual = 1
ligne 4:
	1-10	PROP(25,n)	(INCMAX) cantidad máxima de sub-incremetos
	11-20	PROP(26,n)	(AKS(1)) permeabilidad saturada a la porosidad POROS0 den la dirección $x$ (no se emplea en 2D)
	21-30	PROP(27,n)	(AKS(2)) permeabilidad saturada a la porosidad POROS0 en la dirección $y$
	31-40	PROP(28,n)	(AKS(3)) permeabilidad saturada a la porosidad POROS0 en la dirección $z$
	41-50	PROP(29,n)	(AK0) coeficiente de empuje de tierras inicial ($K_0$)
	51-60	PROP(30,n)	= 0 ; PROP(30) se emplea para guardar el código de inicialización de los mecanismos, se calcula en función de INIPEL.
4	61-70	PROP(31,n)	(SATR) saturación residual ($S_r$) en estática y para un problema acoplado no saturado, donde (IRTNU) número del nivel de agua que cruza al grupo de elementos en estática y en mecánica pura.
4	71-80	PROP(32,n)	(ALFAS) coeficiente de la ley $S(p)$ en estática para un problema acoplado no saturado.
ligne 5:
4	1-10	PROP(33,n)	(DSDPMX) determinado por el programa, máximo de $\frac{d(sat)}{dp}$
5	11-20	PROP(34,n)	(FACNAP) coeficiente de una eventual napa freática inicial
5	21-30	PROP(35,n)	(ZNAPPE) cota de la napa freática
	31-40	PROP(36,n)	(PREF) presión de referencia, por defecto igual a $-1.\times10^6.$
ligne 6:
6	1-5	INIPEL(1)	inicialización mecanismo desviador 1, plano $yz$
	6-10	INIPEL(2)	inicialización mecanismo desviador 2, plano $zx$
	11-15	INIPEL(3)	inicialización mecanismo desviador 3, plano $xy$
	16-20	INIPEL(4)	inicialización mecanismo isótropo

notas

1. El modelo de Hujeux cíclico de 4 mecanismos (versión 1986) es un modelo de comportamiento elasto-plástico para suelos, basado en los concepto siguientes:
   1. emplea el principio de estado crítico modificado para Mohr-Coulomb
   2. la disipación plástica está dividida en 4 mecanismos:
      • 3 mecanismos en deformaciones planas sobre 3 planos ortogonales que coinciden con los ejes de referencia.
        La rotación de estos planos debido al endurecimiento no es considerada en este modelo de comportamiento.
        La superficie de carga para cada mecanismo desviador está definida por:
        $f_k=q_k-\sin \phi p_k \left( 1-b \log \frac{ppp}{p_c} \right) r_k$
        con:
        • $p_c=p_{ci} \exp b e_v^p$ (endurecimiento en densidad)
        • $ppp=p_k$ si IECOUL>0
        • $ppp=p$ si IECOUL<0
        y: $\partial_t r_k = \lambda_k \frac{\left( 1-r_k\right)^2}{a_k}$
        donde: $a_k = a_{cyc} + \left( a_m - a_{cyc} \right) \alpha_k$
        $\alpha_k = \texttt{PROP(11)} \left( \frac{r_k-r_{hys}}{r_{mbl}-r_{hys}} \right)^m$
        La elección de PROP(11)=ALFA es muy importante en relación a la dilatancia del material.
        • PROP(11)=0, dilatancia bloqueada
        • PROP(11)>0, dilatancia
        La ley de dilatancia puede ser seleccionada entre las 3 hipótesis siguiente con el dato de JECOUL (=ABS(IECOUL) ) :
        • si JECOUL=1, $\frac{\partial_t \varepsilon_v^{p_k}}{\partial_t \varepsilon^{p_k}} = \sin \psi - \frac{q_k}{p_k} \cos \theta$
          donde $\theta$ es el ángulo de los esfuerzos y de los incrementos de deformaciones plásticas (=0 en monótono);
        • si JECOUL=2, $\frac{\partial_t \varepsilon_v^{p_k}}{\partial_t \varepsilon^{p_k}} = \sin \psi - \frac{q_k}{p_k}$ (Hipótesis de Roscoe)
        • si JECOUL=3, $\frac{\partial_t \varepsilon_v^{p_k}}{\partial_t \varepsilon^{p_k}} = \sin \psi$ (Hipótesis de Mohr-Coulomb)
        Las dos primera hipótesis difieren solamente en cíclico. La primera fue propuesta por Hujeux y la segunda por Lassoudière.
      • 1 mecanismo en compresión isótropa:
        • La superficie de fluencia está definida por: $f_c=p-d p_c r_4$
        • La regla de endurecimiento está definida por : $\partial_t r_4 = \lambda_c \frac{\left( 1 - r_4 \right)^2}{c}$
   3. La carga cíclica es considerada por el modelo, tanto para los 3 mecanismos desviadores como para el mecanismo isótropo, empleando el concepto de módulo plástico solamente a los módulos de endurecimiento intrínsecos de cada mecanismo, así como la noción de dominio de comportamiento para los mecanismos desviadores.
   4. La parte elástica tiene un comportamiento no lineal isótropo, de la forma:
      • $K=K_i \left( \frac{p}{p_{ref}} \right)^n$ (módulo de compresión isótropo)
      • $G = G_i \left( \frac{p}{p_{ref}} \right)^n$ (module de corte)
      donde $p_{ref}$ es la presión de referencia y $p$ es el esfuerzo isótropo.
2. Los valores límites de $b$ permiten obtener las superficies de fluencia de Cam-Clay y de Mohr-Coulomb:
   • $b = 0$, superficie de fluencia tipo Mohr-Coulomb
   • $b = 1$, superficie de fluencia tipo Cam-Clay
3. XKAUX y XGAUX son los parámetros empleados para calcular la matriz elástica auxiliar, que se emplea para resolver el sistema no lineal.
   De forma general, dar a XKAUX y XGAUX valores iguales o mayores que los valores de XKI y XGI. El incremento de estos valores puede mejorar la convergencia pero puede aumentar significativamente el tiempo de cálculo.
4. La variación del grado de saturación ($S$) y de la permeabilidad ($k$) en función de la porosidad ($n$) y de la presión de poros ($p_{int}$) se define mediante las evoluciones siguientes:
   • $S = S_r + \frac{1 - S_r}{1+\left( \alpha_s \frac{p_{int}}{ \frac{\rho_w}{q_z} } \right)^2 }$
     donde $S_r$ es la saturación residual y $\alpha_s$ caracteriza la curvatura de la función $p(S)$.
   • $K(S,n) = K_{sat}(n) k_r(S)$
     donde $k_r(S)$ es una función dependiente del grado de saturación y que varía entre $0$ y $1$ según la relación de van Genuchten:
     $k_r(S) = \left( \frac{S-S_r}{1-S_r} \right)^3$
   La permeabilidad saturada $K_{sat}$ es una función de la porosidad $n$ según el modelo :
   $K_{sat}(n) = K_{sat}(n_0) \frac{n^3}{n_0^3} \frac{\left(1-n_0 \right)^2}{\left( 1-n \right)^2}$
   donde $K_{sat}(n_0)$ (AKS(i)) es la permeabilidad saturada a la porosidad inicial $n_0$ (POROS0); que es una función de la viscosidad del fluido.
   Observación: en un modelo puramente mecánico, se puede desear modelar la presencia de un nivel de agua que atraviesa un material. En dicho caso, la densidad de masa del material es:
   • $(1-n) \rho_s$ por sobre el nivel de agua
   • $(1-n)(\rho_s-\rho_w)$ por debajo del nivel freático (i.e. masa buoyante)
   El nivel de agua, definido previamente en "cargas" (§4.5), se asigna por su número IRTNU en las propiedades PROP en lugar de la saturación residual que no se emplea en mecánica pura.
   Esta opción no puede ser empleada en dinámica.
5. El empleo de una napa freática sólo sirve para inicializar las presiones de poros PINT en los puntos de integración del grupo de elementos por medio de:
   PINT = FACNAP * DENW * ABS(GZ) * (ZIPT-ZNAPPE)
   donde

   ZIPT es la cota del punto de integración

   DENW es la densidad de masa del fluido

   GZ es la gravedad en la dirección $z$.

   Esto supone que las incógnitas nodales hayan sido inicializadas previamente (con ICON = 10 ó 12)
6. Para la inicialización de mecanismo existen dos posibilidades:
   1. INIPEL(i) = 1: el punto de esfuerzo normalizado se ubica sobre la superficie de fluencia, el programa calcula el radio $r_i$ monótono para ser plásticamente admisible, por lo tanto se ajuste la superficie de fluencia al valor de esfuerzo actual salvo si $r < r_{ela}$ . Entonces, se tiene :
      • Si $r_i \geq r_{ela}$, $\texttt{RAY(i,1)} = r_i$
      • Si $r_i < r_{ela}$, $\texttt{RAY(i,1)} = r_{ela}$ e INIPEL(i) = -1
      se está en rango elástico, por lo tanto no puede haber cambio de sentido (mecanismo cíclico sobrepuesto al monótono) o de carga plástica si se permanece en dicho dominio.
   2. INIPEL(i) = -1: El programa calcula el radio $r_i$ en monótono, correspondiente al estado de esfuerzos:
      • Si $r_i < r_{ela}$, $\texttt{RAY(i,1)} = r_{ela}$
      • Si $r_i \geq r_{ela}$ el programa se detiene
      En todos los caso donde INIPEL(i)=-1, no puede haber cambio de sentido de carga inicial o descarga. En primer lugar se debe extender la superficie de fluencia monótona mediante una carga (i.e. se debe pasar por un estado INIPEL(i).
   Observación:de forma general, dar INIPEL(i) = 1.