3D_Surf_Constraint

{  3D_SURF_CONSTRAINT.PDE

This problem demonstrates the use of constraints on surface integrals in 3D.

( This is a modification of problem  3D_BRICKS.PDE)

}

title '3D Surface Constraint'

select

   regrid=off        { use fixed grid to speed up demonstration }

coordinates

   cartesian3

variables

   Tp

definitions

   long = 1

   wide = 1

   K                           { thermal conductivity -- values supplied later }

   Q = 10*exp(-x^2-y^2-z^2)              { Thermal source }

initial values

   Tp = 0.

equations

   div[k*grad(Tp)] + Q = 0        { the heat equation }

constraints

   sintegral(normal(k*grad(Tp))) = 1        { total surface integral  }

extrusion

surface "bottom" z = -long

layer 'bottom'

surface "middle" z=0

layer 'top'

surface 'top' z= long        { divide Z into two layers }

boundaries

   surface 1 value(Tp)=0        { fix bottom surface temp }

   surface 3 value(Tp)=0        { fix top surface temp }

   Region 1                { define full domain boundary in base plane }

      layer 1 k=1                { bottom right brick }

      layer 2 k=0.1                { top right brick }

      start(-wide,-wide)

        value(Tp) = 0                { fix all side temps }

        line to (wide,-wide)        { walk outer boundary in base plane }

          to (wide,wide)

          to (-wide,wide)

          to close

   Region 2 "Left"                { overlay a second region in left half }

      layer 1 k=0.2                { bottom left brick }

      layer 2 k=0.4                 { top left brick }

      start(-wide,-wide)

        line to (0,-wide)                { walk left half boundary in base plane }

          to (0,wide)

          to (-wide,wide)

          to close

monitors

     contour(Tp) on surface z=0  as "XY Temp"

     contour(Tp) on surface x=0  as "YZ Temp"

     contour(Tp) on surface y=0  as "ZX Temp"

     elevation(Tp) from (-wide,0,0) to (wide,0,0)  as "X-Axis Temp"

     elevation(Tp) from (0,-wide,0) to (0,wide,0)  as "Y-Axis Temp"

     elevation(Tp) from (0,0,-long) to (0,0,long)  as "Z-Axis Temp"

plots

     contour(Tp) on surface z=0  as "XY Temp"

     contour(Tp) on surface x=0  as "YZ Temp"

     contour(Tp) on surface y=0  as "ZX Temp"

     summary

       report("Constraint Validation:")

       report(sintegral(normal(k*grad(Tp))))        { surface integral on total outer surface }

end