The Gold Standard of Free-Boundary Problems

Chapter index

In an ordinary fixed-domain heat problem, the rod or region is supplied in advance. You solve for temperature inside it. Boundary temperatures or heat fluxes are imposed on known locations. In a melting problem, the material may occupy a fixed container, yet the region occupied by each phase changes. The container boundary is known; the phase boundary is not.

Write the moving interface in one dimension as x=s(t). The temperature may satisfy one equation for 0<x<s(t) and another for s(t)<x<a. The function s is an additional unknown. Its location affects the temperature gradients, while those gradients affect its speed. This feedback is what makes the problem coupled.

A free boundary need not move in time. In a static contact problem, the region where a membrane touches an obstacle is also unknown. Its edge separates two mathematical regimes, even though the final configuration is stationary. Conversely, a moving boundary is not necessarily free: a piston whose position is prescribed by an external motor has known motion in the mathematical model.

Three questions that expose the distinction

Ask whether the geometry is supplied, inferred from the solution or prescribed as a function of time. Then ask which condition determines it. Finally ask which quantities are solved on either side. These questions prevent the phrase “moving boundary” from hiding important differences between models.

For example, prescribing s(t)=√t and solving a heat equation inside that region is a moving-domain calculation. It becomes a Stefan free-boundary problem only when s must also satisfy an energy balance derived from phase change. A visually plausible front can fail that balance and therefore fail the physical model.

The fixed-domain equations remain essential. The parabolic PDE guide explains diffusion and heat evolution. Free-boundary reasoning adds the missing geometry and interface law. It does not remove the need for initial data, external boundary conditions or a clear definition of the unknown field.

Contents · Next chapter

Consider melting from the left. Liquid occupies 0<x<s(t), solid lies to the right, and positive x points into the solid. Assume equal constant densities, a sharp interface, negligible fluid motion and an equilibrium melting temperature Tm. Let q be heat flux measured positive to the right, with Fourier’s law q=−kT_x.

During a short time Δt, an interface of unit area advances by approximately s′Δt. The mass melted is ρs′Δt, so the latent energy required is ρLs′Δt. The incoming heat from the liquid side is q_leftΔt and the outgoing heat into the solid side is q_rightΔt. Their difference supplies the phase change.

The resulting condition is ρLs′=q_left−q_right. Equivalently, ρLs′=k_s T_x(s+,t)−k_l T_x(s−,t). The signs follow from the chosen coordinate direction. Rather than memorising a jump formula, redraw the interface and identify energy entering and leaving it.

A numerical flux balance

Suppose q_left=900 W/m² and q_right=300 W/m². Let ρ=1,000 kg/m³ and L=200,000 J/kg. Then s′=600/200,000,000=3×10⁻⁶ m/s, or 10.8 mm per hour if those fluxes remain constant. This is a local instantaneous prediction; it does not establish that the fluxes will stay fixed as the region grows.

If the outgoing flux rises to 900 W/m², the interface is stationary under this balance. If it exceeds the incoming flux, s′ becomes negative and the assumed interface orientation describes retreat. A solver that clips all negative velocities to zero would change the model unless an independent physical constraint justified that operation.

The dimensional check is useful: ρL has units J/m³ and s′ has units m/s, so their product is W/m². If a calculation uses total heating power in watts on the right, divide by area first or write the corresponding total-energy equation consistently.

Previous chapter · Contents · Next chapter

A useful idealisation keeps the solid at Tm and solves temperature only in the liquid. Heat conduction in the liquid satisfies ρcT_t=kT_xx, or T_t=αT_xx with α=k/(ρc). At x=0 impose Th>Tm. At the unknown front impose T(s(t),t)=Tm and ρLs′=−kT_x(s−,t).

The last expression is the preceding energy balance with zero conductive flux on the solid side. It is not an arbitrary front-speed rule. The temperature gradient determines how much heat reaches the front, and latent heat determines how much new liquid that heat can produce.

For the canonical semi-infinite problem, material begins at Tm with no liquid layer, so s(0)=0. That ideal start introduces a singular initial gradient under a suddenly imposed boundary temperature. It is mathematically useful, but it does not reproduce every detail of switching on a real heater.

What the model leaves out

The assumptions exclude convection in the melt, changing density, variable material properties, a melting range, external heat-transfer resistance and interfacial kinetic effects. A pure substance under suitable conditions may make a sharp equilibrium interface a useful model. Mixtures or alloys can require a region containing both phases rather than one sharply located boundary.

A two-phase model also solves heat transfer in the solid. If that solid begins below Tm, some supplied energy warms it before or during melting. Treating its flux as zero can overpredict interface advance. The correct simplification depends on the scale of that contribution, not on the convenience of having one equation instead of two.

Oxford’s continuum-methods notes formulate the Stefan problem through this combination of heat equations and latent-energy balance. Use that formulation as a starting point, then state explicitly which terms your own application retains. A list of constants without a statement of the model assumptions does not define a complete problem.

Previous chapter · Contents · Next chapter

Suppose the liquid temperature adjusts rapidly compared with the advance of the front. Approximate its profile as linear: T(x,t)=Th−ΔT x/s(t), where ΔT=Th−Tm. The heat flux reaching the front is then kΔT/s. Substitution into the Stefan condition gives s′=kΔT/(ρLs).

Multiply by s and integrate: s²(t)=s₀²+2kΔTt/(ρL), where s₀ is the initial liquid-layer thickness. With s₀=0, the front grows like √t. The reason is mechanistic: a thicker layer offers more conductive resistance, so the front slows as it moves farther from the heated surface.

Original calculation

Take k=2 W/(m K), ΔT=10 K, ρ=1,000 kg/m³ and L=200,000 J/kg. Then 2kΔT/(ρL)=2×10⁻⁷ m²/s. Starting from zero thickness, the approximate melted depth at t=500 s is √(10⁻⁴)=0.01 m, or 10 mm.

At t=2,000 s, four times the duration, the depth is 20 mm rather than 40 mm. The instantaneous speed at 500 s is kΔT/(ρLs)=10⁻⁵ m/s. The average speed over the first 500 s is 2×10⁻⁵ m/s, twice as large, because the front was faster earlier in this idealised calculation.

Why this is not the exact transient temperature solution

The linear profile has T_xx=0, but its time derivative is generally nonzero because s changes. It therefore does not satisfy the full heat equation exactly. It is a quasi-steady approximation. Calling it exact would hide the neglected sensible heating of the liquid.

A useful scale ratio is the Stefan number Ste=cΔT/L. Small Ste supports a regime in which latent energy dominates sensible energy, making this approximation more plausible. It does not excuse ignoring the initial singularity, finite boundaries or other physics. Use the ratio to motivate an approximation, then test the error against a better model or benchmark.

Previous chapter · Contents · Next chapter

The canonical one-phase melting problem admits a similarity form s(t)=2λ√(αt). The temperature profile is T=Th−ΔT erf[x/(2√(αt))]/erf(λ), where erf is the error function. The constant λ is determined by the interface energy balance rather than chosen to fit a desired front location.

Differentiating the profile at x=s gives −kT_x=kΔT exp(−λ²)/[√(παt) erf(λ)]. Since s′=λ√α/√t, the Stefan condition reduces to √π λ exp(λ²) erf(λ)=Ste. This scalar equation is the remaining nonlinear task after the similarity reduction.

A controlled numerical comparison

Keep the previous k, ρ, L and ΔT, and choose c=2,000 J/(kg K). Then α=10⁻⁶ m²/s and Ste=0.1. Solving the scalar equation gives λ approximately 0.2200. At t=500 s, the front depth is approximately 2×0.2200×√(0.0005)=0.00984 m.

The quasi-steady prediction was 0.01000 m. It is slightly larger because it directs all the conducted heat in that approximation toward latent melting, whereas the transient model also accounts for sensible heating in the liquid. The difference here is about 0.16 mm, or roughly 1.6 percent relative to the transient depth.

For small λ, erf(λ) is approximately 2λ/√π and exp(λ²) is approximately one. The scalar equation then gives 2λ²≈Ste. Substituting into s²=4λ²αt recovers s²≈2kΔTt/(ρL). This connects the approximate law to a controlled limiting calculation rather than a coincidence.

Addison, Howison and King’s study includes the classical similarity solution within a broader analysis of free-boundary asymptotics. Different limits of the Stefan number produce different balances. Always check which parameter is small before importing an approximation from another regime.

Previous chapter · Contents · Next chapter

A numerical mesh can follow the liquid interval by introducing y=x/s(t), so 0≤y≤1. Define θ(y,t)=T(s(t)y,t). The temperature now lives on a fixed computational interval, but the transformed equation acquires terms describing mesh motion.

At fixed physical position x, y_t=−ys′/s. Therefore T_t=θ_t−(ys′/s)θ_y, while T_x=θ_y/s and T_xx=θ_yy/s². The heat equation becomes θ_t−(ys′/s)θ_y=(α/s²)θ_yy. The front condition becomes ρLs′=−(k/s)θ_y(1,t).

The term involving s′ is sometimes mistaken for physical fluid advection. Here it comes from the moving coordinate system. Omitting it does not merely choose a simpler mesh; it changes the differential equation being solved.

A chain-rule test

Take an artificial field T(x,t)=x and an artificial boundary s(t)=1+t. Then θ(y,t)=(1+t)y. At fixed y, θ_t=y and θ_y=1+t. The combination θ_t−ys′θ_y/s equals y−y=0, matching T_t=0 at fixed x.

This manufactured test is deliberately not a melting solution. Its purpose is to isolate the coordinate transformation. A code can pass a complicated physical test accidentally through compensating errors; a small identity like this makes a missing mesh term unmistakable.

The transformation also exposes a difficulty near s=0, where coefficients involving 1/s and 1/s² become singular. Starting directly at zero with an ordinary finite-difference update can be ill-conditioned. One option is to initialise at a small positive time using a known similarity solution. Another is a formulation that avoids a vanishing computational interval. Either choice requires documenting its effect on the early-time result.

Previous chapter · Contents · Next chapter

An alternative is to keep the whole material on a fixed spatial grid and evolve enthalpy. Enthalpy includes energy stored as sensible heat and latent phase change. For a simplified equal-density material, a volumetric expression is h=ρc(T−Tm)+ρLf, where f is the liquid fraction, with the exact relation between T and f specified by the phase model.

Energy conservation can be written h_t+∇·q=0 when there are no sources or advective contributions. The difficult part shifts from explicitly moving an interface to defining and inverting the relation between enthalpy, temperature and phase fraction. For a sharp isothermal transition, added energy at Tm changes f rather than temperature.

A cell-energy calculation

Consider a cell of volume 10⁻⁶ m³, density 1,000 kg/m³ and latent heat 200,000 J/kg. Its mass is 0.001 kg, so fully melting it at Tm requires 200 J. If it starts 20 percent liquid at Tm and receives a net 50 J, the liquid fraction rises by 50/200=0.25 to 0.45, assuming no other energy change.

A temperature-only update would miss this plateau if it converted all 50 J into sensible heating. Conversely, if the cell begins below Tm, the energy needed to warm it must be supplied before using the entire increment for melting. The order and accounting follow from the enthalpy relation.

Enthalpy methods can handle complex phase regions without explicitly tracking each front. However, the interface may be spread over several cells numerically, and a regularised melting interval introduces another scale. The reported front position must then be defined: for example, by a chosen phase-fraction contour, with sensitivity to that choice assessed.

Conservation of total enthalpy is a valuable check, but it is not sufficient by itself. A method can conserve energy while placing the phase boundary inaccurately because of excessive numerical smearing. Compare both the energy budget and the interface location against a known test problem.

Previous chapter · Contents · Next chapter

Imagine a taut membrane that cannot pass below a rigid obstacle. Its shape minimises an energy subject to the inequality u≥ψ, where ψ is the obstacle. The contact set is where u=ψ; the noncontact set is where u>ψ. Their boundary inside the domain is the free boundary.

For a simple model with energy J(u)=∫[|∇u|²/2−fu] dx and fixed outer boundary values, unconstrained minimisation gives −Δu=f. With the obstacle, the solution instead satisfies u−ψ≥0, −Δu−f≥0 and (u−ψ)(−Δu−f)=0 in an appropriate weak sense. The extra nonnegative quantity acts as a contact reaction.

Why complementarity matters

Where the membrane is separated from the obstacle, u−ψ>0, so the reaction must vanish and the ordinary equation applies. Where it touches, the reaction may be positive to prevent penetration. Solving the unconstrained equation everywhere and clipping the answer afterward generally does not reproduce this coupled equilibrium.

The geometry is determined by the constraint and the equilibrium together. There is no separate clock or interface velocity in this static problem. This makes it an important counterexample to the assumption that every free boundary is a moving material surface.

The exact sign of the reaction inequality follows from the energy convention. If a source or load is defined with the opposite sign, the written equations change. Derive the first variation or test a simple example before comparing formulas from different references.

Xavier Ros-Oton’s overview explain obstacle problems and the regularity of their free boundaries. That theory studies when the boundary is smooth and how singular points can occur. It does not justify assuming that every computed contact edge is smooth merely because the obstacle itself is smooth.

Previous chapter · Contents · Next chapter

Use a one-dimensional model on −2≤x≤2. Set u(−2)=u(2)=1, impose the flat obstacle u≥0 and minimise J(u)=∫[(u′)²/2+2u] dx. The positive linear term favours smaller values, while the derivative term penalises steep changes. Without the obstacle, the equilibrium would satisfy u″=2.

By symmetry, seek a central contact interval [−a,a]. In contact, u=0. On the right noncontact interval, solve u″=2 with u(a)=0 and u′(a)=0. The result is u(x)=(x−a)². The endpoint condition u(2)=1 gives a=1.

Thus the solution is u=0 for |x|≤1 and u=(|x|−1)² for 1<|x|≤2. The two free-boundary points are x=−1 and x=1. They were determined from the equations and constraints, not supplied in advance.

Verify every part

The function is nonnegative and has the prescribed endpoint values. Its first derivative joins continuously to zero at the contact edges. In the noncontact region, −u″+2=0. In the contact interior, −u″+2=2≥0. The product of this reaction and u is zero everywhere apart from the usual weak interpretation at transitions.

The energy is also calculable. On the right, let z=x−1, so 0≤z≤1. The integrand is (2z)²/2+2z²=4z². Its integral is 4/3. Doubling for the left half gives J=8/3. The constant admissible competitor u=1 has energy eight, so it is worse.

The matching slope is not decorative. Joining the free parabola to the obstacle with a slope jump would introduce a concentrated second-derivative contribution. The equilibrium and admissible reaction would then require separate examination. In this simple setting, tangential contact produces a valid weak solution and gives a clean benchmark for an obstacle solver.

Prove optimality instead of checking only one competitor

The energy comparison with a constant function is reassuring, but it does not by itself prove that our candidate beats every admissible function. Let v be any other admissible function with the same endpoint values and write w=v−u. Expanding the squares in the energy gives J(v)−J(u)=∫[u′w′+2w] dx + (1/2)∫(w′)² dx.

Integrating the first term by parts produces ∫(−u″+2)w dx, with no outer endpoint contribution because w vanishes there. There is also no interior slope-jump contribution because our candidate has continuous first derivative. In the noncontact region the coefficient is zero. In the contact region it equals two, while w=v≥0 because u=0 and v must stay above the obstacle.

Both terms in the energy difference are therefore nonnegative. This proves that the candidate minimises the energy over all admissible competitors in the appropriate energy space. Equality forces w′=0; combined with the zero endpoint difference, this gives w=0. The minimiser is unique for this example.

This proof supplies a stronger diagnostic than a plot or a pointwise differential-equation residual. If a numerical candidate penetrates the obstacle, the sign argument fails. If its reaction is negative in contact, the same argument fails. If it has an unexplained slope jump, integration by parts reveals an additional contribution that must be assessed. The calculation connects feasibility, equilibrium and optimality in one place.

It also clarifies what a contact reaction means. It is not an arbitrary extra load chosen to make an answer fit. It is the quantity that enforces the inequality while remaining compatible with minimum energy. That interpretation is useful when extending from a one-dimensional teaching example to a membrane or a more complicated constrained field.

Previous chapter · Contents · Next chapter

A fluid surface carried by the fluid obeys a kinematic condition: its normal speed equals the fluid’s normal velocity when no mass crosses it. A phase boundary can move relative to the material because melting or solidification changes which phase occupies a location. These two mechanisms should not be assigned the same speed law automatically.

If an interface is represented by φ(x,t)=0 and moves with velocity v, differentiation along it gives φ_t+v·∇φ=0. If only normal speed Vn is specified, the corresponding form is φ_t+Vn|∇φ|=0, using a consistent normal orientation. The equation transports geometry; a separate physical law must determine v or Vn.

A volume check for a circular region

Suppose a two-dimensional circular fluid region receives area at a constant rate Q, measured in square metres per second. If area is πR², conservation gives 2πRR′=Q. Hence R²(t)=R₀²+Qt/π. With R₀=1 m and Q=π/2 m²/s, the radius at t=6 s is two metres.

The instantaneous speed then is Q/(2πR)=1/8 m/s. Prescribing a constant radial speed of 1/4 m/s would match the initial speed but not the later area balance. Geometry changes the relationship between volume input and interface speed.

The level-set guide develops implicit interface representation, while the Hamilton–Jacobi root explains its first-order PDE connection. Neither representation supplies missing physics. Before selecting a numerical method, write the condition that determines the normal motion in the actual problem.

Previous chapter · Contents · Next chapter

A fitted mesh places nodes or elements along the interface. It can represent boundary conditions sharply, but large deformations may distort elements and require remeshing. When an interface splits or merges, its connectivity becomes an additional computational challenge.

An implicit level-set field can represent complicated geometry through a zero contour. It handles changes in topology naturally at the representation level, but maintaining geometric accuracy and conserving volume may need additional treatment. Reinitialising the field can improve its distance-like structure while also introducing interface error if done carelessly.

A volume-fraction method tracks how much of each cell is occupied by a phase. Conservation can be built directly into the update, but reconstructing a sharp interface and its curvature adds complexity. An enthalpy method instead focuses on energy and phase fraction, which is particularly natural for melting and solidification.

Coupled structures require more than front tracking

The immersed boundary method couples moving structures to a fluid grid through force spreading and velocity interpolation. The existing immersed-boundary specialist article explains that mechanism. Charles Peskin’s methodological paper develops its mathematical structure for fluid–structure interaction.

Do not choose a method solely because it draws an attractive moving curve. Identify whether the important output is latent-energy consumption, phase volume, contact force, curvature, interface position or a coupled structural response. Different errors dominate those outputs.

For a melting benchmark, compare total supplied energy with sensible and latent storage. For the obstacle benchmark, check nonpenetration, reaction sign and complementarity. For a transported fluid region, inspect volume and boundary position. A single scalar convergence tolerance cannot replace these problem-specific checks.

Previous chapter · Contents · Next chapter

Start with a sign test. Does positive net heat into a melting interface advance it in the intended direction? Does a downward load increase contact rather than lift the membrane away? Does an outward normal speed expand a circle? These tests are simple enough to perform before running a large simulation.

Then check closure. A heat equation plus an interface temperature condition does not determine both temperature and an unknown front without an additional law. Similarly, prescribing both a front path and its speed from a conflicting energy balance can overconstrain a model. Count the unknowns and identify which independent conditions determine them.

Separate four kinds of error

Model error comes from omitted physics or inappropriate assumptions. Discretisation error comes from approximating the equations on a grid or mesh. Algebraic error comes from incompletely solving the discrete equations. Geometric error comes from locating or representing the interface. These categories interact, but reducing one does not automatically reduce the others.

A nonlinear solver may converge tightly while the grid is too coarse to resolve the liquid layer. A fine grid may accurately solve a model that ignores important convection. A good temperature field may coexist with a poor front location if latent heat is mishandled. Report the checks that correspond to the claim being made.

For refinement, compare results at the same physical times and locations, not merely at equal step counts. Reduce the timestep when the spatial mesh is refined if the algorithm requires it. Track whether the front position, conservation residual and relevant field norms stabilise together.

Unexpected roughness also needs interpretation. Some interface models genuinely amplify perturbations; others become ill-posed without missing regularisation. Smoothing every rough edge can conceal the mechanism rather than repair the model. Compare behaviour against the equations, known stability properties and physical scales before treating all irregularity as numerical noise.

Previous chapter · Contents · Next chapter

Return to the obstacle problem, but replace the load coefficient two by a positive parameter q. The energy is J(u)=∫[(u′)²/2+qu] dx on [−2,2], with endpoint values one and obstacle zero. In a noncontact region, u″=q.

If a central contact interval exists, the right-hand solution is u=(q/2)(x−a)². The endpoint condition gives 1=(q/2)(2−a)², so a=2−√(2/q). This is physically admissible as a nonnegative contact half-width only when q≥1/2.

Test the threshold rather than trusting the formula blindly

For q=2, recover a=1. For q=8, obtain a=1.5, so the contact interval grows to [−1.5,1.5]. Stronger downward loading increases contact, as expected. But for q=1/8 the formula gives a=−2, which is not a meaningful central contact half-width.

The correct weak-load solution has no contact. Solve u″=q with the endpoint values to get u=1+(q/2)(x²−4). Its minimum at x=0 is 1−2q. When q=1/8, that minimum is 3/4>0, so the obstacle remains inactive. At q=1/2, the minimum reaches zero and contact first occurs at a point.

This exercise tests whether you understand the assumed regime. A derived formula can be algebraically correct and still cease to describe the problem when its own geometric assumptions fail. Check inequalities and feasibility after solving, rather than treating them as background information.

Now transfer the same habit to melting. If the initial layer thickness is s₀>0, retain s₀² in the square-root law. If the solid carries a nonzero heat flux, restore q_right in the Stefan balance. If the wall supplies constant heat flux rather than constant temperature, the flux law changes. The boundary condition is part of the mechanism, so changing it can change the growth law.

Previous chapter · Contents · Next chapter

Begin with conservation laws, elementary heat conduction and the chain rule. Derive the one-dimensional Stefan condition with arrows before using compact jump notation. Calculate the quasi-steady front law, compare it with the similarity benchmark and then solve the obstacle example. Those calculations reveal the main ideas without requiring a large simulation.

A useful small project contains three independent tests. First, reproduce the moving-coordinate identity using T=x. Second, calculate the similarity front for Ste=0.1 and compare it with the approximate square-root law. Third, recover the obstacle contact interval for q=2 and verify the no-contact result for q=1/8. Each test isolates a different source of mistakes.

Does “free” mean unconstrained?

No. The boundary is unknown in advance, but its location is constrained by the equations, interface laws and data. It may be more tightly determined than a prescribed boundary because several conditions must hold together. “Free” describes its status as an unknown, not an absence of physical or mathematical restrictions.

Is a free boundary always the exterior edge?

No. The interface between solid and liquid can lie entirely inside a container. The boundary of a contact set can lie inside a fixed membrane domain. A model may therefore contain both fixed external boundaries and internal free boundaries, with different conditions on each.

Can a correct field still give a wrong boundary?

Yes. An interface may be inferred from a threshold, phase fraction or gradient balance that is more sensitive than the field itself. Check the inferred geometry directly. In phase change, also verify latent-energy accounting; in contact, verify the inequality and reaction conditions.

Is a sharp interface always preferable?

No. A sharp interface is appropriate only when it represents the relevant physics at the chosen scale. A material with a substantial mixed-phase region may need a different model. Numerical smearing and a physical mixed-phase region are also different things: one is an approximation error, the other may be a feature of the material.

What belongs in a defensible result?

State the unknown fields and geometry, the governing equations, normal orientation, interface conditions, initial and outer boundary data, parameter values and assumptions. Then show a benchmark, conservation or complementarity checks, and a refinement study appropriate to the desired output. Explain the limitations that would change the conclusion.

The linked Oxford notes provide a route into phase change, Ros-Oton’s overview into obstacle problems, and Peskin’s paper into fluid–structure coupling. Read them for their particular mathematical settings rather than combining their interface laws indiscriminately. The common discipline is to determine the geometry and field together, then test the condition that makes the boundary part of the solution.

Previous chapter · Contents · Gold Standard master