3
$\begingroup$

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["NDSolve`FEM`"]
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.

$\endgroup$

1 Answer 1

1
$\begingroup$

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

Needs["NDSolve`FEM`"]
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}] 

enter image description here

Still have an issue with the boundaries. I will formulate that question better before I ask it.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.