mjhjmtu

Dear wolfram community,

I hope my problem is clear and easy to solve.

I have already solved the following heat equation over a domain:

Clear["Global`*"]
Needs["NDSolve`FEM`"];

pde = 1/r D[r*D[T[t, r, z], r, 1], r, 1] + 
 D[D[T[t, r, z], z, 1], z, 1] - D[T[t, r, z], t];

[CapitalDelta]z = 50*10^(-4);(*[m]*)
R = 5*10^(-2);(*[m]*)
[Tau] = 3.5*10^(-3)(*[s]*);
Tw = 200;(*[K]*)
[CapitalOmega] = 
 ImplicitRegion[0 <= r <= R && -2*[CapitalDelta]z <= z <= 0, r, z];

bc = DirichletCondition[T[t, r, z] == Tw, z == -2*[CapitalDelta]z], 
 DirichletCondition[T[t, r, z] == Tw, r == R];
(*Note that the unspecified boundaries are set to Neuman zero by 
default*)

ic = T[0, r, z] == 300 (*[K]*);


sol = NDSolve[pde == 0, bc, ic, 
 T, t, 0, [Tau], r, z [Element] [CapitalOmega]]

I optained an Interpolationfunction.

Manipulate[(*Plot mit relativen Farbspectrum*)
 Plot3D[T[t, r, z] /. sol, r, 0, R, z, -2*[CapitalDelta]z, 0,
 PlotRange -> 180, Tliq*1.1, Mesh -> None,
 Axes -> True, AxesLabel -> "r", "z", "T" , 
 AxesStyle -> Thickness[0.001],
 ColorFunction -> "TemperatureMap",
 (*ColorFunction[Rule](ColorData["TemperatureMap"][#3]&),*)
 ColorFunctionScaling -> True,
 ImageSize -> Large,
 PlotLegends -> Automatic],
 t, 0, [Tau], [Tau]*0.01]

Now I want to use the temperature distribution T[[Tau],r,z] as part of the initial value (T[0,r,z]) for a new domain (omeganew ) and then solve the heat equation over omeganew. The new domain is essentially the old domain but enlarged by an amount of 4*Δz:

 omeganew = ImplicitRegion[0 <= r <= R && -4*Δz <= z <= 0, r, z];

Here is the tricky part:
For the new solution I want the Temperature distribution from the previous solution to be "shifted" to the lower part of omeganew and a different function fc[r,z] (lets assume for the purpose of simplicity that fc[r,z]= constant) to be valid within the region that is newly added. So I need a shifted version of the previous solution with respect to the z-coordinate.

In a nutshel, I would like to extract the information (Temperaturevalues at an abitrary time t) obtained by NDsolve (e.g. sol) to construct a Unitstep function (is there a better alternative?) that I would then use as the initial condition for the new domain (omeganew).

Can anyone please help me?

Here is a way to do it. Let' set up the model:

Δz = 50*10^(-4);(*[m]*)R = 5*10^(-2);(*[m]*)τ = 
 3.5*10^(-2);(*[s]*)
pde = 1/r D[r*D[T[t, r, z], r, 1], r, 1] + D[D[T[t, r, z], z, 1], z, 1] - D[T[t, r, z], t];
Ω = ImplicitRegion[0 <= r <= R && -2*Δz <= z <= 0, r, z];
bc = DirichletCondition[T[t, r, z] == Tw, z == -2*Δz], DirichletCondition[T[t, r, z] == Tw, r == R];

If you now call NDSolveValue you will get a solution (looks like it's zero but that is a different issue)

sol = NDSolveValue[pde == 0, bc, T[0, r, z] == 0, T, t, 0, τ, r, z ∈ Ω];

If you evaluate the inteprolating function out side of the region you will get a warning and an Indetermiante as an answer.

sol[0, -1, 3]

Indeterminate

To change that you can use:

sol = NDSolveValue[pde == 0, bc, T[0, r, z] == 0, 
 T, 
 t, 0, τ, 
 r, z ∈ Ω, 
 "ExtrapolationHandler" -> 5 &, "WarningMessage" -> False
 ];

Now you will get the extrapolation value specified (5) and no warning:

sol[0, -1, 3]
5

With this you can then call NDSolveValue on a different domain with a different initial value like so:

sol2 = NDSolveValue[pde == 0, bc, T[0, r, z] == sol[τ, r, z], 
 T, 
 t, 0, τ, 
 r, z ∈ ImplicitRegion[ 0 <= r <= R && -10*Δz <= z <= 0, r, z]
 ];