What is the shape of the perfect coffee cup for heat retention assuming coffee is being drunk at a constant rate?

Question

Find the optimal shape of a coffee cup for heat retention. Assuming

1. A constant coffee flow rate out of the cup.
2. All surfaces radiate heat equally, i.e. liquid surface, bottom of cup and sides of cup.
3. The coffee is drunk quickly enough that the temperature differential between the coffee and the environment can be ignored/assumed constant.

So we just need to minimise the average surface area as the liquid drains

I have worked out the following 2 alternative equations for the average surface area over the lifetime of the liquid in the cup (see below for derivations):

$$ S_{ave} =\pi r_0^2+ \frac{\pi^2}{V}\int_{0}^{h}{{r(s)}^4ds}+\frac{2\pi^2}{V}\int_{0}^{h}{\int_{0}^{s}{r\sqrt{1+\left(\frac{dr}{dz}\right)^2}\ dz\ }{r(s)}^2ds\ } \tag{1}$$

$$S_{ave}=\pi r_0^2+\frac{\pi^2}{V}\int_{0}^{h}r\left(s\right)^4ds+\frac{2\pi^2}{V}\int_{0}^{h}r(s){\underbrace{\int_{s}^{h}{{r\left(z\right)}^2dz\ }}_{\text{Volume Drunk}}}\sqrt{1+\left(\frac{dr}{ds}\right)^2}ds \tag{2}$$

If the volume of the cup is constant

$$ V=\pi\int_{0}^{h}{{r(z)}^2dz\ }$$

Can the function, $r(z)$, be found that minimises the average surface area $S_{ave}$?

If r is expressed as a parametric equation in the form $r=f(t), z=g(t)$ and $f,g$ are polynomials then a genetic search found the best function of parametric polynomials to be: $r\left(z\right)=\sqrt{\frac{3}{2}}z^\frac{1}{2}-\frac{\sqrt6}{9}z^\frac{3}{2}, f\left(t\right)=\sqrt{\frac{3}{2}}t-\sqrt{\frac{3}{2}}t^3, g\left(t\right)=\frac{9}{2}t^2$

This parametric shape has a maximum radius of 1, height of 4.5, starting volume of $\frac{3^4}{2^5}\pi$ and is shown here:

I can't prove that there is (or is not) a better $r(z)$ but...

the average surface area of this surface turns out to be $12.723452r^2$ or $4.05\pi r_{max}^2$. I suspect that the optimal surface will have the same surface area as a sphere, i.e. $4\pi r_{max}^2$ $(12.5664)$

Conjecture: The optimally shaped coffee cup has the same average surface area as a sphere of the same maximum radius. Shown to be false by this answer

Derivation of Surface Area Formula:

Surface area when surface of liquid is at level s is the sum of the areas of the top disc, bottom disc and the sides.

$S(s)=\pi r_0^2+\pi r_s^2+2\pi\int_{0}^{s}{r(z)dldz}$

$S(s)=\pi r_0^2+\pi r_s^2+2\pi\int_{0}^{s}{r(z)\sqrt{1+\left(\frac{dr}{dz}\right)^2}\ dz\ }$

The average surface area will be the sum of all the As’s times the time spent at each surface area.

$S_{ave}=\frac{1}{T}\int_{t_0}^{t_h}{S(s)dt\ }$

In order to have the drain rate constant we need to set the flow rate Q to be constant i.e. the rate of change volume is constant and $Q=dV/dt =V/T$

Time spent at a particular liquid level $dt\ =\frac{T}{V}dV$ and $ dV={\pi r}^2ds$

$dt=\frac{T\pi r^2}{V}ds$

$S_{ave}=\int_{s=0}^{s=h}{S(s)\frac{T\pi{r(s)}^2}{V}ds\ }$ $S_{ave}=\frac{\pi}{V}\int_{s=0}^{s=h}{(r_0^2+r(s)^2+2\int_{z=0}^{z=s}{r(z)\sqrt{1+(\frac{dr(z)}{dz})^2}\ dz\ })\pi{r(s)}^2ds\ }$

$S_{ave}=\frac{\pi}{V}r_0^2\int_{0}^{h}{\pi{r(s)}^2ds}+\frac{\pi}{V}\int_{s=0}^{s=h}{\left(r\left(s\right)^2+2\int_{z=0}^{z=s}{r(z)\sqrt{1+(\frac{dr(z)}{dz})^2}\ dz\ }\right)\pi{r(s)}^2ds\ }$

$S_{ave}=\pi\ r_0^2+\frac{\pi^2}{V}\int_{s=0}^{s=h}{\left(r\left(s\right)^2+2\int_{z=0}^{z=s}{r(z)\sqrt{1+\left(\frac{dr(z)}{dz}\right)^2}\ dz\ }\right){r(s)}^2ds\ }$

Alternative Formula Derivation:

Surface area of highlighted ribbon in the diagram is:

$S_{ribbon}=2\pi rdl$

And the contribution towards the average surface area lasts for the ratio of volume of the liquid above the current level to the total volume.

$$S_{sides}=2\pi\frac{\pi}{V}{\underbrace{\int_{s}^{h}{{r\left(z\right)}^2dz\ }}_{\text{Volume Drunk}}}rdl=\frac{2\pi^2}{V}{\underbrace{\int_{s}^{h}{{r\left(z\right)}^2dz\ }}_{\text{Volume Drunk}}}r\left(s\right)\sqrt{1+\left(\frac{dr}{ds}\right)^2}ds$$

Integrate the contribution of all such sections.

$$S_{sides}=\frac{2\pi^2}{V}\int_{0}^{h}{r\left(s\right){\underbrace{\int_{s}^{h}{{r\left(z\right)}^2dz\ }}_{\text{Volume Drunk}}}\sqrt{1+\left(\frac{dr}{ds}\right)^2}ds}$$

Contribution of top surfaces to average surface area is area of top by proportion of volume that area x dz is:

$$S_{tops}=\frac{1}{V}\int_{0}^{h}{\pi{r(s)}^2{\pi r(s)}^2}ds$$

Contribution of bottom surface is constant $\pi r_0^2$ so adding together all three gives:

$$S_{ave}=\pi r_0^2+\frac{\pi^2}{V}\int_{0}^{h}r\left(s\right)^4ds+\frac{2\pi^2}{V}\int_{0}^{h}r(s){\underbrace{\int_{s}^{h}{{r\left(z\right)}^2dz\ }}_{\text{Volume Drunk}}}\sqrt{1+\left(\frac{dr}{ds}\right)^2}ds$$

---

Answer 1

This isn't really a solution - just rewriting it in terms of a fourth-order differential equation. I've also cross-posted this answer on MathOverflow.

---

Let $$V(t):=\pi\int_0^tr(s)^2ds$$

We can then rewrite $r(s)$ as $$r(s)=\sqrt{\frac{V'(s)}{\pi}}$$

If we rewrite $S_{ave}$ using just $V$, we get $$S_{ave}=\pi\sqrt{\frac{V'\left(0\right)}{\pi}}+\frac{\pi^{2}}{V\left(h\right)}\int_{0}^{h}\left(\left(\frac{V'\left(s\right)}{\pi}\right)^{2}+2\sqrt{\frac{V'\left(s\right)}{\pi}}\int_{s}^{h}\frac{V'\left(z\right)}{\pi}dz\sqrt{1+\frac{1}{4\pi}\cdot\frac{V''\left(s\right)^{2}}{V'\left(s\right)}}\right)ds$$

The inner integral simplifies to $\frac{1}{\pi}(V(h)-V(s))$, so this simplifies to (getting rid of some of the $\pi$ terms as well) $$\sqrt{\pi V'(0)}+\frac{1}{V(h)}\int_{0}^{h}\left(V'(s)^2+2\cdot(V(h)-V(s))\sqrt{\pi V'(s)+\frac{V''(s)^2}{4}}\right)ds$$

We want to minimize this value assuming a fixed $V(h)$. Assume we have a fixed $V'(0)$ and $h$ as well. We then want to minimize $\frac{\sqrt{\pi V'(0)}}{h}+\frac{1}{V(h)}\mathcal{L}(s,V, V', V'')$, where $\mathcal{L}(s,V, V', V'')$ is given by $$(V')^2+2\cdot(V_h-V)\sqrt{\pi V'+\frac{(V'')^2}{4}}$$

and $V_h=V(h)$. We can then use the Euler-Lagrange equation to get that the stationary points of the average surface area (with respect to $V(s)$) would be given by $$\frac{\partial\mathcal{L}}{\partial V}-\frac{d}{ds}\left(\frac{\partial \mathcal{L}}{\partial V'}\right)+\frac{d^2}{ds^2}\left(\frac{\partial \mathcal{L}}{\partial V''}\right)=0$$

This ends up being a fourth-order differential equation with a long form (I started writing it out before realizing that the last term would make it be very long).

---

Edit: Using the Beltrami identity, which TheSimpliFire mentioned in a comment, and this answer, we can write $$\mathcal{L}-V'\frac{\partial\mathcal{L}}{\partial V'}+V'\frac{d}{ds}\frac{\partial \mathcal{L}}{\partial V''}-V''\frac{\partial \mathcal{L}}{\partial V''}=C$$

Plugging in $\mathcal{L}$ and simplifying, we get $$-V'(s)^{2}\left(1+\frac{V''(s)\left(4\pi V'(s)+V''(s)^{2}\right)+4\pi\left(V(s)-V_{h}\right)\left(2\pi+V^{(3)}(s)\right)}{\left(4\pi V'(s)+V''(s)^{2}\right)^{\frac{3}{2}}}\right) = C$$

Answer 2

Solution:

For a coffee cup of volume $\frac{81\pi}{32}$ the shape of maximum heat retention is the surface of revolution of $$r(z)=\sqrt\frac{3}{2}z^\frac{1}{2}-\frac{\sqrt{6}}{9}z^{\frac{3}{2}}$$.

The height of this cup is $4.5$ and the max radius is $1$. For other starting volumes, V, the radius and height scale by $(\frac{32V}{81\pi})^\frac{1}{3}.$

Proof:

In the question I show that the average surface area for a coffee cup is:

$$S_{ave}=\pi r_0^2+\frac{\pi^2}{V}\int_{0}^{h}r\left(s\right)^4ds+\frac{2\pi^2}{V}\int_{0}^{h}r(s){\underbrace{\int_{s}^{h}{{r\left(z\right)}^2dz\ }}_{\text{Volume Drunk}}}\sqrt{1+\left(\frac{dr}{ds}\right)^2}ds \tag{2}$$ and ask what function of $r$ minimises this. Recall that I give $(1)$ as the best equation I have found using a genetic search algorithm.

TheSimpliFire showed in this answer and Varun Vejalla show in this answer (with a scaling of $\pi$ difference) that the optimal shape must be a solution to the third-order ODE: $$\frac{4(P(h)-P)(2+P''')-P''(4P'+P''^2)}{(4P'+P''^2)^{3/2}}=1+\frac C{P'^2}\tag{3}$$

where $P'=r^2(z)$

Let's plug $(1)$ into $(3)$ and see if it is a solution.

Using $P=\frac{3\ z^2}{4}-\frac{2\ z^3}{9}+\frac{z^4}{54}+K$ (Where $K$ is the constant of integration).

The numerator of the LHS evaluates to: $$\frac{81}{4}-8\left(\frac{3\ z^2}{4}-\frac{2\ z^3}{9}+\frac{z^4}{54}\right)+4\left(\frac{81}{32}-\left(\frac{3\ z^2}{4}-\frac{2\ z^3}{9}+\frac{z^4}{54}\right)\right)\left(\frac{4\ z}{9}-\frac{4}{3}\right)-4\left(\frac{3\ z}{2}\ -\frac{2\ z^2}{3}\ +\frac{2\ z^3}{27}\right)\left(\frac{3}{2}-\frac{4\ z}{3}+\frac{2\ z^2}{9}\right)-\left(\frac{3}{2}-\frac{4\ z}{3}+\frac{2\ z^2}{9}\right)^3$$

$$=\frac{27}{8}+\frac{9z}{2}+\frac{z^2}{2}-\frac{28z^3}{27}-\frac{2z^4}{27}+\frac{8z^5}{81}-\frac{8z^6}{729}$$

$$=-\frac{(-27-12z+4z^2)^3}{5832}$$

The denominator is:

$$\left(4\left(\frac{3\ z}{2}\ -\frac{2\ z^2}{3}\ +\frac{2\ z^3}{27}\right)+\left(\frac{3}{2}-\frac{4\ z}{3}+\frac{2\ z^2}{9}\right)^2\right)\frac{3}{2}$$

$$=\frac{\left(\left(-27\ -\ 12\ z\ +\ 4\ z^2\right)^6\right)^\frac{1}{2}}{5832}$$

Which is always equal to either the numerator or its negative

$$\frac{-(-27-12z+4z^2)^3}{\sqrt{(27 + 12 z - 4 z^2)^6}} = \pm1\tag{4}$$ But $P'=r^2 \implies r=\pm \sqrt {P'}$

So we can rewrite the ODE as $$\frac{4(P(h)-P)(2+P''')-P''(4P'+P''^2)}{\pm(4P'+P''^2)^{3/2}}=1+\frac C{P'^2}\tag{3a}$$

squaring both sides of (3a) and (4), we see that $(1)$ satisfies the ODE with $C=0$

Q.E.D.

Well, almost. The function is at a stationary point. Any small change in the coefficients increases the surface area so it’s not a maximum, so it must be a minimum. Now I just need to see if it is also a global minimum.

Answer 3

We are interested in minimising $$I(h)=\int_0^hr(s)^4+2r(s)\sqrt{1+r'(s)^2}\int_s^hr(z)^2\,dz\,ds$$ such that $V=\pi\int_0^hr(s)^2\,ds$. The substitution $P'=r^2$ with $P(h)=V/\pi$ means that $r'=P''/2r$ so $$I(h)=\int_0^hP'^2+\sqrt{4P'+P''^2}\left(P(h)-P\right)\,ds.$$ Let $L(P,P',P'')$ denote the integrand. As this is independent of $s$, we can derive the second-order Beltrami identity as follows.

Recall the total derivative $L'=L_s+L_PP'+L_{P'}P''+L_{P''}P'''$ which rearranges to $$L_PP'=L'-L_{P'}P''-L_{P''}P'''$$ as $L_s=0$. The second-order Euler-Lagrange equation is $$L_P-L_{P'}'+L_{P''}''=0\implies L_PP'-L_{P'}'P'+L_{P''}''P'=0.$$ Substituting the first equation gives $L'-L_{P'}P''-L_{P''}P'''-L_{P'}'P'+L_{P''}''P'=0$ which can be rewritten as $(L+(L_{P''}'-L_{P'})P'-L_{P''}P'')'=0$. Thus $$L+(L_{P''}'-L_{P'})P'-L_{P''}P''=C$$ where $C$ is a constant, $$L_{P'}=2P'+\frac{2(P(h)-P)}{\sqrt{4P'+P''^2}},\quad L_{P''}=\frac{P''(P(h)-P)}{\sqrt{4P'+P''^2}}$$ and $$L_{P''}'=\frac{(P(h)-P)(4P'P'''-2P''^2)}{(4P'+P''^2)^{3/2}}-\frac{P'P''}{\sqrt{4P'+P''^2}}.$$ Plugging all of these expressions in and simplifying gives $$\frac{4(P(h)-P)(2+P''')-P''(4P'+P''^2)}{(4P'+P''^2)^{3/2}}=1+\frac C{P'^2}$$ which is now a third-order ODE.

Equivalently, substituting $Q=P(h)-P$ yields $$\frac{4Q}{Q''}R'(Q)+Q''R(Q)=1+\frac C{Q'^2}$$ where $Q''^2-4Q'=R^{-2}$. Differentiating both sides gives $$Q'''=2+\frac{(R^{-2})'}{2Q''}$$ but this is of Abel type with no known closed form solution.

A minor result is that if $r(s)$ is a signomial then $\deg r\ge1$, as the ODE would imply the equality $\max\{2\deg P-3,3\deg P-6\}=\max\{\deg P-1,2\deg P-4\}$ forcing $\deg P\ge3$.

Answer 4

A variational calculus approach.

$\int 2 \pi r ds + \pi r^2 $ is to be minimised for constant $ \int \pi r^2 dz$. In Cartesian coordinates

$$2 \pi\int y \sqrt {1+y^{'2}}dx + \lambda \pi y^2 -\pi\int y ^2 dx \tag 1 $$

The second term

$$ 2 \lambda \pi \int y dy = 2 \lambda \pi \int y y' dx \tag 2 $$

The third term is modified with an arbitrary constant $2H $ of reciprocal linear dimension turning volume term also to second degree.

Dividing by $2 \pi$ the Lagrangian takes the form

$$ y\left( \sqrt {1+y^{'2}} + \lambda y' - H y \right) =F \tag 3$$

From Euler Lagrange Equation Beltrami pde integration without explicit $x$ term

$$ F -y'~\frac{\partial F}{\partial y'} = c $$

$$ y ~ ( \sqrt {1+y^{'2}} + \lambda ~y' - Hy )- y'\cdot y \left( \frac{y'}{\sqrt {1+y^{'2}}} +0 \right)=c \tag 4 $$

$$ \frac{ y}{ \sqrt {1+y^{'2} }}+ y^2 H = c \tag 5 $$

The quantity with radical sign is $ \sec \phi $; simplifying,

$$ \cos \phi = H y +\frac {c}{y} ~~ \tag 6$$

Differentiate w.r.t meridional arc $s$

$$ -\sin \phi \frac{d \phi}{ds} = ( H-\frac{c}{y^2})\sin \phi\tag 7$$

Simplifying and introducing principal curvatures with proper sign convention

$$ \kappa_2=\frac{\cos \phi}{y}; ~ \kappa1=\frac{d \phi}{ds}; \tag 8$$

$$ \frac{\kappa_1+\kappa_2 }{2}= H \tag 9 $$

This results in CMC constant mean curvature DeLaunay surface as a solution surface of revolution.

A sphere is a particularly simple case, $c=0$.

As is known, a CMC surface encloses maximum volume for given surface area.

Among three CMC surface of revolution varieties is the following unduloid of wavy /corrugated shape. Two other toroidal types are not likely to be really usable as a coffee cup... although a workaround might still be possible.

EDIT1:

Did not take the calculation to its full logical conclusion including drinking rate/time etc.( reg. physics of formulation). For large number of unduloid waves an asymptotic value of $ Vol^{1/3}/Area^{1/2} \approx 0.44 $ seems attained that is less than OP's value 0.56;