10
$\begingroup$

My question

Consider the following circuit. A lossless long transmission lines ($R' = G' = 0 \text{ } \Omega$), of length $\ell$ and inductance per unit length $L'$ and capacitance per unit length $C'$, with one of its ports short-circuited by ideal wires. In the other port of the transmission line, we connect an ideal independent sinusoidal voltage source of voltage $v_\text{s}(t) = V_\text{m} \cos{(\omega t)}$ in series with an ideal normally-open switch. The switch is closed at $t = 0 \text{ s}$.

Before that instant, the circuit has been for a long time in the configuration shown, meaning it is operating in steady-state, so the inductors and capacitors act as short-circuits and open-circuits, respectively; all voltages and currents are zero.

The circuit is shown in the following figure.

Figure 1

Figure 1. Image source: own.

Given the above assumptions, I'd like to determine the current $i(z,t)$ flowing through the transmission line as a function of space $z$ and time $t$, the voltage $v(z,t)$ across the transmission line as a function of space and time. I want to solve them analytically (if possible) or numerically. For simplicity, let's assume that $L' = 1 \text{ H}$, $C' = 1 \text{ F}$, $\ell = 5 \text{ m}$, $V_\text{m} = 1 \text{ V}$, $\omega = 1 \text{ rad/s}$, and the maximum time we're interested on is $t_\text{max} = 20 \text{ s}$.


My attempt

As we know, the uncoupled telegrapher's equations for a lossless long transmission line in the time domain (which are pretty much the one-dimensional wave equation, a hyperbolic PDE) are:

$\left\{ \begin{align} \dfrac{\partial^2 v}{\partial z^2} &= L' C' \dfrac{\partial^2 v}{\partial t^2} \\ \dfrac{\partial^2 i}{\partial z^2} &= L' C' \dfrac{\partial^2 i}{\partial t^2} \end{align} \right. \quad , \, 0 \text{ m} \le z \le \ell, \, t \ge 0 \text{ s} \tag 1$

As for the boundary conditions:

  • At the sending-end (where the voltage source is) of the transmission:

$v(0, t) = V_\text{m} \cos{(\omega t)} \tag 2$

  • At the receiving-end (short-circuited port) of the transmission:

$v(\ell, t) = 0 \text{ V} \tag 3$

  • Before the switch is closed, the voltage and current at all locations is zero:

$v(z, 0 \text{ s}) = 0 \text{ V} \tag 4$

$i(z, 0 \text{ s}) = 0 \text{ V} \tag 5$

Is the above correct?

I tried to solve them in Mathematica with NDSolve using the following code:

LL = 1;
CC = 1;
L = 5;
tmax = 20;
Clear[v, i, z, t];
eqs = {
   D[v[z, t], {z, 2}] == LL*CC*D[v[z, t], {t, 2}],
   D[i[z, t], {z, 2}] == LL*CC*D[i[z, t], {t, 2}]
   };
bcs = {
   v[0, t] == Cos[t]*UnitStep[t],
   v[L, t] == 0,
   v[z, 0] == 0,
   i[z, 0] == 0
   };
sol = NDSolve[Join[eqs, bcs], {v, i}, {z, 0, L}, {t, 0, tmax}]

But Mathematica can't solve them:

Figure 2

Figure 2. Image source: own.

Is it because there are too many or too few boundary conditions, or because there are too few equations, or because Mathematica can't simply solve them (which doesn't seem likely since it can solve much more complex equations), or because of something else?

$\endgroup$
1
  • 1
    $\begingroup$ One thing that is missing are the derivatives of the initial conditions. $\endgroup$
    – user21
    Dec 18 '21 at 2:31
5
$\begingroup$

The model of transmission line includes two equations: $$\frac{\partial V}{\partial x}=-L_0\frac{\partial I}{\partial t},\frac{\partial I}{\partial x}=-C_0\frac{\partial V}{\partial t}$$ This system can be written also as one wave equation as follows $$\frac{\partial^2 V}{\partial x^2}=C_0L_0\frac{\partial^2 V}{\partial t^2}$$

In a case of bcs by Alejandro Nava we have

LL = 1;
CC = 1;
L = 5;
tmax = 20;
Clear[v, i, z, t];
eqs = {D[v[z, t], {z, 2}] == LL*CC*D[v[z, t], {t, 2}]}; eq1 = 
 D[i[z, t], {z, 2}] == LL*CC*D[i[z, t], {t, 2}];
bcs = {DirichletCondition[v[z, t] == Cos[t], z == 0], 
  DirichletCondition[v[z, t] == 0, z == L]}; ics = {v[z, 0] == 0, 
  Derivative[0, 1][v][z, 0] == 0};
sol = NDSolve[Join[eqs, bcs, ics], v, {z, 0, L}, {t, 0, tmax}]

We can check mesh using in this example

v["ElementMesh"] /. sol

Out[]= {NDSolve`FEM`ElementMesh[{{0., 
    5.}}, {NDSolve`FEM`LineElement["<" 20 ">"]}]}

Therefore, mesh consists of 20 line elements only. Visualization

Plot3D[Evaluate[v[z, t] /. sol], {z, 0, L}, {t, 0, tmax}, Mesh -> None, 
  ColorFunction -> Hue, PlotPoints -> 50, 
  AxesLabel -> Automatic] // AbsoluteTiming

Figure 1

We also can solve system of equations, but since it consists of first order equations, we need to remove one bc on the right side at $x=L$, then we have

Needs["NDSolve`FEM`"]

LL = 1;
CC = 1;
L = 5;
tmax = 20;
Clear[v, i, z, t];
eqs = {D[v[z, t], {z, 1}] == -LL*D[i[z, t], {t, 1}], 
   D[i[z, t], {z, 1}] == -CC*D[v[z, t], {t, 1}]};
bcs = {DirichletCondition[v[z, t] == Cos[t], z == 0], 
  DirichletCondition[i[z, t] == 0, z == L]}; ics = {v[z, 0] == 0, 
  i[z, 0] == 0};
sol = NDSolve[Join[eqs, bcs, ics], {v, i}, {z, 0, L}, {t, 0, tmax}]

{Plot3D[Evaluate[v[z, t] /. sol[[1]]], {z, 0, L}, {t, 0, tmax}, 
  Mesh -> None, ColorFunction -> Hue, PlotPoints -> 50, 
  AxesLabel -> Automatic], 
 Plot3D[Evaluate[i[z, t] /. sol[[1]]], {z, 0, L}, {t, 0, tmax}, 
  Mesh -> None, ColorFunction -> Hue, PlotPoints -> 50, 
  AxesLabel -> Automatic]}

Figure 2

Also we have a message

NDSolve::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

It means that we need to add some resistance to our model taken into account physical effects like electromagnetic wave generation and heating.

$\endgroup$
6
  • 2
    $\begingroup$ Use Evaluate for the plot arguments to seed things up. $\endgroup$
    – user21
    Dec 19 '21 at 6:40
  • $\begingroup$ @user21 Thank you very much! With Evaluate we speed up 100 times Plot3D in v.13. But why Evaluate is not Automatic for Plot3D and what Mathematica doing during 112 s before plotting without Evaluate? :) $\endgroup$ Dec 19 '21 at 8:52
  • 3
    $\begingroup$ The plot functions are HoldAll - see Attributes. Your rule is evaluated for every plot point. This nothing new. $\endgroup$
    – user21
    Dec 19 '21 at 10:19
  • $\begingroup$ @user21 Ah, yes, you are right, my answer has been corrected. $\endgroup$ Dec 19 '21 at 15:45
  • 1
    $\begingroup$ @AlejandroNava See update to my answer with solution of the system of equations. $\endgroup$ Dec 20 '21 at 4:36
7
$\begingroup$

Give consistent initial and boundary conditions and derivatives of the initial condition:

LL = 1;
CC = 1;
L = 1;
tmax = 1;
Clear[v, i, z, t];
eqs = {D[v[z, t], {z, 2}] == LL*CC*D[v[z, t], {t, 2}], 
  D[i[z, t], {z, 2}] == LL*CC*D[i[z, t], {t, 2}]}
bcs = {v[0, t] == Cos[t](**UnitStep[t]*), v[L, t] == 0,
  v[z, 0] == 1 - z,
  Derivative[0, 1][v][z, 0] == 0,
  i[z, 0] == 0,
  Derivative[0, 1][i][z, 0] == 0}
sol = NDSolve[Join[eqs, bcs], {v, i}, {z, 0, L}, {t, 0, tmax}]

Plot3D[Evaluate[v[z, t] /. sol], {z, 0, L}, {t, 0, tmax}, 
 AxesLabel -> {"z", "time", "v(z,t}"}]

enter image description here

$\endgroup$
2
  • $\begingroup$ Thanks! It seems to work. Why did you use v[z, 0] == 1 - z instead of v[z, 0] == 0? Both seem to work and give warning messages, with different solutions of course. Comparison (using L = 5 and tmax = 20) Using v[z, 0] == 1 - z gives a more distorted waveform. $\endgroup$
    – alejnavab
    Dec 19 '21 at 22:26
  • $\begingroup$ I do not get a message with 1-z. This is to make the boundary and initial condition consistent. You have to get rid of these messages, the result will not be good if you do not do that. $\endgroup$
    – user21
    Dec 20 '21 at 7:52
3
$\begingroup$

Your equations seem to misbehave at large time and/or Large L.

Reducing time duration and L one can get values that makes it work.

It uses method of lines to march in time. At each time step it needs to solve the system. It seems at one point the system becomes numerically not solvable by LinearSolve as the matrix becomes singular. Something to do with nature of your equations for large time. Are you sure the Physical model you used is correct?

Strange thing, is that on V 13 it seems to give the pivot error when trying same values again. Here is one case where the error did not show up.

LL = 1;
CC = 1;
L = 1;
tmax = 1;
Clear[v, i, z, t];
eqs = {D[v[z, t], {z, 2}] == LL*CC*D[v[z, t], {t, 2}], 
  D[i[z, t], {z, 2}] == LL*CC*D[i[z, t], {t, 2}]}
bcs = {v[0, t] == Cos[t]*UnitStep[t], v[L, t] == 0, v[z, 0] == 0, i[z, 0] == 0}
sol = NDSolve[Join[eqs, bcs], {v, i}, {z, 0, L}, {t, 0, tmax}]


Plot3D[Evaluate[v[z, t] /. sol], {z, 0, L}, {t, 0, tmax}, 
 AxesLabel -> {"z", "time", "v(z,t}"}]

Mathematica graphics

Trying the above again, now it gives the pivot error!. So if you get the error, try again from clean kernel. Eventually it will work with no error. I have no idea why.

It seems NDSolve remembers something, or trying it one more time causes different result to show up. This should not happen as there is nothing random here. Do others see the same issue?

edit

I tried this in V 11.3 and it gives different looking solution

LL = 1;
CC = 1;
L = 5;
tmax = 10; (*20;*)
Clear[v, i, z, t];
eqs = {D[v[z, t], {z, 2}] == LL*CC*D[v[z, t], {t, 2}], 
  D[i[z, t], {z, 2}] == LL*CC*D[i[z, t], {t, 2}]}
bcs = {v[0, t] == Cos[t]*UnitStep[t], v[L, t] == 0, v[z, 0] == 0, i[z, 0] == 0}
sol = NDSolve[Join[eqs, bcs], {v, i}, {z, 0, L}, {t, 0, tmax}]

Plot3D[Evaluate[v[z, t] /. sol], {z, 0, L}, {t, 0, tmax}, 
 AxesLabel -> {"z", "time", "v(z,t}"}]

Mathematica graphics

ps. Very strange. After restarting V 13, now I get that pivot error again on the same code. I have no idea why it worked one time only and now it no longer works on V 13. Will try to find out.

$\endgroup$
1
  • $\begingroup$ "Are you sure the Physical model you used is correct?" I'm not sure but it looks correct. I've seen textboks model the voltage source/generator with a series resistance (to more correctly represent a physical/real generator). That changes the equations. Maybe in that there're no errors. I'll do more research and reply back. Thanks for the answer BTW! $\endgroup$
    – alejnavab
    Dec 18 '21 at 2:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.