@Rouben Rostamian 

Yes, Rouben ... This modeling is complex and sometimes difficult to understand ...

Answering your additional question: each phase corresponds to the time interval where some behavior domain (elastic, plastic, viscoplastic) is formed.

In the attached figure, I present the compaction law H(t) (the evolution of the position of the upper boundary). The phases are highlighted in solid colors below this curve. At the end of each phase, the profile of the basin and the respective behavior domains stand out.

For the modeling data, at the beginning of phase 1 we have t=0My and H(0)=0m and at the end of it we have t=Te=5.1My and H(5.1)=He=509.2m. Phase 2 starts at t=Te=5.1My with H(5.1)=He=509.2m and ends at t=Tp=13.3My with H(13.3)=Hp=1097.4m.

Anyway, I think you must have tired of reading all my explanations. I appreciate your patience. Thank you

@Rouben Rostamian  

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

Thank you

@Rouben Rostamian  

Yes, it works. Thank you for your help.

@tomleslie 

Yes, you are right. Thanks for the good remarks about the problem.

@Carl Love 

I did not succeed . Maybe the Heaviside function can help.
Anyway, thank you very much for your help.

@Carl Love 

Well, after fixing some errors, I have one more problem. Now in relation to the amount of boundary conditions. I believe this is related to an earlier discussion (https://mapleprimes.com/questions/225219-Why-Does-Pdsolve-In-Numeric-Mode-Return).

Thank you for your help.

Test_1.mw

@Rouben Rostamian  

Yes, most expressions have an error. Thanks for your comments and for your help.

@Carl Love 

Yes, you're right. Thanks for your comments and for your help.

@acer 

It is true. You're absolutely right. I was careless in my placements.

Your routine works perfectly.

Thanks a lot for the help.

@acer 

Your response exceeded my expectations. Thank you.
However, in my version 18 the shadebetween function does not seem to exist. I checked the online help and saw that this function has been available since version 15. Can this be a bug in my program?

@Preben Alsholm 
Perfect ... Everything is working ... Thank you very much.

@Preben Alsholm 

I think the problem is still not resolved.

The boundary conditions of the ode system depend on H, which is updated within the loop (H [i + 1] = H [i] + dH).

So, with each loop we have the same ode system with different ics.

Thank you

@Preben Alsholm 

Yes, now it works very well... Perfect result, with a very short processing time...

Well, just to clarify a point... As you well observed, the loop closes at i = 742... This is natural, since the purpose of the routine is to find a value of t that fulfills the condition "if fvp> = 0 then "... After observation on i, the value of dt (specifically Ts) can be adjusted, causing the loop to be executed close to the value 20000, which improves the discretization of the problem ...

Thanks a lot for the help...

@Preben Alsholm 
As always, you are helping me a lot. Not only with the maple but also with the math. And for that, I thank you.

Well, I'm looking at your code suggestion, line by line. However, when executing it, the following error occurs in the command "resNew: = dsolve ({sys_ode} ...": "Error, (in dsolve / numeric / process_input) input system must be an ODE system, got independent variables { x, H-47.3135421221118} ".

Any tips on this?

@mmcdara 

The variable nn is equal to 20, which causes this loop to execute at this speed ... But I usually need nn = 10000 or nn = 20000 ... In these situations the time processing is much larger ...

Also, I still deal with 8 more loops of this type, with the same problem ... So processing time really is my primary issue ...

About the Bcs you mentioned, they are correct ...

Thank you