8
$\begingroup$

Statement of problem

I'm writing a script to calculate the temperature of a 2D system with time-dependent heat deposition. The heat deposition is a square wave pulse of duration w = 10^-6 seconds, occurring at a rate of rep = 150Hz (every 6.67*10^-3 seconds). Without using options, NDSolveValue misses these temporal pulses. If I set MaxStepFraction -> w/MaxTime, the calculation takes much too long (and runs out of memory on my system).

Is there a way to tell NDSolveValue to use shorter time steps around the points in time at which the pulses turn on? I've tried using WhenEvent[] but have not had any success.

Background information

The system comprises a thin metal sheet with bottom face at y=0 with water flowing over the top of the metal sheet. The top 1/8 of the metal sheet is an interface region with enhanced thermal transport properties. The heat pulse has a Lorentzian (Cauchy) shape centered on x=0.

If the system performs the calculation with adequate time steps, the temperature in the solid should look something like this: enter image description here This is calculated using MaxTimeStep->w and th1=0.02

Here's the code (modified from here, running on V 11.1):

Constants

    cs = 134.;(* Tungsten Specific Heat (J/kg K)*)
    cf = 4187.; (* Water Specific Heat (J/kg K)*)
    ps = 0.0193;(* Tungsten Density (kg/cm^3)*)
    pf = 0.001; (* Water density (kg/cm^3)*)
    ks = 1.; (* Tungsten Thermal Conductivity (W/cm K)*)
    kf = 650.*10^-5; (* Water Thermal Conductivity (W/cm K)*)
    f = 0.1; (* Thermal Transfer Coefficient (Metal/Water) (W/cm^2 K)*)

Adjustable Parameters

    (* Geometric Sizes*)
    th1 = 0.02; (* Thickness of solid layer (cm) *)
    th2 = 0.02; (* Thickness of liquid layer considered (cm) *)
    th3 = 0.2; (* Actual total thickness of liquid layer (cm) *)
    thick = th1/8; (* Thickness of interface (cm) *)
    len = .3; (* Width of water channel (cm) *)
    ht = th1 + th2; (* total height of system *)

    (* Temperature Parameters *)
    T0 = 276.; (* Inlet/outlet temperature (K) *)
    Tstarts = 350.; (* Starting solid temperature (K) *)
    Tstartf = 276.; (* Starting liquid temperature (K) *)

    (* Water Flow Parameters *)
    windia = 0.1; (* Water inlet diameter (cm) *)
    r = 0.2; (* scale factor for flow curve *)
    v0 = 50.; (* Water velocity in channel (cm/s) *)

    (* Beam Parameters *)
    w = 2*^-6; (*Pulse Temperal length (s)*)
    rep = 150;(*Pulse Rep Rate (Hz)*)
    MaxT = 1/rep; (* Max exposure time (s) *)

    (*Parameters for heat deposition (MeV) *)
    as = 10^7;
    bs = 0.0368;
    ns = 1.662;

    af = 10^6;
    bf = 0.03923;
    nf = 1.703;

    (* Time shape of pulse *)
    pulse[t_] = Piecewise[{{1, Abs[t] < w}}, 0];
    pulsemod[t_] = pulse[Mod[t - w, 1/rep, -1/rep]];

    (* Instantaneous power deposited in each layer (J/s) *)
    qsb[x_, y_, t_] = as/(bs^ns + Abs[x]^ns)*pulsemod[t];
    qfb[x_, y_, t_] = af/(bf^nf + Abs[x]^nf)*pulsemod[t];

Build Mesh

    Needs["NDSolve`FEM`"]

    (* Coordinates of edges *)
    top = ht;
    bot = 0;
    left = -len/2;
    right = len/2;
    interfacef = th1;
    interfaces = th1 - thick;
    buffery = 1.5 thick; (* Thickness of modified mesh around interface *)
    bufferx = len/10; (* Thickness of modified mesh around beam *)

    (* Mesh scale constants (larger values makes finer mesh) *)
    meshf = 1;
    meshs = 1;
    meshint = 1;
    meshbuf = 2;
    (*Use associations for clearer assignment later*)
    bounds = <|inlet -> 1, hot -> 2, outlet -> 3|>;
    regs = <|solid -> 10, fluid -> 20, interface -> 15|>;

    (*Meshing Definitions*)
    (*Coordinates*)
    crds = {{left, bot}(*1*), {right, bot}(*2*), {right, top}(*3*), {left, top}(*4*), {left, interfacef}(*5*), {right, interfacef}(*6*), {left, interfaces}(*7*), {right, interfaces}(*8*)};
    (*Edges*)
    lelms = {{5, 4}(*left edge*)(*1*), {1, 2}(*bottom edge*)(*2*), {6, 
        3}(*3*), {2, 8}, {8, 6}, {3, 4}, {5, 6}, {7, 8}, {1, 7}, {7, 
        5}(*4*)};
    boundaryMarker = {bounds[inlet], bounds[hot], bounds[outlet], 4, 4, 4,
       4, 4, 4, 4};(*4 will be a default boundary*)
    bcEle = {LineElement[lelms, boundaryMarker]};
    bmesh = ToBoundaryMesh["Coordinates" -> crds, "BoundaryElements" -> bcEle];

    (*Identify Center Points of Different Material Regions*)
    fluidCenter = {(left + right)/2, th1 + th2/2};
    fluidReg = {fluidCenter, regs[fluid], (th2/meshf)^2};

    interfaceCenter = {(left + right)/2, interfaces + thick/2};
    interfaceReg = {interfaceCenter, regs[interface], (thick/meshint)^2};

    solidCenter = {(left + right)/2, bot + th1/2};
    solidReg = {solidCenter, regs[solid], (th1/meshs)^2};

    (* Create and refine mesh *)
    meshRegs = {fluidReg, interfaceReg, solidReg};
    mesh = ToElementMesh[bmesh, "RegionMarker" -> meshRegs, MeshRefinementFunction -> Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
          If[y > (interfaces + interfacef)/2 - buffery && 
            y < (interfaces + interfacef)/2 + buffery, 
           area > (thick/meshbuf)^2, area > (th2/meshf)^2]]]];

    (* Plot Mesh *)
    (* Show[{mesh["Wireframe"["MeshElementStyle" -> {FaceForm[Blue], FaceForm[Yellow], 
          FaceForm[Red]}, ImageSize -> Large]]}, PlotRange -> {{-20 thick, 
        20 thick}, {(interfaces + interfacef)/2 - 
         2 buffery, (interfaces + interfacef)/2 + 2 buffery}}] *)

Region Values

    (*Region Dependent Properties with Piecewise Functions*)
    k = Evaluate[
      Piecewise[{{kf, ElementMarker == regs[fluid]}, {ks, 
         ElementMarker == regs[interface] || 
          ElementMarker == regs[solid]}, {0, True}}]];
    p = Evaluate[
       Piecewise[{{pf, ElementMarker == regs[fluid]}, {ps, 
          ElementMarker == regs[interface] || 
           ElementMarker == regs[solid]}, {0, True}}]];
    c = Evaluate[
       Piecewise[{{cf, ElementMarker == regs[fluid]}, {cs, 
          ElementMarker == regs[interface] || 
           ElementMarker == regs[solid]}, {0, True}}]];
    vp = Evaluate[
       Piecewise[{{v0 (1 - ((y - (th1 + ht)/2)/r)^2), 
          ElementMarker == regs[fluid]}, {0, True}}]];

    qsp[x_, y_, t_] = 
      Evaluate[Piecewise[{{qsb[x, y, t], 
          ElementMarker == regs[interface] || 
           ElementMarker == regs[solid]}, {0, True}}]];
    qfp[x_, y_, t_] = 
      Evaluate[Piecewise[{{qfb[x, y, t], 
          ElementMarker == regs[fluid]}, {0, True}}]];

    (*fac increases heat transfer coefficient in interface layer*)
    fac = Evaluate[If[ElementMarker == regs[interface], f/thick, 0]];

    (*Neumann Conditions *)
    nvsolid = 0;
    nvfluid = 
      NeumannValue[-(tf[x, y, t] - T0)*v0*th2*len*pf*cf, 
       ElementMarker == bounds[outlet]];

    (*Dirichlet Conditions for the Left Wall*)
    dcfluid = 
      DirichletCondition[tf[x, y, t] == T0, 
       ElementMarker == bounds[inlet]];

    (*Balance Equations for Fluid and Solid Temperature*)
    fluideqn = 
      p c (D[tf[x, y, t], t] + vp D[tf[x, y, t], x]) - 
        k Inactive[Laplacian][tf[x, y, t], {x, y}] - 
        fac (ts[x, y, t] - tf[x, y, t]) - qfp[x, y, t] == nvfluid;
    solideqn = 
      p c D[ts[x, y, t], t] - k Inactive[Laplacian][ts[x, y, t], {x, y}] -
         fac (tf[x, y, t] - ts[x, y, t]) - qsp[x, y, t] == nvsolid;
    ics = ts[x, y, 0] == Tstarts;
    icf = tf[x, y, 0] == Tstartf;

Solve the System

    (* Setup timer for convenience *)
    MSz = w; (* Max time step *)
    tp = 0;
    t0 = AbsoluteTime[];
    rate := Quiet[tp/(AbsoluteTime[] - t0)];
    ProgressIndicator[Dynamic[tp/MaxT]]
    Print[Dynamic[ToString[tp] <> " / " <> ToString[N[MaxT]]]]
    Print[Dynamic[
      "Time Remaining: " <> ToString[Round[(MaxT - tp)/rate]] <> " s"]]

    (* Execute Solving *)
    ifun = NDSolveValue[{fluideqn, solideqn, dcfluid, ics, icf}, {tf, 
        ts}, {t, 0, MaxT}, {x, y} \[Element] mesh
       , StepMonitor :> (tp = t)
       , MaxStepSize -> MSz];

    (* Plot Result *)
    (* Plot[ifun[[1]][0, th1 + thick/2, t], {t, 0, MaxT}
     , PlotRange -> All]
    Plot[ifun[[2]][0, th1 - thick/2, t], {t, 0, MaxT}
     , PlotRange -> All] *)

Additional Notes:

The following is independent of the primary purpose for this post.

Internal Boundary Conditions

It doesn't seem obvious that I need to have two separate temperature functions - one for the solid and one for the liquid. However, MMA cannot handle internal boundary conditions. If it could, I would add a Robin boundary condition, namely

    rc = NeumannValue[-(temp[x, y, t] - temp[0, th1, t])*
       f/thick, ElementMarker == bounds[interface]]

and the single differential equation to solve would be 

    tempeqn = 
     D[temp[x, y, t], t] == 
      1/(p c) (Inactive[Div][
           k Inactive[Grad][temp[x, y, t], {x, y}], {x, y}] + 
          qp[x, y, t] + nv + rc) - vp D[temp[x, y, t], x]

pursuant to initial condition

ic = temp[x, y, 0] == 
   Tstartf + (Tstarts - 
       Tstartf) (1 - 
        Tanh[(y - (interfacef + interfaces)/2)/(.25 thick)])/2;

However, when attempting to do this, MMA produces the error 

    NDSolveValue::delpde: Delay partial differential equations are not currently supported by NDSolve.

Since the boundary condition expression is 

    k d/dx(u(x,y,t)) = h(u(x,y,t)-u(x,th1,t))

one could imagine taking a linear approximation to the temperature within the interface region. That is, adding a term to the differential equation that looks like

    (k-h(y-th1))d/dx(u(x,y,t)) = 0

However, I think such an assumption is not justified and will produce incorrect results.

WhenEvent[]

According to the help file, WhenEvent[] allows the user to change the value of a state variable. A state variable seems to mean the dependent variable. The following is something that one might expect to work, but does not:

    ifun = NDSolveValue[{fluideqn, solideqn, dcfluid, ics, icf, 
      WhenEvent[Abs[t - 1/rep] <= w, MSz -> w], 
      WhenEvent[Abs[t - 1/rep] > w, MSz -> 10 w]}, 
      {tf, ts}, {t, 0, MaxT}, {x, y} \[Element] mesh
      , MaxStepSize :> MSz];

...because 'delayed rule (:>)' cannot be used as a MaxStepSize. The above attempt produces the error

    NDSolveValue::ndmss: Value of option MaxStepSize -> MSz is not a positive real number or Infinity.

Replacing the 'delayed rule (:>)' with just 'rule (->)', we get the errors:

    NDSolveValue::wenset: Warning: the rule MSz->w will not directly set the state because the left-hand side is not a list of state variables.
    NDSolveValue::wenset: Warning: the rule MSz->10 w will not directly set the state because the left-hand side is not a list of state variables.

This is telling us that the WhenEvent[] event will not be evaluated because MSz is not a state variable. WhenEvent[] does not work with 'set (=)', so I cannot actually change the value of MSz. And even if I could, I'd need to be able to use 'delayed rule (:>)' i.e. MaxStepSize :> MSz. I think the solution to my question is not found in using WhenEvent[], or if it is, then it is not in using WhenEvent[] with MaxStepSize or MaxStepFraction

Improve this question
$\endgroup$
11
  • $\begingroup$ How can I check that a solution is good? In other words what are you expecting the result to look like? $\endgroup$
    – user21
    Aug 30, 2019 at 7:30
  • $\begingroup$ The fluid velocity in the sheet metal and in the enhanced interface region appear to be solids as there is no fluid velocity in those regions. If so, you should only need one temperature variable equation. The referenced case was a porous media problem where the fluid flowed through the porous region and the fluid and solid were at different temperatures. $\endgroup$
    – Tim Laska
    Aug 30, 2019 at 12:40
  • $\begingroup$ Hi @TimLaska! Thanks for commenting. I suppose I do only need one temperature variable. How might I tackle the thermal transfer? Keeping 'fac' as a parameter that is nonzero only in the interface, I think I would include the term '-fac (ts[x, interfaces, t] - ts[x, interfacef, t])' in the differential equation. This would increase the time rate of change of the temperature in the interface region, assuming the solid is always hotter than the fluid. $\endgroup$
    – BohemianTapestry
    Aug 30, 2019 at 16:07
  • $\begingroup$ @user21, I edited the question to show how the temperature should evolve if the calculation is accurate. $\endgroup$
    – BohemianTapestry
    Aug 30, 2019 at 16:24
  • $\begingroup$ @TimLaska, unfortunately Mathematica doesn't like this system. Trying to combine the two temperature variables into one as you recommended, the solver complains that it can't solve a delay-differential equation. Of course, this is not a DDE. The only thing that might make MMA think it is a DDE is the function Mod[t-w] in the temporal definition of the pulse. So now I'm trying to work through MMA's quirks... $\endgroup$
    – BohemianTapestry
    Aug 30, 2019 at 16:29

1 Answer 1

Reset to default
12
+50
$\begingroup$

Your question and your system is fairly complex and I would consider breaking it up into more manageable chunks. It is easier to get help that way. Your system contains multiple materials, thin layers, liquids and solids, convection-diffusion, transient pulses, etc., so there are a lot of interactions to sort out. I also recommend that you conduct a dimensional analysis as it can help you sort out the dominant regimes that are present in the system. That aside, this is not a complete answer, rather it shows some building blocks that might be useful.

The following shows how I broke up the tasks into four steps:

  • Pulsed heating using WhenEvent on a 0D model.
  • Structured quad meshing to reduce model size.
  • Combine structured quad meshing with WhenEvent on a layered conduction problem.
  • Add convection.

Perhaps the following concepts can be used to at least reduce the model size so that concepts can be tested on a shorter cycle.

Pulsed heating using WhenEvent on a 0D model

I do not use WhenEvent enough to be a pro with its usage. Therefore, I always start with a simple model to make sure my WhenEvent construction behaves as intended. Consider the following simple model of a flow tank heated by a pulsed coil as shown by the equation below.

$$\frac{{du}}{{dt}} = - u(t) + q(t)$$

In the following Mathematica code, I introduce a unit heat load with a period of one time unit with a duty cycle of 0.025.

duty = 1/40;
period = 1;
{sol} = NDSolve[{u'[t] == -u[t] + q[t], q[0] == 0, u[0] == 0, 
    WhenEvent[{Mod[t, period], 
      Mod[t + period duty, period]}, {q[t] -> 
       If[q[t] == 0, 1/duty, 0]}]}, {u, q}, {t, 0, 10}, 
   DiscreteVariables -> q];
Plot[{Evaluate[{u[t], q[t]} /. sol], 0, 1/duty}, {t, 0, 10}, 
 PlotTheme -> "Web", PlotStyle -> {Thick, Thick, Dashed, Dashed}, 
 PlotPoints -> 500]
Row[{
  Column[{
    Plot[{Evaluate[q[t] /. sol], 0, 1/duty}, {t, 0, 10}, 
     PlotTheme -> "Web", 
     PlotStyle -> {Directive[Thick, Green], Dashed, Dashed}, 
     PlotPoints -> 500, ImageSize -> Medium],
    Plot[{Evaluate[u[t] /. sol]}, {t, 0, 10}, PlotTheme -> "Web", 
     PlotStyle -> {Directive[Thick, Red]}, PlotPoints -> 500, 
     ImageSize -> Medium]
    }], Column[{
    Plot[{Evaluate[q[t] /. sol], 0, 1/duty}, {t, 0, 2.1}, 
     PlotTheme -> "Web", 
     PlotStyle -> {Directive[Thick, Green], Dashed, Dashed}, 
     PlotPoints -> 500, ImageSize -> Medium],
    Plot[{Evaluate[u[t] /. sol]}, {t, 0, 2.1}, PlotTheme -> "Web", 
     PlotStyle -> {Directive[Thick, Red]}, PlotPoints -> 500, 
     ImageSize -> Medium]
    }]}]

0D WhenEvent Pulsed Heating

The results looks similar to the OP so this looks like a working representation of a pulse sequence with WhenEvent.

Structured quad meshing to reduce model size

A good computational mesh is necessary for accurate simulation results. For a model such as this that contains thin layers and potentially very thin thermal boundary layers, one generally uses an anisotropic mesh that is fine in the direction of steep gradients and coarser in the direction of shallow gradients.Using this approach, you will have a much smaller mesh and potentially longer time steps due to CFL considerations thus substantially reducing your computational requirements.

Unfortunately, Mathematica does not provide a GUI to construct these types of mapped structured meshes. Fortunately, Mathematica provides lots of geometric computation that should allow us to slap something together to construct layered structured meshes. In fact, I was inspired by the RegionProduct documentation that shows how one can simply construct a tensor product grid with a graded mesh. This combined with the two Element Mesh Tutorial should give us what we need to construct a valid FEM mesh.

I apologize in advance for the following code. It is hastily constructed, but it appears to function and will allow us to construct structured layered meshes on rectangular domains with a few lines of code.

Mathematica code for structured meshes

Needs["NDSolve`FEM`"]
ex = {1, 0};
ey = {0, 1};
eleft = -ex;
eright = ex;
etop = ey;
ebot = -ey;
ebi = ElementIncidents[#["BoundaryElements"]][[1]] &;
ebm = ElementMarkers[#["BoundaryElements"]][[1]] &;
ei = ElementIncidents[#["MeshElements"]][[1]] &;
em = ElementMarkers[#["MeshElements"]][[1]] &;
epi = Flatten@ElementIncidents[#["PointElements"]] &;
epm = Flatten@ElementMarkers[#["PointElements"]] &;

(* Shortand *)
FP = Flatten@Position[#, True] &;
UF = Union@Flatten[#, Infinity] &;

gidx = Flatten@Position[#, True] &;
gelm = #1[[gidx[#2]]] &;
ginc = Union@Flatten@gelm[#1, #2] &;
getBoundaryNodes = ginc[#["pureBoundaries"], #[dirs[#2]]] &;

lineElms2Nodes[lelms_, mask_] := 
 Union@Flatten@lelms[[Flatten@Position[mask, True]]]
pfn[ei_, em_, marker_] := Pick[ei, # == marker & /@ em]
in1dMask[l1_, l2_] := MemberQ[l1, #] & /@ l2
in2dMasks[l1_, l2_] := in1dMask[l1, #] & /@ Transpose[l2]
inBothMask[l1_, l2_] := Inner[And, #1, #2, List] & @@ in2dMasks[l1, l2]
regBothMask[assoc_, marker_] := 
 inBothMask[assoc["regIncAssoc"][marker], assoc["pureBoundaries"]]
meshinfo[mesh_] := Module[{crd, nCrd, elms, nElms, markers, nMarkers,
   uniqueMarkers, boundaries, boundaryNormals, bndNodes, bndMarkers, 
   regInc, regIncAssoc},
  crd = mesh["Coordinates"];
  nCrd = Dimensions[crd][[1]];
  elms = ei[mesh];
  nElms = Dimensions[elms][[1]];
  markers = em[mesh];
  nMarkers = Dimensions[markers][[1]];
  uniqueMarkers = Union@markers;
  boundaries = ebi[mesh];
  boundaryNormals = mesh["BoundaryNormals"][[1]];
  bndNodes = epi[mesh];
  bndMarkers = epm[mesh];
  regInc = pfn[elms, markers, #] & /@ uniqueMarkers;
  regIncAssoc = 
   AssociationThread[uniqueMarkers -> (Union[Flatten@#] & /@ regInc)];
  <|
   "crd" -> crd,
   "nCrd" -> nCrd,
   "elms" -> elms,
   "nElms" -> nElms,
   "markers" -> markers,
   "nMarkers" -> nMarkers,
   "uniqueMarkers" -> uniqueMarkers,
   "boundaries" -> boundaries,
   "boundaryNormals" -> boundaryNormals,
   "bndNodes" -> bndNodes,
   "bndMarkers" -> bndMarkers,
   "regIncAssoc" -> regIncAssoc
   |>
  ]
extinfo[mesh_] := 
 Module[{flat, flatinfo , assoc, regBndList, regBoundMasks,
   pureBoundaryNormals, pureNorth, pureEast, pureSouth, pureWest},
  assoc = meshinfo[mesh];
  flat = flatMesh[mesh];
  flatinfo = meshinfo[flat];
  AppendTo[assoc, "pureBoundaries" -> flatinfo["boundaries"]];
  AppendTo[assoc, 
   "pureBoundaryMarkers" -> 
    First@ElementMarkers@flat["BoundaryElements"]];
  AppendTo[assoc, 
   "nPureBoundaries" -> Dimensions[flatinfo["boundaries"]][[1]]];
  AppendTo[assoc, "pureBndNodes" -> flatinfo["bndNodes"]];
  AppendTo[assoc, "pureBndMarkers" -> flatinfo["bndMarkers"]];
  pureBoundaryNormals = flat["BoundaryNormals"][[1]];
  AppendTo[assoc, "pureBoundaryNormals" -> pureBoundaryNormals];
  pureNorth = (0.9999 < ey.#) & /@ pureBoundaryNormals;
  pureEast = (0.9999 < ex.#) & /@ pureBoundaryNormals;
  pureSouth = (0.9999 < -ey.#) & /@ pureBoundaryNormals;
  pureWest = (0.9999 < -ex.#) & /@ pureBoundaryNormals;
  AppendTo[assoc, "pureNorth" -> pureNorth];
  AppendTo[assoc, "pureEast" -> pureEast];
  AppendTo[assoc, "pureSouth" -> pureSouth];
  AppendTo[assoc, "pureWest" -> pureWest];
  regBndList = regBothMask[assoc, #] & /@ assoc["uniqueMarkers"];
  regBoundMasks = 
   AssociationThread[assoc["uniqueMarkers"] -> regBndList];
  AppendTo[assoc, "regBoundMasks" -> regBoundMasks]
  ]
meshGrowth[x0_, xf_, n_, ratio_] := Module[{k, fac, delta},
  k = Log[ratio]/(n - 1);
  fac = Exp[k];
  delta = (xf - x0)/Sum[fac^(i - 1), {i, 1, n - 1}];
  N[{x0}~Join~(x0 + 
      delta Rest@
        FoldList[(#1 + #2) &, 0, 
         PowerRange[fac^0, fac^(n - 3), fac]])~Join~{xf}]
  ]
adjust[l_, assoc_] := 
 Module[{itest, newlist, nodesfound, newmarks, pos, ll},
  newlist = l["pbm"];
  itest = Inner[And, assoc["reg"], assoc["dir"], List];
  pos = Flatten@Position[itest, True];
  newlist[[pos]] = assoc["marker"];
  nodesfound = UF@assoc["lelm"][[pos]];
  ll = assoc["lnodes"];
  newmarks = l["pbnm"];
  newmarks[[Flatten@(Position[ll, #] & /@ nodesfound)]] = 
   assoc["marker"];
  <|"pbm" -> newlist, "pbnm" -> newmarks|>]
adjustMarkers[mesh_, adjustments_] :=
 Module[{itest, extmi, assocs, l, bcEle},
  extmi = extinfo[mesh];
  assocs = 
   AssociationThread[{"lelm", "lnodes", "reg", "dir", "marker"},
      {extmi["pureBoundaries"], extmi["pureBndNodes"], 
       extmi["regBoundMasks"][#["region"]], 
       extmi[#["dir"]], #["marker"]}] & /@ adjustments;
  l = <|"pbm" -> extmi["pureBoundaryMarkers"], 
    "pbnm" -> extmi["pureBndMarkers"]|>;
  l = Fold[adjust, l, assocs];
  bcEle = {LineElement[extmi["pureBoundaries"], l["pbm"]]};
  (*l=extmi["pureBndMarkers"];
  l=Fold[adjust,l,assocs];*)
  pEle = {PointElement[Transpose@{extmi["pureBndNodes"]}, l["pbnm"]]};
  {bcEle, 
   ToElementMesh["Coordinates" -> mesh["Coordinates"], 
    "MeshElements" -> mesh["MeshElements"], 
    "BoundaryElements" -> bcEle, "PointElements" -> pEle]}]
pointsToMesh[data_] :=
  MeshRegion[Transpose[{data}], 
   Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
rp2Mesh[rh_, rv_, marker_] := Module[{sqr, crd, inc, msh, mrkrs},
  sqr = RegionProduct[rh, rv];
  crd = MeshCoordinates[sqr];
  inc = Delete[0] /@ MeshCells[sqr, 2];
  mrkrs = ConstantArray[marker, First@Dimensions@inc];
  msh = ToElementMesh["Coordinates" -> crd, 
    "MeshElements" -> {QuadElement[inc, mrkrs]}]
  ]
combineMeshes[mesh1_, mesh2_] := 
 Module[{crd1, crd2, newcrd, numinc1, inc1, inc2, mrk1, mrk2, melms},
  crd1 = mesh1["Coordinates"];
  crd2 = mesh2["Coordinates"];
  numinc1 = First@Dimensions@crd1;
  newcrd = crd1~Join~ crd2;
  inc1 =  ElementIncidents[mesh1["MeshElements"]][[1]];
  inc2 =  ElementIncidents[mesh2["MeshElements"]][[1]];
  mrk1 = ElementMarkers[mesh1["MeshElements"]][[1]];
  mrk2 = ElementMarkers[mesh2["MeshElements"]][[1]];
  melms = {QuadElement[inc1~Join~(numinc1 + inc2), mrk1~Join~mrk2]};
  ToElementMesh["Coordinates" -> newcrd, "MeshElements" -> melms]
  ]
markerSubsets[mesh_] := With[
  {crd = mesh["Coordinates"], 
   bids = Flatten[ElementIncidents[mesh["PointElements"]]],
   ei = ei[mesh], em = em[mesh]},
  {crd, bids, ei, em, pfn[ei, em, #] & /@ Union[em]}]
incidentSubs[mesh_] := 
 Module[{coords, ei, em, boundaryIDs, pureboundaryIDs, mei, 
   interiorIDs, interfaceNodes},
  {coords, boundaryIDs, ei, em, mei} = markerSubsets[mesh];
  interiorIDs = Complement[Range[Length[coords]], boundaryIDs];
  interfaceNodes = 
   Flatten[Intersection @@ (Flatten[#] &) /@ # & /@ 
     Partition[mei, 2, 1]];
  pureboundaryIDs = Complement[boundaryIDs, interfaceNodes];
  {pureboundaryIDs, interfaceNodes, interiorIDs}
  ]
flatMesh[mesh_] := 
 ToElementMesh["Coordinates" -> mesh["Coordinates"], 
  "MeshElements" -> {QuadElement[
     ElementIncidents[mesh["MeshElements"]][[1]]]}]
nodeTypes[mesh_] := 
 Module[{mtemp, pureboundaryIDs, interfaceNodes, intIDs, 
   tpureboundaryIDs, tinterfaceNodes, tintIDs, boundaryInts, 
   interiorInterfaceNodes, bool},
  mtemp = flatMesh[mesh];
  {pureboundaryIDs, interfaceNodes, intIDs} = incidentSubs[mesh];
  {tpureboundaryIDs, tinterfaceNodes, tintIDs} = incidentSubs[mtemp];
  boundaryInts = Complement[tpureboundaryIDs, pureboundaryIDs];
  interiorInterfaceNodes = Complement[interfaceNodes, boundaryInts];
  bool = ContainsAll[tpureboundaryIDs, #] & /@ ebi[mesh];
  {bool, tpureboundaryIDs, interiorInterfaceNodes, intIDs}]
(*Use associations for clearer assignment later*)
bounds = <|"inlet" -> 1, "hot" -> 2, "outlet" -> 3, "cold" -> 4, 
   "default" -> 0|>;
regs = <|"solid" -> 10, "fluid" -> 20, "interface" -> 15, 
   "insulation" -> 100|>;
dirs = <|"north" -> "pureNorth", "east" -> "pureEast", 
   "south" -> "pureSouth", "west" -> "pureWest"|>;
bcadj = <|"region" -> regs[#1], "dir" -> dirs[#2], 
    "marker" -> bounds[#3]|> &;

The following constructs a thin ${\color{Red} {Red}}$ solid region with a uniform mesh and a thicker ${\color{Green} {Green}}$ fluid region with a boundary layer mesh to capture the solid fluid interface. I also marked certain edges by what I think there boundary conditions are going to be later. If they are not used, they default to Neumann value of zero or that of an insulated wall condition.

(* Model Dimensions *)
lf = 0;
rt = 5;
th1 = 2;
th2 = 8;
bt = -th1;
tp = th2;
(* Horizontal Flow Dir Region *)
rh = pointsToMesh[Subdivide[lf, rt, 10]];
(* Thin Metal Region Uniform Mesh*)
rv = pointsToMesh[Subdivide[bt, 0, 10]];
(* Thick Fluid Region Geometric Growth Mesh *)
rv2 = pointsToMesh@meshGrowth[0, tp, 40, 16];
(* Build Element Meshes From Region Products *)
m1 = rp2Mesh[rh, rv, regs["solid"]];
m2 = rp2Mesh[rh, rv2, regs["fluid"]];
(* Combine the solid and fluid mesh *)
mesh = combineMeshes[m1, m2];
(* Define a series of BC adjustments *)
(* Last assignement takes precedence with PointElement *)
adjustments = {bcadj["solid", "south", "hot"]};
AppendTo[adjustments, bcadj["fluid", "north", "cold"]];
AppendTo[adjustments, bcadj["fluid", "west", "inlet"]];
(* Adjust the mesh with new boundary and point elements *)
{bcEle, mesh} = adjustMarkers[mesh, adjustments];
(* Display the mesh and bc's *)
Column[{Row@{mesh[
     "Wireframe"["MeshElement" -> "BoundaryElements", 
      "MeshElementMarkerStyle" -> Blue, 
      "MeshElementStyle" -> {Black, Green, Red}, ImageSize -> Medium]],
    mesh["Wireframe"[
      "MeshElementStyle" -> {FaceForm[Red], FaceForm[Green]}, 
      ImageSize -> Medium]]},
  Row@{mesh[
     "Wireframe"["MeshElement" -> "PointElements", 
      "MeshElementIDStyle" -> Black, ImageSize -> Medium]],
    mesh["Wireframe"["MeshElement" -> "PointElements", 
      "MeshElementMarkerStyle" -> Blue, 
      "MeshElementStyle" -> {Black, Green, Red}, 
      ImageSize -> Medium]]}}]

Structure Quad Mesh

The images show that I constructed the mesh as I intended.

Combine structured quad meshing with WhenEvents on a layered conduction problem

Now, we are ready to combine the WhenEvent, structured mesh, and heat equation example from the finite element tutorial into an example where we pulse the solid layer with heat and watch it transfer into the fluid layer. For simplicity, we are considering conduction only and I have set the top of the model to be a cold wall at the initial starting temperature condition.

duty = 1/32;
period = 0.5;
fac = Evaluate[
   Piecewise[{{0.1, ElementMarker == regs["solid"]}, {0, True}}]];
k = Evaluate[
   Piecewise[{{285, ElementMarker == regs["solid"]}, {1, True}}]];
op = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[t, x, y]\)\) - 
   Inactive[
     Div][(-{{k, 0}, {0, k}}.Inactive[Grad][u[t, x, y], {x, y}]), {x, 
     y}] - fac q[t];
Subscript[\[CapitalGamma], D2] = 
  DirichletCondition[u[t, x, y] == 0, ElementMarker == bounds["cold"]];
ufunHeat = 
  NDSolveValue[{op == 0, u[0, x, y] == 0 , Subscript[\[CapitalGamma], 
    D2], q[0] == 0, 
    WhenEvent[{Mod[t, period], 
      Mod[t + period duty, period]}, {q[t] -> 
       If[q[t] == 0, 1/duty, 0]}, 
     "DetectionMethod" -> "Interpolation"]}, {u, q}, {t, 0, 
    5}, {x, y} \[Element] mesh, DiscreteVariables -> q, 
   MaxStepFraction -> 0.001];

This code should run in a few seconds. Due to the discretization differences between the layers, I find it is usually best to plot each layer separately and combined them with Show.

plrng = {{lf, rt}, {bt, tp}, {0, 0.320}};
SetOptions[Plot3D, PlotRange -> plrng, PlotPoints -> Full, 
  ColorFunction -> 
   Function[{x, y, z}, Directive[ColorData["DarkBands"][3 z]]], 
  ColorFunctionScaling -> False, MeshFunctions -> {#3 &}, Mesh -> 20, 
  AxesLabel -> Automatic, ImageSize -> Large];
plts = Plot3D[ufunHeat[[1]][#, x, y], {x, y} \[Element] m1, 
    MeshStyle -> {Black, Thick}] &;
pltf = Plot3D[ufunHeat[[1]][#, x, y], {x, y} \[Element] m2, 
    MeshStyle -> {Dashed, Black, Thick}] &;
showplot = 
  Show[{plts[#], pltf[#]}, 
    ViewPoint -> {3.252862844243345`, 0.28575764805522785`, 
      0.8872575066569075`}, 
    ViewVertical -> {-0.2612026545717462`, -0.022946143077719586`, 
      0.9650112163920842`}, ImageSize -> 480, 
    Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False] &;
ListAnimate[showplot /@ Evaluate@Subdivide[0, 5, 80]]

Pulsed Conduction

The results appear to be reasonable.

Add convection

Now, we are in a position to add the convective term to the fluid layer. I will begin by making the flow length four times longer and I will increase the resolution at the fluid-solid interface using the following code. The fluid enters via the inlet at the initial conditions.

(* Model Dimensions *)
lf = 0;
rt = 20;
th1 = 2;
th2 = 8;
bt = -th1;
tp = th2;
(* Horizontal Region *)
rh = pointsToMesh[Subdivide[lf, rt, 40]];
(* Thin Metal Region Uniform Mesh*)
rv = pointsToMesh[Subdivide[bt, 0, 10]];
(* Thick Fluid Region Geometric Growth Mesh *)
rv2 = pointsToMesh@meshGrowth[0, tp, 80, 32];
(* Build Element Meshes From Region Products *)
m1 = rp2Mesh[rh, rv, regs["solid"]];
m2 = rp2Mesh[rh, rv2, regs["fluid"]];
(* Combine the solid and fluid mesh *)
mesh = combineMeshes[m1, m2];
(* Define a series of BC adjustments *)
(* Last assignement takes precedence with PointElement *)
adjustments = {bcadj["solid", "south", "hot"]};
AppendTo[adjustments, bcadj["fluid", "north", "cold"]];
AppendTo[adjustments, bcadj["fluid", "west", "inlet"]];
(* Adjust the mesh with new boundary and point elements *)
{bcEle, mesh} = adjustMarkers[mesh, adjustments];
(* Display the mesh and bc's *)
Column[{Row@{mesh[
     "Wireframe"["MeshElement" -> "BoundaryElements", 
      "MeshElementMarkerStyle" -> Blue, 
      "MeshElementStyle" -> {Black, Green, Red}, ImageSize -> Medium]],
    mesh["Wireframe"[
      "MeshElementStyle" -> {FaceForm[Red], FaceForm[Green]}, 
      ImageSize -> Medium]]},
  Row@{mesh[
     "Wireframe"["MeshElement" -> "PointElements", 
      "MeshElementIDStyle" -> Black, ImageSize -> Medium]],
    mesh["Wireframe"["MeshElement" -> "PointElements", 
      "MeshElementMarkerStyle" -> Blue, 
      "MeshElementStyle" -> {Black, Green, Red}, 
      ImageSize -> Medium]]}}]
(* Simulation *)
duty = 1/32;
period = 0.5;
v = Evaluate[
   Piecewise[{{{0.1 (y/th2)^2 {1, 0}}, 
      ElementMarker == regs["fluid"]}, {{{0, 0}}, True}}]];
fac = Evaluate[
   Piecewise[{{0.2, ElementMarker == regs["solid"]}, {0, True}}]];
k = Evaluate[
   Piecewise[{{285, ElementMarker == regs["solid"]}, {1, True}}]];
op = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[t, x, y]\)\) + 
   v.Inactive[Grad][u[t, x, y], {x, y}] - 
   Inactive[
     Div][(-{{k, 0}, {0, k}}.Inactive[Grad][u[t, x, y], {x, y}]), {x, 
     y}] - fac q[t];
Subscript[\[CapitalGamma], D1] = 
  DirichletCondition[u[t, x, y] == 0, 
   ElementMarker == bounds["inlet"]];
Subscript[\[CapitalGamma], D2] = 
  DirichletCondition[u[t, x, y] == 0, ElementMarker == bounds["cold"]];
ufunHeat = 
  NDSolveValue[{op == 0, u[0, x, y] == 0 , Subscript[\[CapitalGamma], 
    D1], Subscript[\[CapitalGamma], D2], q[0] == 0, 
    WhenEvent[{Mod[t, period], 
      Mod[t + period duty, period]}, {q[t] -> 
       If[q[t] == 0, 1/duty, 0]}, 
     "DetectionMethod" -> "Interpolation"]}, {u, q}, {t, 0, 
    5}, {x, y} \[Element] mesh, DiscreteVariables -> q, 
   MaxStepFraction -> 0.001];
plrng = {{lf, rt}, {bt, tp}, {0, 0.22}};
(* Movie Generation *)
SetOptions[Plot3D, PlotRange -> plrng, PlotPoints -> Full, 
  ColorFunction -> 
   Function[{x, y, z}, Directive[ColorData["DarkBands"][5 z]]], 
  ColorFunctionScaling -> False, MeshFunctions -> {#3 &}, Mesh -> 20, 
  AxesLabel -> Automatic, ImageSize -> Large];
plts = Plot3D[ufunHeat[[1]][#, x, y], {x, y} \[Element] m1, 
    MeshStyle -> {Black, Thick}] &;
pltf = Plot3D[ufunHeat[[1]][#, x, y], {x, y} \[Element] m2, 
    MeshStyle -> {Dashed, Black, Thick}] &;
showplot = 
  Show[{plts[#], pltf[#]}, 
    ViewPoint -> {-2.9775556124522455`, 0.6436172037401853`, 
      1.473064652282362`}, 
    ViewVertical -> {0.4255034386507697`, -0.09197522028503674`, 
      0.9002707273647687`}, ImageSize -> 400, 
    Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False] &;
ListAnimate[showplot /@ Evaluate@Subdivide[0, 5, 80]]

The above code should produce the following animation.I have made no attempts at validation, but the model seems to be behaving reasonably well.

Pulsed Convection

Here is a plot of the temperature taken at the vertical middle and the horizontal beginning, middle, and end of the strip.

Plot[{ufunHeat[[1]][t, 0.05 rt, -th1/2], 
  ufunHeat[[1]][t, 0.5 rt, -th1/2], 
  ufunHeat[[1]][t, 0.95 rt, -th1/2]}, {t, 0, 5}, PlotPoints -> {200}, 
 WorkingPrecision -> 20, MaxRecursion -> 10, PlotRange -> {0, 0.280}, 
 ImageSize -> 600, PlotTheme -> "Web", 
 Filling -> {2 -> {{3}, {LightGreen}}, 1 -> {{2}, {LightYellow}}}, 
 PlotLegends -> 
  Placed[SwatchLegend[{"Beg", "Mid", "End"}, 
    LegendFunction -> "Frame", LegendLayout -> "Column", 
    LegendMarkers -> list[[-1]]], {{0.1, 0.75}, {0.15, 0.75}}]]

Temperature at mid strip

It looks similar to the graph provided in the OP.

I do not precisely know the inner workings of WhenEvent, but other solvers will tighten their time steps around explicit events. I would presume that the same happens in Mathematica. Because it is a physical system with finite diffusivity, the square pulses will most likely be convoluted With a broadening function and will manifest itself as a Gaussian or Lorentzian type shape.

Inlet Boundary Condition Sensitivity

At the liquid-solid inlet interface, the model appears to be pinned. This is due to the Dirichlet condition at the shared node. Local heat transfer coefficients are infinite at the entrance for constant temperature or constant flux prescribed boundary conditions. This pinning would be required if one wanted to compare to analytical solutions. However, in real systems, although local heat transfer coefficients can be very high at the entrance, they are not infinite. Depending on your need, you may want to make adjustments to the inlet boundary condition.

As stated previously, we can override that condition by adjusting the west-solid boundary after the inlet assignment. Alternatively, we can extend the model by adding a solid insulation layer before the heated solid. I also adjusted the equations and the domain a bit, but we still should be able to observe if the model is still pinned at the interface.

Adjusting the Inlet Interface Node to be a Default Insulating Neumann Value

We can adjust the model and simulate with the following code:

(* Model Dimensions *)
th1 = 1;
th2 = 2 th1;
lf = 0;
rt = 5 th1;
bt = -th1;
tp = th2;
(* Horizontal Region *)
rh = pointsToMesh@meshGrowth[lf, rt, 80, 8];
(* Thin Metal Region Uniform Mesh*)
rv = pointsToMesh[Subdivide[bt, 0, 10]];
(* Thick Fluid Region Geometric Growth Mesh *)
rv2 = pointsToMesh@meshGrowth[0, tp, 80, 32];
(* Build Element Meshes From Region Products *)
m1 = rp2Mesh[rh, rv, regs["solid"]];
m2 = rp2Mesh[rh, rv2, regs["fluid"]];
(* Combine the solid and fluid mesh *)
mesh = combineMeshes[m1, m2];
(* Define a series of BC adjustments *)
(* Last assignement takes precedence with PointElement *)
adjustments = {bcadj["solid", "south", "hot"]};
AppendTo[adjustments, bcadj["fluid", "north", "cold"]];
AppendTo[adjustments, bcadj["fluid", "west", "inlet"]];
AppendTo[adjustments, bcadj["solid", "west", "default"]];
(* Adjust the mesh with new boundary and point elements *)
{bcEle, mesh} = adjustMarkers[mesh, adjustments];
(* Display the mesh and bc's *)
Column[{Row@{mesh[
     "Wireframe"["MeshElement" -> "BoundaryElements", 
      "MeshElementMarkerStyle" -> Blue, 
      "MeshElementStyle" -> {Black, Green, Red}, ImageSize -> Medium]],
    mesh["Wireframe"[
      "MeshElementStyle" -> {FaceForm[Red], FaceForm[Green]}, 
      ImageSize -> Medium]]},
  Row@{mesh[
     "Wireframe"["MeshElement" -> "PointElements", 
      "MeshElementIDStyle" -> Black, ImageSize -> Medium]],
    mesh["Wireframe"["MeshElement" -> "PointElements", 
      "MeshElementMarkerStyle" -> Blue, 
      "MeshElementStyle" -> {Black, Green, Red}, 
      ImageSize -> Medium]]}}]
duty = 1/6000 (*6000*);
period = 1;
w = 1/period;
tmax = 10;
v = Evaluate[
   Piecewise[{{{16.6 (y/th2)^2 {1, 0}}, 
      ElementMarker == regs["fluid"]}, {{{0, 0}}, True}}]];
fac = Evaluate[
   Piecewise[{{1, ElementMarker == regs["solid"]}, {0, True}}]];
gamma = Evaluate[
   Piecewise[{{1, ElementMarker == regs["solid"]}, {1.64, True}}]];
k = Evaluate[
   Piecewise[{{0.446, ElementMarker == regs["solid"]}, {50 0.0021, 
      True}}]];
op = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[t, x, y]\)\) + 
   v.Inactive[Grad][u[t, x, y], {x, y}] - 
   Inactive[
     Div][(-{{k, 0}, {0, k}}.Inactive[Grad][u[t, x, y], {x, y}]), {x, 
     y}] - fac q[t];
Subscript[\[CapitalGamma], D1] = 
  DirichletCondition[u[t, x, y] == 0, 
   ElementMarker == bounds["inlet"]];
Subscript[\[CapitalGamma], D2] = 
  DirichletCondition[u[t, x, y] == 0, ElementMarker == bounds["cold"]];
ufunHeat = 
  NDSolveValue[{op == 0, u[0, x, y] == 0 , Subscript[\[CapitalGamma], 
    D1], Subscript[\[CapitalGamma], D2], q[0] == 0, 
    WhenEvent[{Mod[t, period], 
      Mod[t + period duty, period]}, {q[t] -> 
       If[q[t] == 0, 1/duty, 0]}, 
     "DetectionMethod" -> "Interpolation"]}, {u, q}, {t, 0, 
    tmax}, {x, y} \[Element] mesh, DiscreteVariables -> q, 
   MaxStepFraction -> 0.001];

Adjusted Dirichlet Condition

We can see that the pinning effect is reduced, but it has not been eliminated.

Adding an Insulated Entrance Region

Next we will try adding an insulated entrance region before the heated element to mitigate the pinning at the entrance. I added a yellow insulated region to obtain the mesh shown below.

Insulated Layer

Unfortunately, I have reached the character limit, but when the above mesh is simulated, it produces the following:

Padded Simulation

We have substantially mitigated the pinning issue by extending the boundary.

Summary

  • Created a 0D model of pulse heating with WhenEvent that behaves reasonably well.
  • Developed some prototype code that facilitates the building of structured quad meshes and assignment of boundary conditions.
  • Demonstrated that the quad mesh works reasonably well on a two layer heat equation system.
  • Created a prototype with convection that works reasonably well.
  • Made no attempts at validation and the code should be used at your own risk.
Improve this answer
$\endgroup$
6
  • 2
    $\begingroup$ Great answer. You can verify that WhenEvent indeed does tighten the time steps around events by looking at them: ListPlot[Evaluate[(u /. sol)["Coordinates"][[1]]], AxesLabel -> {"Steps", "Time"}] shows that when the events happen more time steps are taken. $\endgroup$
    – user21
    Sep 9, 2019 at 6:44
  • $\begingroup$ @user21 Thank you! I thought there should be a way to show. Also, I think this exercise illustrates the point we discussed about point elements here. For certain classes of problems (e.g., conjugate heat transfer), you probably need a special BC beyond just DCs or NVs, but those too can have their quirks. I found that extending the model is often easier than figuring out the quirks. $\endgroup$
    – Tim Laska
    Sep 9, 2019 at 13:35
  • $\begingroup$ @TimLaska. Great answer indeed! The key is in creating a state variable (q[t]) that is a discrete variable used to turn on and off the heat pulse. I'll need weeks to parse through the rest. Thank you! $\endgroup$
    – BohemianTapestry
    Sep 10, 2019 at 18:17
  • $\begingroup$ Tim, would you mind if we used this example in the documentation? $\endgroup$
    – user21
    Mar 10, 2020 at 9:04
  • 1
    $\begingroup$ @BohemianTapestry, Tim made a proposal to add anisotropic meshing capabilities to the FEM mesh generation here. If you think this is worthwhile, consider giving this proposal an upvote $\endgroup$
    – user21
    Feb 22, 2021 at 10:34

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.