Thank you for your attention. I'll try to explain it better.
The problem is the modeling of particle compaction in sedimentary basins (reference: [Lemos, P. S. B .; Brüch, A .; Maghous, S.] Formulation of reference solutions for compaction process in sedimentary basins). In summary, a flow of sedimentary particles is deposited in a given region over geological time. These particles accumulate and being compressed by the action of gravity.
As a basic mechanical hypothesis, a constant rate of mass of sediments deposited per area (variable Mp) is considered. This corresponds to a small value since the geological time scale is in the order of millions of years.
As a basic thermal hypothesis, it is considered that the particles deposited at the top of the basin have an initial temperature T0, which corresponds to the approximate temperature of the seabed. At the base of the basin, we have a boundary condition of constant heat flow q0.
From a mechanical point of view, the sedimentary material may have an elastic behavior (phase 1), elastoplastic (phase 2) and elasto-viscoplastic (phase omitted in the attached file). This behavior depends on the degree of compaction of the particles. The boundary between these domains of behavior is determined by yield criterion (in phase 2, the viscoplastic criterion is verified in the variable fvpx0).
The variables that characterize the mechanical behavior are the vertical stretch of the particles (variables L, L1, LA), the irreversible component of the Jacobian of the particles (variables Jir, Jir1, JIR), the velocity of the particles (variables u, u1, U), and the horizontal and vertical stress components (variables sh, SH, sh1, sv, SV, sv1). Other mechanical characteristics are derived from these functions, such as porosity (variables phi, phi_t), plastic consolidation pressure (variables pc, pc1) and density of the particles (variables rho, rho1, rho_t). The temperature field of the particles corresponds to the variable T, which is a function of the spatial x1 and temporal t1 variables (T(x1,t1)).
For this model, I am considering a partial coupling between mechanical and thermal phenomena. The thermal properties (thermal conductivity k, specific heat c and density rho) are porosity dependent. The next step in the research is to put some mechanical properties dependent on temperature.
The nature of the sedimentation and compaction process makes the height of the basin (variable H) dependent on time (a function of Mp) and the thermo-mechanical behavior of the particles (expressed by the velocity of the particles) . Thus, an Eulerian approach is used for the material description of the problem.
The solution of the mechanical problem of phase 1, which corresponds to the elastic behavior of the material, is obtained analytically (analytical functions H, L, sh, sv, Jir and u) as well as the limit values at the end of this phase (variables He, Le, she , sve, Jire and ue). At the end of phase 1, the basin profile is recalculated and the mechanical functions are stored in L1, sh1, sv1, Jir1, u1, phi1, pc1, rho1.
The thermal problem of phase 1 is nonlinear and so I look for a numerical solution (k, c and rho vary with porosity). As a result, there is a time loop, where at each time t[i]=t[i-1]+dt the temperature field is recalculated. At time t[i], the thermal boundary conditions are: fixed temperature at the top of the basin T(H,t1)=T0 and constant flow at the base T(0,t1)=q0. The initial condition at each time t[i] is given by the function f1. This piecewise function determines that the particles located between the base and the height of the basin in the previous time (H[i-1] at t[i-1]) have as an initial condition the function f2, which corresponds to the temperature field calculated in the previous time (T[i-1]). Particles above this level (variable limit1) have an initial temperature equal to T0.
Phase 2 has no analytical solution for thermo-mechanical phenomena. Thus, a methodology similar to that of phase 1 is used, including the EDO system that characterizes the mechanical problem. At each time t[i] the mechanical EDO system is solved as well as the PDE of heat conduction. The height of the basin is then recalculated for the next time H=H+dH, where dH is the height increase (a function of time and the velocity of the particles).
I added comments to the routine. Problem_1.mw