I'm trying to solve Poisson's equation for the temperature difference ($\Delta T$) between two layers with known conduction between them. The layers have conductive elements spaced in a hexagonal pattern, so I'm solving the equation over a hexagonal unit cell with periodic boundary conditions on all sides. I set a simple Dirichlet condition that $\Delta T=0$ at the center.
(*A function to define the hexagonal domain of the calculation, with a scaling parameter.*)
\[CapitalOmega]Hex[scale_] := Polygon[{{1, 0}, {1/2, Sqrt[3]/2}, {-1/2, Sqrt[3]/2}, {-1, 0}, {-1/2, -Sqrt[3]/2}, {1/2, -Sqrt[3]/2}}*scale]
Needs["NDSolve`FEM`"]
(*Now the main function.*)
TempSolver[scale_, rpost_] := (
(*First define the vertical conduction through the air and the solid posts.*)
parameters = {cond1 -> 5000, cond2 -> 15000, extraHeat -> 100, condxy -> 10^-7};
(*Define where each conduction term applies in the domain.*)
cond =
Piecewise[{{cond1,
Sqrt[x^2 + y^2] > rpost && Sqrt[(x - scale)^2 + y^2] > rpost && Sqrt[(x + scale)^2 + y^2] > rpost &&
Sqrt[(x - 0.5*scale)^2 + (y - Sqrt[3]/2*scale)^2] > rpost && Sqrt[(x - 0.5*scale)^2 + (y + Sqrt[3]/2*scale)^2] > rpost &&
Sqrt[(x + 0.5*scale)^2 + (y - Sqrt[3]/2*scale)^2] > rpost && Sqrt[(x + 0.5*scale)^2 + (y + Sqrt[3]/2*scale)^2] > rpost}},
cond2];
mesh = ToElementMesh[\[CapitalOmega]Hex[scale], "IncludePoints" -> {{0, 0}}];
eps = 10^-8;
pbcs = {PeriodicBoundaryCondition[\[CapitalDelta]T[x, y], Abs[y/(x + scale) - Sqrt[3]]<= eps, TranslationTransform[{3/2 scale, -Sqrt[3]/2*scale}]],(*top left to bottom right*)
PeriodicBoundaryCondition[\[CapitalDelta]T[x, y], Abs[y + Sqrt[3]/2*scale] <= eps, TranslationTransform[{0, Sqrt[3]*scale}]],(*bottom to top*)
PeriodicBoundaryCondition[\[CapitalDelta]T[x, y], Abs[y/(x + scale) + Sqrt[3]] <= eps, TranslationTransform[{3/2 scale, Sqrt[3]/2*scale}]](*bottom left to top right*)};
NDSolve[{extraHeat + condxy*Laplacian[\[CapitalDelta]T[x, y], {x, y}] - 2*\[CapitalDelta]T[x, y]*cond == 0, pbcs, DirichletCondition[\[CapitalDelta]T[x, y] == 0, x == 0 && y == 0]} /. parameters, \[CapitalDelta]T, {x, y} \[Element] mesh])
(*A function to plot the result*)
TempPlot[scale_, rpost_] := (soln = TempSolver[scale, rpost];
Show[DensityPlot[Evaluate[(\[CapitalDelta]T[x/10^6, y/10^6] /. soln[[1]])], {x, -scale*10^6, scale*10^6}, {y, -scale*10^6, scale*10^6}, PlotRange -> All, PlotLegends -> Automatic, Frame -> True, GridLines -> Automatic, ColorFunction -> "TemperatureMap", FrameStyle -> Directive[Black, 12], FrameLabel -> {Style["x (\[Mu]m)", 12, Black], Style["y (\[Mu]m)", 12, Black], Style["\[CapitalDelta]T", 12, Black]}]])
This returns several errors, such as:
NDSolve::fembcbadpred: Applying the predicate Abs[-Sqrt[3]+y/(1/10000+x)]<=1/100000000 from boundary condition PeriodicBoundaryCondition[\[CapitalDelta]T,Abs[-Sqrt[3]+y/(1/10000+x)]<=1/100000000,TransformationFunction[{{1,0,3/20000},{0,1,-(Sqrt[3]/20000)},{0,0,1}}]] on the discretized region failed. The predicate must evaluate to machine precision numbers at all coordinates of the discretized boundary. As an alternative, try to specify the boundary predicate with an ElementMarker.
ReplaceAll::rmix: Elements of {220-2 (\[Piecewise] 5000 And[<<7>>] 15000 True) \[CapitalDelta]T[x,y]+((\[CapitalDelta]T^(0,2))[x,y]+(\[CapitalDelta]T^(2,0))[x,y])/5000000==0,{PeriodicBoundaryCondition[\[CapitalDelta]T[x,y],Abs[-Power[<<2>>]+Power[<<2>>] y]<=1/100000000,TransformationFunction[{{1,0,3/20000},{0,1,-(Power[<<2>>]/20000)},{0,0,1}}]],PeriodicBoundaryCondition[\[CapitalDelta]T[x,y],Abs[Power[<<2>>]/20000+y]<=1/100000000,TransformationFunction[{{1,0,0},{0,1,Power[<<2>>]/10000},{0,0,1}}]],PeriodicBoundaryCondition[\[CapitalDelta]T[x,y],Abs[Sqrt[3]+Power[<<2>>] y]<=1/100000000,TransformationFunction[{{1,0,3/20000},{0,1,Power[<<2>>]/20000},{0,0,1}}]]},DirichletCondition[\[CapitalDelta]T[x,y]==0,x==0&&y==0]} are a mixture of lists and nonlists.
This worked on a previous version of Mathematica (now using 14.1). Any ideas? Also, is there a way to simplify my periodic boundary conditions?
scale, rpost? $\endgroup$