Solving second order non-linear partial differential in three independent variables

Temperature(T) is function of r, teta, z (cylindrical coordinates)

V1(r)*(1/r)* + V2(r)* = a* ( (1/r)* + )

where

a is constant

V1 and V2 are function of r

can we solve this differential equation in matlab

Answer 1

[William Rose]

@Mahesh Nallamothu, Yes you can.

[Editing because I made a mistake in analyzing your postiing.]

Matlab's PDE toolbox is very powerful. You will have to invest some effort to master its powers.

If you do not have access to the PDE toolbox, or you do not have time to learn it, you can solve your PDE directly. Write the discretized version of the PDEs, i.e. the difference equations. Initialize the array elements that correspond to the boundary conditions. Then update other elements, using the difference equations.

See the answer I posted here, to a similar question.

https://www.mathworks.com/matlabcentral/answers/1641800-explicit-heat-conduction-in-3d

I notice that your equations do not involve time. Therefore I assume you are trying to fnd a steady state solution that is consistent with V1(r) and V2(r). Can you give an example of what V1(r) and V2(r) might be?

It may be necessary to solve this by a relaxation method. In this approach, you create a 4-dimensional array, where the 4th dimension is for time. You write the difference equations for the PDE for heat diffusion. You initialize the system with some arbitrary distribution of temperature that is consistent with V1 and V2 but probably not consitent with the PDEs you gave. Then you enforce the conditions on V1 and V2 and choose a small time step size. (The time step, delta T, will have to be determined by trial and error. If delta T is not small enough, you will get oscillations. If it is too large, you will need to take many steps before the system reaches equilibrium.) You apply the difference equaitons at each new time step. Let the system evovle over time until it steops changing significantly.

Good luck!

[Torsten]

@William Rose

There is no time coordinate in the equation.

It seems like the 3d-Laplace equation in cylindrical coordinates.

Forget it - "theta" seems to be (maybe dimensionless) time.

[Mahesh Nallamothu]

@William Rose

v2(r) = a*( - r2*r2 - (b*log(r/r2)));

v1(r) = w*r1*((r2/r) - (r/r2))/(c - (1/c));

r1, r2, w, c, a, b are constants

can i get a code

[William Rose]

What is the range of values for r? V1(r), V2(r) become undefined at r=0.

[William Rose]

Do we know if r1>r2, or vice versa?

[William Rose]

@Mahesh Nallamothu,

Would you please check your orignal differential equation, especially the term related to θ? I am asking because the term , on the left hand side, surprises me. What is the range for θ? , or something less than that? What is the value of when θ goes from to 0?

[William Rose]

@Mahesh Nallamothu,

[edited: I replaced the pdf file because the initial verison had a mistake.]

[edited again to correct a mistake in equation 1, which propagated into other equations, and edited to reflect comments by @Mahesh Nallamothu which explained the units for some quantitites.]

The relaxation method is described nicely in Numerical Recipes, section 20.5 of the 3rd edition, which you can read online here.

http://numerical.recipes/book/book.html

I recommend using Jacobi's method, equation 20.5.5, even though the authors say it is not used because it is inefficient.

My attempt to convert the relaxation method to a set of difference equations is attached as a word doc. Check my work carefully, and if you agree, implement it.

[Mahesh Nallamothu]

@William Rose

[Editing this because i had given wrong inputs of r1 and r2 and ]

r1= 0.0127

r2 = 0.0205

r1<r<r2

In a double pipe counter current heat exchanger inner tube is rotated with angular velocity w, we are trying to derive temperature profile in annular part (r1<r<r2). We assumed that the system is in steady state, no heat source, no viscous dispersion and conduction is negligible in teta and z direction, on solving the energy equation in cylindrical coordinates with above assumptions we will get V1(r)*(1/r)* + V2(r)* = a* ( (1/r)* + ) where V1(r), V2(r) are velocity profiles in teta and z direction respectively, no velocity in r direction.

[Mahesh Nallamothu]

@William Rose

In the file relaxation method pdf, we can use that if it is unsteady state but our case it is steady state

V1 and V2 equations are attached please refer the attached pdf

Answer 2

[William Rose]

@Mahesh Nallamothu,

An analytical solution exists: Consider a solution which is independent of z and θ . In other words, T(r) is the same for all values of z and for all values of θ . Then . Then the differential equation simplifies to . Then we remember the derivative of natural log, from calculus class a million years ago (it seems). We observe that

is a solution to the diferential equation. k1 and k2 can be chosen to meet boundary conditions at the inner and outer radii. I assume the volume of interest is an annular cylinder, because several terms in the equations diverge at r=0.

[Mahesh Nallamothu]

@William Rose

we have already done by assuming = 0 and now we are trying out by increasing complexity. so, cannot be assumed to be zero

So to solve this we need a best sutible matlab code or numerical method

Answer 3

[Torsten]

Interpret teta as time.

Discretize the dT/dz, 1/r*d/dr(r*dT/dr) expressions in space.

You'll arrive at a system of ordinary differential equations of the form

dT/dteta = A(r,z)*T + b(r,z)

Advance the solution in time e.g. using the Crank-Nicolson method.

If it's difficult to program this on your own, maybe a google search for

"instationary Laplace equation in cylindrical coordinates & matlab"

might help.

[Mahesh Nallamothu]

@Torsten

we cannot use Crank-Nicolson method since our differential equation is independent of time.

Can please discretize and write a code (didn't find any code in google )

[Torsten]

But you can interprete teta as time, and then the usual Crank-Nicolson comes into play.

Can please discretize and write a code (didn't find any code in google )

This is overcharged, but if you make an attempt, the forum usually is willing to help.

[William Rose]

@Mahesh Nallamothu,

When @Torsten makes a comment, I pay attention, becase he knows a great deal about math and about Matlab. I recommend trying anything he recommends.

I have doubts about the original equaitons which you posted. I am afraid there is a probem with dimensions. When I made the dimensions comment in my PDF file, I did not have all the info which you have since posted. Now I understand that v1 and v2 are linear velocities (units=L/T) and w is angular velocity (units=1/T). It follows from the equaiton for velocity v1(r) that the constant "a" has units=. , and that b has units=. I assume the variiable T is temperature and I will call its units K. The right hand side of original differential equation is

which has units=.

The first term on the left hand side of the original differential equation is

which has units of . This is incompatible with the units of the right hand side of the differential equation, analyzed above. The other terms on the left hand side of the orignal differential equation is

which has units of . This matches the units for the other term on the left hand side, but is incompatible with the units on the right hand side. Since the units don't match, there must be a mistake in the original differential equation.

[Torsten]

@William Rose

I think a = lambda/(rho*cp) with unit m^2/s.

Then the right-hand side has unit K/s.

This corresponds to v2*dT/dz with unit m/s * K/m = K/s.

To make the last unit be correct, theta has to be dimensionless which is correct if it's the angular direction.

So I must correct my previous comment:

It's better to interprete z as time in the Crank-Nicolson approach.

[William Rose]

@Mahesh Nallamothu, @Torsten,

My conclusion that a has units of is based on the equations in the PDF posted by @Mahesh Nallamothu:

If v1 is linear velocity, then a must have units . But if one uses those units for a in the original differential equation, the right and left hand side units do not match. If a is what @Torsten says it is, then v1 , given by the equation above, does not have units of linear velocity. I suspect that the solution to this problem is that there are two different constants with different units. @Mahesh Nallamothu, what do you think?

[Mahesh Nallamothu]

@William Rose @Torsten

I had attached a document of differential equation. All the constants had been defined in that document.

[William Rose]

@Mahesh Nallamothu, The handwritten note you posted shows that "a" in the definition of v1(r) is not the same as "alpha" on the right hand side of the differential equation. In your original posting, you used "a" for both, which created a problem with units. That issue is now resolved.

[William Rose]

@Mahesh Nallamothu,

I have a question about your equation for v2(r). This is the velocity in the z direction, parallel to the pipe axis.

The domain is R1<=r<=R2. Notice that v2=0 when r=R2, as expected, if there is no slip at the wall.

However, v2 is not = 0 when r=R1. If we invert the fraction that is circled in red, then v2=0 when r=R1. Do you think the fraction circled in red shoud be inverted?

[William Rose]

@Mahesh Nallamothu, I think there is another problem with the equation

Consider the case when R1 goes to zero and R2 becomes R. This equaiton should simplify to the standard equaiton for Poiseuille flow in a pipe with radius R. However, it looks to me like it simplifies to

which is wrong by a factor of -1.

[Mahesh Nallamothu]

@William Rose

In our case R1 is finite not equal to zero and also r=R1 v2 is not equal to zero