# Solving Piecewise Differential Equation using NDSolve (coupling at BC)

I am having some issues in dealing with a system of differential equations. I would like to solve a 1D diffusive heat equation across several regions with different material properties. I now have a working solution for a given region (thank you) but am not sure on how to take into account the interfaces and boundary conditions.

From other posts, it appears to be possible to use a piecewise function to take into account different materials, however, I am unable to get this to work. Also, I am unsure if the mesh has to be modified to account for potential boundary instabilities.

Needs["NDSolveFEM"]
parameters1 = {v -> 0.001, DT -> 0.000000143, rho -> 1000, c -> 2000,
eta -> 0.001002}; (*Region 1*)
parametersWall = {v -> 0, DT -> 0.000000143, rho -> 2000, c -> 10,
eta -> 0.001002}; (*Region 2*)
parameters2 = {v -> 0.001, DT -> 0.000000143, rho -> 1000, c -> 4200,
eta -> 0.001002};(*Region 3*)
outerR = 0.002; (*outer position*)
innerR = 0.001; (*inner position*)
WalinPos = 0.0012 ;(*inner position of wall*)
WaloutPos = 0.0014; (*outer position of wall*)
phi1 = 0.2;(*heat flux on the inner position*)
phi2 = 0.2 ;(*heat flux on outer position*)
Tmax = 0.1; (*Temporal Endpoint for the simulation*)
HEATIMP = (-D[u[r, t], t] -
v D[u[r, t], r] + (DT/r^2) D[r^2. D[u[r, t], r],
r] + (12. eta)/(rho c) (v/r) (v/r));
HeatTrue =
Piecewise[{{HEATIMP /. parameters1,
r < WalinPos}, {HEATIMP /. parametersWall,
WalinPos <= r <= WaloutPos}, {HEATIMP /. parameters2,
r > WaloutPos}}]
Neum1 = NeumannValue[phi1, r == innerR]
Neum2 = NeumannValue[phi2, r == outerR]
ic = {u[r, 0] == 0};
heatdist =
NDSolve[{HeatTrue == Neum1 + Neum2, ic},
u, {r, innerR, outerR}, {t, 0, Tmax},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {MaxCellMeasure -> 0.000001}}}}]
Plot3D[Evaluate[First[u[r, t] /. %]], {r, innerR, outerR}, {t,
0.00000, Tmax}, PlotRange -> Full,
AxesLabel -> {Position, Time, Temperature}]


My other issue is with dealing with coupled differential equations that are only coupled at a specific point and effectively serve as a boundary condition. The "Neum1" Boundary condition in the example would be given by another differential equation (coupling to other differential equation).

D[P[t],t] + P[t](R[t] - D[u[r,t],r]|r=innerR+ )


where P[t] and R[t] will be given by other differential equations (and the system solved together). Without going into detail how would I go about coupling only at that particular point r=innerR+?

If it helps, the heat part doesn't couple to anything else (P[t] is effectively generating a heat flow at the boundary) and is just used to model heat being removed from the system.

Figured it out. Easiest way is to use IF to specify values for particular regions.

Needs["NDSolveFEM"]
outerR = 0.002; (*outer position*)
innerR = 0.001; (*inner position*)
WalinPos = 0.0011;(*inner position of wall*)
WalOutPos = 0.0013; (*outer position of wall*)

vel1 = 0.001; vel2 = 0; vel3 = -0.001;
DT1 = 0.000000143; DT2 = 0.0143; DT3 = 0.000000143;
rho1 = 1000; rho2 = 1000; rho3 = 1000;
c1 = 4200; c2 = 0.000001; c3 = 4200;
eta1 = 0.001002; eta2 = 0.001002; eta3 = 0.001002;

vel[r_] = If[r < WalinPos, vel1, If[r < WalOutPos, vel2,vel3]];
DT[r_] = If[r < WalinPos, DT1, If[r < WalOutPos, DT2, DT3]] ;
rho[r_] = If[r < WalinPos, rho1, If[r < WalOutPos, rho2, rho3]] ;
c[r_] = If[r < WalinPos, c1, If[r < WalOutPos, c2, c3]] ;
eta[r_] = If[r < WalinPos, eta1, If[r < WalOutPos, eta2, eta3]] ;
phi1 = -0.2;(*heat flux on the inner wall*)
phi2 = 0 ;(*heat flux on outer wall*)
Tmax = 0.5; (*Temporal Endpoint for the simulation*)
HEATIMP = (-D[u[r, t], t] - vel[r] D[u[r, t], r] + (DT[r]/r^2) D[r^2. D[u[r, t], r],
r] + (12. eta[r])/(rho[r] c[r]) (vel[r]/r) (vel[r]/r))
Neum1 = NeumannValue[phi1, r == innerR]
Neum2 = NeumannValue[phi2, r == outerR]
ic = {u[r, 0] == 0};
heatdist = NDSolve[{HEATIMP == Neum1 + Neum2, ic}, u, {r, innerR, outerR}, {t, 0, Tmax},Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {MaxCellMeasure -> 0.0000005}}}}]
Plot3D[Evaluate[First[u[r, t] /. %]], {r, innerR, outerR}, {t,0.00000, Tmax}, PlotRange -> Full, AxesLabel -> {Position, Time, Temperature}] Still have an issue with the boundaries. I will formulate that question better before I ask it.