Analytical simulation of a cooling channel in an engine block - temperature evolution in time

I would start out by looking at some simplified limiting cases, just to get some results under my belt.

One case I would look at is the case where the coolant is flowing so fast that its temperature does not have a chance to change during the time it is in the block. In this case, the temperature that the coolant presents would always be To. This would represent the fastest that the cooling could take place. It would also simplify the math by removing the axial dependence, reducing the problem to transient radial conduction. I'm sure there is an analytic solution to this.

Next, I would look at a case where the coolant in the channel is very well mixed axially. So it would be like a continuous stirred tank, with coolant entering at $T_0$, and exiting at the well mixed temperature of T(t). And, here again, the temperature it would present to the block would be T(t), which changes with time, but not axial position. This problem too would have no axial dependence for the block, and the temperature would vary only with r and t.

There are other simplifying limiting cases I can think of, but I think I'll stop here for now.

As a simulation I'm using a $2D$ rectangle with a $1D$ cooling channel (of radius $R_{in}$) A uniform, constant heat source $q$ represents engine heat:

The principal equation here is:

$$alpha(u_{xx}+u_{yy})+q=0$$ with boundary conditions (BC): $$u_x(0,y)=u_x(L,y)=0$$ $$u_y(x,+H/2)=u_y(x,-H/2)=0$$ These BCs represent no heat loss from the surface of the engine block.

We're interested in the steady-state regime and by definition this means the coolant carries of all the engine heat, or: $$HLq=dot{m}c_p(T_1-T_0)Rightarrow T_1=T_0+frac{HLq}{dot{m}c_p}$$ For an infinitesimal cooling channel element $text{d}x$, with Newton's law of cooling: $$text{d}dot{q}=-htext{d}A[u(x,0)-T(x)]$$ $$frac{text{d} T}{text{d} x}=-frac{2pi R_{in} h}{dot{m}c_p}[u(x,0)-T(x)]$$ Now we make another assumption (that will get rid of the coupling): $$T(x)approx T_0+frac{HLq}{dot{m}c_p}x$$ $$frac{text{d} T}{text{d} x}=-frac{2pi h R_{in}}{dot{m}c_p}[u(x,0)-T(x)]=-frac{HLq}{dot{m}c_p}$$ $$[u(x,0)-T(x)]=-frac{HLq}{2pi h R_{in}}$$ $$h[u(x,0)-T(x)]=-ku_y(x,0)$$ $$u_y(x,0)=frac{HLq}{2pi hk R_{in}}$$

---

Briefly recapping: $$u_t=0 Rightarrow alpha(u_{xx}+u_{yy})+q=0tag{1}$$ $$u_x(0,y)=u_x(L,y)=0$$ $$u_y(x,+H/2)=u_y(x,-H/2)=0$$ $$u_y(x,0)=frac{HWq}{2pi hkR_{in}}$$

---

$(1)$ is of course non-homogeneous, so separation of variables won't work. But a comrade from math.SE came to the rescue with the following substitution: $$u(x,y)=w(x,y)-frac{q}{4alpha}(x^2+y^2)$$ $$alphaBig(w_{xx}-frac{q}{2alpha}+w_{yy}-frac{q}{2alpha}Big)+q=0$$ $$alpha w_{xx}-frac{q}{2}+alpha w_{yy}-frac{q}{2}+q=0$$ $$w_{xx}+w_{yy}=0$$ Neat! Now we transcribe the BCs: $$u_x(0,y)=w_x(0,y)=0$$ $$u_x(L,y)=w_x(L,y)=0$$ $$u_y(x,-H/2)=w_y(x,-H/2)=0$$ $$u_y(x,+H/2)=w_y(x,+H/2)=0$$ $$u_y(x,0)=w_y(x,0)=frac{HLq}{2pi hkR_{in}}$$ Ansatz: $$w(x,y)=X(x)Y(y)$$ Separation: $$frac{X''}{X}+frac{Y''}{Y}=0Rightarrow frac{X''}{X}=-frac{Y''}{Y}=-m^2$$ The $X$ ODE: $$X''(x)+m^2 X(x)=0$$ $$X(x)=Asin mx +Bcos mx$$ $$X'(0)=mAcos 0-m Bsin 0=0 Rightarrow A=0$$ $$X(x)=Bcos mx$$ $$X'(L)=-mBsin mL=0$$ $$mL=npiRightarrow m=frac{npi}{L}$$ for $n=1,2,3,...$ So: $$X_n(x)=B_ncosBig(frac{npi x}{L}Big)$$

---

The $Y$ ODE: $$Y''(y)-m^2 Y(y)=0$$ $$Y(y)=Ce^{-my}+De^{+my}$$ $$Y'(y)=-mCe^{-my}+mDe^{+my}$$ $$Y'(H/2)=-mCe^{-mH/2}+mDe^{+mH/2}=0$$ $$-Ce^{-mH/2}+De^{+mH/2}=0$$ $$Y'(0)=-mCe^{0}+mDe^{0}=frac{HLq}{2pi hk R_{in}}$$ $$-mC+mD=frac{HLq}{2pi hk R_{in}}$$ This looks messy but develops quite easily to: $$Y_n(y)=frac{HLq}{2pi hk R_{in}m(1-e^{mH})}Big(e^{m(H-y)}+e^{my}Big)$$

---

So that: $$w_n(x,y)=B_nfrac{HLq}{2pi hk R_{in}m(1-e^{mH})}Big(e^{m(H-y)}+e^{my}Big)cosBig(frac{npi x}{L}Big)$$

---

Using the usual Fourier expansion: $$B_n=frac{2}{L}frac{1}{H}int_0^L int_0^{H/2}T_i(x)X_n(x)Y_n(y)text{d}xtext{d}y$$ $$B_n=frac{2}{L}frac{1}{H}frac{HLq}{2pi hk R_{in}m(1-e^{mH})}int_0^L int_0^{H/2}T_1(x)Big(e^{m(H-y)}+e^{my}Big)cosBig(frac{npi x}{L}Big)text{d}xtext{d}y$$ $$pi hk R_{in}=mu$$ $$B_n=frac{q}{mu m(1-e^{mH})}int_0^L int_0^{H/2}T_i(x)Big(e^{m(H-y)}+e^{my}Big)cosBig(frac{npi x}{L}Big)text{d}xtext{d}y$$ Let: $T_i(x)=T_i(1-frac{1}{2L}x)$ Then: $$int_0^{H/2}Big(e^{m(H-y)}+e^{my}Big)text{d}y$$ $$=frac{e^{mH}-1}{m}=frac{L(e^{mH}-1)}{npi}$$ $$int_0^L[T_i(1-frac{1}{2L}x)]cosBig(frac{npi x}{L}Big)text{d}x$$ $$=frac{2LT_i}{pi^2 n^2}(1-(-1)^n)$$ $$B_n=frac{2qL^2T_i}{mu n^4 pi^4}(1-(-1)^n)$$ $$w_n(x,y)=B_nfrac{HL^2q}{2mu npi (1-e^{mH})}Big(e^{m(H-y)}+e^{my}Big)cosBig(frac{npi x}{L}Big)$$ $$w_n(x,y)=frac{q^2HL^5 T_i}{2mu^2n^5pi^5(1-e^{mH})}(1-(-1)^n)Big(e^{m(H-y)}+e^{my}Big)cosBig(frac{npi x}{L}Big)$$ $$f(n)=frac{(1-(-1)^n)}{n^5(1-e^{mH})}$$ This function very quickly decays to $approx 0$, so a 'first term approximation' is in order. $$u(x,y)=frac{q^2HL^5T_i}{2mu^2pi^5}displaystylesum_{n=1}^{infty}f(n)Big(e^{m(H-y)}+e^{my}Big)cosBig(frac{npi x}{L}Big)-frac{q}{4alpha}(x^2+y^2)$$ $$u(x,y)approx - 0.903 frac{q^2HL^5T_i}{2mu^2pi^5}Big(e^{m(H-y)}+e^{my}Big)cosBig(frac{pi x}{L}Big)-frac{q}{4alpha}(x^2+y^2)$$ where $m=frac{pi}{L}$. For $L=H=1$ we can write: $$u(x,y)=alpha Big(e^{pi(1-y)}+e^{pi y}Big)cos(pi x)+beta(x^2+y^2)$$ For $alpha=24.7$ and $beta=944$, then $u(1,0)=373$ and $u(0,frac{1}{2})=473$ and we obtain the following plots (Wolfram alpha):

A plot for $y=0$:

And one for $x=0$:

I have another approximation for you. Suppose that the mass flow rate of coolant is high enough that its temperature does not change much in passing through your model block. Then we can use the average of To and the exit temperature to calculate the heat flux from the engine block to the coolant. Since you are neglecting the thermal inertia of the holdup coolant within the duct, we can then write, $$dot{m}C_p(v_f -T_0)=2pi R_{in}hleft(u(R_{in},t)-frac{(T_0+v_f)}{2}right)L$$The solution to this equation for $v_f$ is: $$v_f=frac{lambda-1}{lambda+1}T_0+frac{2}{lambda+1}u(R_{in},t)$$where $lambda=frac{dot{M}C_p}{pi R_{in}hL}$ That gives an average coolant temperature of $$frac{(T_0+v_f)}{2}=frac{lambda}{lambda+1}T_0+frac{1}{lambda+1}u(R_{in},t)$$and an average temperature driving force of $$u(R_{in},t)-frac{(T_0+v_f)}{2}=frac{(u(R_{in},t)-T_0)}{1+1/lambda}$$From this, it follows that the boundary condition at $r=R_{in}$ becomes: $$-ku_r=frac{h}{1+1/lambda}(u(R_{in},t)-T_0)$$This effectively eliminates the coolant from further consideration. And, the solution to all this is a function only of r and t. This is equivalent to a case in which the coolant is at To throughout, and the heat transfer coefficient is lower by a factor of $1+1/lambda$. Note that, as the coolant mass flow rate becomes larger, $1/lambda$ approaches zero and the reduction factor approaches 1.