| Author | Message | ||
     Tim McDaniel (tim) New member Username: tim Post Number: 1 Registered: 04-2007 |
Hi, I am trying to model the thermo-elastic PDE. The situation I am trying to model is a 500 ps laser pulse with a spot size of 1e-6 m incident in a piece of glass. The laser pulse serves as the heat source and I am concerned with the behavior of the density waves on the time scale of picoseconds. The problem is spherically symmetric so I am using the sphere1 coordinate system. My unit system is in meters and seconds. The problem that I am having is that there seems to be a discontinuity in the temperature while running the script. I am looking at a time range from 0 to 2 nanoseconds. The temperature rises to the max temperature and is steady state (what I expect) until 1.5 ns, when it suddenly drops about 100 K. Also, the max temperature that glass rises to is about 200K more than what it should given the temperature equation. What is going on here? Any help would be much appreciated. Here is the script: COORDINATES sphere1 VARIABLES u (threshold = 2e-9) {displacement, (t,r) } temp (threshold = 1300) {temperature distribution, (t,r)} p {density, (t,r) } v (threshold = 0.1) { v = dt(u) } p_ratio DEFINITIONS p0 = 2.51e3 { initial density, kg/m^3 } Y = 72.9e9 { Young's Modulus, Pa } s = 0.21 { Poisson's ratio } beta = 72e-7 { themal expansion coefficient, K^-1 } tempi = 300 { initial temp, K } tempd = 1000 { temp increase, K } w = 1e-6 { spatial width of temperature distribution, m } tau = 0.5e-12 { pulse width of laser, s } A = ( Y * (1 - s) ) / ( (1 + s) * (1 - 2*s) ) B = Y / ( (1 + s) * (1 - 2*s) ) C = ( Y * beta ) / ( 3 * (1 - 2*s) ) Ra = 0 Rb =1e-4 p_ratio_percent = p_ratio * 100 INITIAL VALUES u = 0 temp = tempi p = p0 EQUATIONS v: dt(u) = v u: p*dt(v) = ( A * ( drr(u) + (1/r)*( dr(u) ) ) - B * (u/((r^2))) - C * dr(temp) ) p: p = p0 / ( dr(u) + 1 ) temp: temp = tempi + ( tempd ) *( EXP ( - (r^2 / w^2 ) ) ) * (1 - EXP (-( t / tau ) ) ) p_ratio: p_ratio = (p - p0) / p0 BOUNDARIES REGION 1 START(Ra) POINT VALUE(u) = 0 LINE TO (Rb) POINT LOAD(p_ratio)= 0 TIME 0 BY 1e-12 TO 2e-9 MONITORS FOR TIME 0 BY 1e-12 TO 2e-9 ELEVATION (temp) FROM (Ra) to (6e-6) RANGE(0,1400) ELEVATION (u) FROM (Ra) TO (30e-6) RANGE(0,1e-9) ELEVATION (p) FROM (Ra) TO (12e-6) ELEVATION(p_ratio) FROM (Ra) TO (12e-6) RANGE(-4e-3,1e-3) HISTORY(temp) AT (0) HISTORY(p_ratio_percent) AT (0) RANGE(-0.4,0.1) PLOTS FOR TIME 0 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 0' RANGE(-4e-3,1e-3) FOR TIME 100e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 100 ps' RANGE(-4e-3,1e-3) FOR TIME 200e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 200 ps' RANGE(-4e-3,1e-3) FOR TIME 400e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 400 ps' RANGE(-4e-3,1e-3) FOR TIME 600e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 600 ps' RANGE(-4e-3,1e-3) FOR TIME 800e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 800 ps' RANGE(-4e-3,1e-3) FOR TIME 1200e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 1200 ps' RANGE(-4e-3,1e-3) FOR TIME 1600e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 1600 ps' RANGE(-4e-3,1e-3) FOR TIME 2000e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 2000 ps' RANGE(-4e-3,1e-3) END Thanks, Tim | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 832 Registered: 06-2003 |
I'm not sure what causes the temperature jump, but TEMP should not be a variable anyway. It is a simple algebraic function of time and space, and should therefore be a definition. P and P_RATIO are also algebraic functions, and should be definitions. This leaves you two equations instead of five. In the modified script below, I have made these changes, and also re-labeled the equations so that the label follows the time derivative (I don't know if this matters). I have also applied the boundary condition POINT VALUE(v)=0, since v=dt(u) and u=0 on the left boundary. Without this there is instability at the left boundary. I have changed the MONITOR control to FOR CYCLE=1, so it doesn't limit the timestep. The new script follows: COORDINATES sphere1 VARIABLES u (threshold = 1e-11) {displacement, (t,r) } ! temp (threshold = 1300) {temperature distribution, (t,r)} ! p {density, (t,r) } v (threshold = 0.1) { v = dt(u) } ! p_ratio DEFINITIONS p0 = 2.51e3 { initial density, kg/m^3 } Y = 72.9e9 { Young's Modulus, Pa } s = 0.21 { Poisson's ratio } beta = 72e-7 { themal expansion coefficient, K^-1 } tempi = 300 { initial temp, K } tempd = 1000 { temp increase, K } w = 1e-6 { spatial width of temperature distribution, m } tau = 0.5e-12 { pulse width of laser, s } A = ( Y * (1 - s) ) / ( (1 + s) * (1 - 2*s) ) B = Y / ( (1 + s) * (1 - 2*s) ) C = ( Y * beta ) / ( 3 * (1 - 2*s) ) Ra = 0 Rb =1e-4 p = p0 / ( dr(u) + 1 ) p_ratio = (p - p0) / p0 p_ratio_percent = p_ratio * 100 temp = tempi + ( tempd ) *( EXP ( - (r^2 / w^2 ) ) ) * (1 - EXP(-( t / tau ) ) ) INITIAL VALUES u = 0 ! temp = tempi ! p = p0 EQUATIONS ! v: dt(u) = v u: dt(u) = v ! u: p*dt(v) = ( A * ( drr(u) + (1/r)*( dr(u) ) ) - B * (u/((r^2))) - C* dr(temp) ) v: p*dt(v) = ( A * ( drr(u) + (1/r)*( dr(u) ) ) - B * (u/((r^2))) - C* dr(temp) ) ! p: p = p0 / ( dr(u) + 1 ) ! temp: temp = tempi + ( tempd ) *( EXP ( - (r^2 / w^2 ) ) ) * (1 - EXP(-( t / tau ) ) ) ! p_ratio: p_ratio = (p - p0) / p0 BOUNDARIES REGION 1 START(Ra) POINT VALUE(u) = 0 POINT VALUE(v) = 0 LINE TO (Rb) ! POINT LOAD(p_ratio)= 0 TIME 0 BY 1e-12 TO 2e-9 MONITORS ! FOR TIME 0 BY 1e-12 TO 2e-9 FOR cycle=1 ELEVATION (temp) FROM (Ra) to (6e-6) RANGE(0,1400) ELEVATION (u) FROM (Ra) TO (30e-6) RANGE(0,1e-9) ELEVATION (v) FROM (Ra) TO (30e-6) ELEVATION (p) FROM (Ra) TO (12e-6) ELEVATION(p_ratio) FROM (Ra) TO (12e-6) RANGE(-4e-3,1e-3) HISTORY(temp) AT (0) HISTORY(p_ratio_percent) AT (0) RANGE(-0.4,0.1) PLOTS FOR TIME 0 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 0' RANGE(-4e-3,1e-3) FOR TIME 100e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 100 ps' RANGE(-4e-3,1e-3) FOR TIME 200e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 200 ps' RANGE(-4e-3,1e-3) FOR TIME 400e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 400 ps' RANGE(-4e-3,1e-3) FOR TIME 600e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 600 ps' RANGE(-4e-3,1e-3) FOR TIME 800e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 800 ps' RANGE(-4e-3,1e-3) FOR TIME 1200e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 1200 ps' RANGE(-4e-3,1e-3) FOR TIME 1600e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 1600 ps' RANGE(-4e-3,1e-3) FOR TIME 2000e-12 ELEVATION(p_ratio) FROM (Ra) TO (12e-6) AS 'p_ratio at t = 2000 ps' RANGE(-4e-3,1e-3) END |