<< Click to Display Table of Contents >> complex_sinusoidal_heat |
{ COMPLEX_SINUSOIDAL_HEAT.PDE
This example demonstrates the use of COMPLEX variables and ARRAY definitions
to compute the time-sinusoidal behavior of a rod in a box.
The heat equation is
div(k*grad(temp)) = cp*dt(temp)
If we assume that the sources and solutions are in steady oscillation at a frequency
omega, then we can write
temp(x,y,t) = phi(x,y)*exp(i*omega*t) = phi(x,y)*(cos(omega*t) + i*sin(omega*t))
Substituting this into the heat equation and dividing the exp(i*omega*t) out of the
result leaves
div(k*grad(phi)) - i*omega*cp*phi = 0
The temperature temp(x,y,t) can be reconstructed at any time by expanding the above
definition.
In this example, we construct an array of sample times and the associated arrays
of sine and cosine factors. These arrays are then used to display a time history of
temperature at various points in the domain.
}
TITLE 'Time Sinusoidal Heat flow around an Insulating blob '
VARIABLES
! define the complex amplitude function phi and its real and imaginary components
phi=complex(phir,phii)
DEFINITIONS
k=1
ts = array (0 by pi/20 to 2*pi) ! an array of sample times
fr = cos(ts) ! sine and cosine arrays
fi = sin(ts)
! define a function for evaluating a time array of temp(px,py,t) at a point
temp(px, py) = eval(phir,px,py)*fr + eval(phii,px,py)*fi
EQUATIONS
phi: Div(k*grad(phi)) - complex(0,1)*phi= 0
BOUNDARIES
REGION 1 'box'
START(-1,-1)
VALUE(Phi)=complex(0,0) LINE TO(1,-1) { Phi=0 in base line }
NATURAL(Phi)=complex(0,0) LINE TO (1,1) { normal derivative =0 on right side }
VALUE(Phi)=complex(1,0) LINE TO (-1,1) { Phi = 1 on top }
NATURAL(Phi)=complex(0,0) LINE TO CLOSE { normal derivative =0 on left side }
REGION 2 'rod' { the embedded circular rod }
k=0.01
START 'ring' (1/2,0)
ARC(CENTER=0,0) ANGLE=360 TO FINISH
PLOTS
CONTOUR(Phir) ! plot the real part of phi
REPORT(k) REPORT(INTEGRAL(Phir, 'rod'))
CONTOUR(Phii) ! plot the imaginary part of phi
REPORT(k) REPORT(INTEGRAL(Phii, 'rod'))
! reconstruct the temperature distribution at a few selected times
REPEAT tx=0 by pi/2 to 2*pi
SURFACE(phir*cos(tx)+phii*sin(tx)) as "Phi at t="+$[4]tx
ENDREPEAT
! plot the time history at a few selected positions
ELEVATION(temp(0,0), temp(0,0.2), temp(0,0.4), temp(0,0.5)) vs ts as "Phi(t) at X=0, Y=(0, 0.2 ,0.4, 0.5)"
VECTOR(-k*grad(Phir))
! plot a lineout of phir and phii through the domain
ELEVATION(Phi) FROM (0,-1) to (0,1)
! plot the real component of flux on the surface of the rod
ELEVATION(Normal(-k*grad(Phir))) ON 'ring'
END