| Author | Message | ||
     DukeVN (dukevn) New member Username: dukevn Post Number: 1 Registered: 04-2008 |
When solving div(grad(u)) = 0 and plot the results, contour(u) draws the equipotential lines, vector(grad(u)) draws electric vectors laying on electric line of force. Anybody knows how to draw the electric line of force? In matlab, there is a function can do this job called streamslice, but I do not see any similar function in FlexPDE. Thanks. | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 1103 Registered: 06-2003 |
There is currently no built-in field line plot in FlexPDE, but you can often generate them as contour lines of the adjoint equation. See "Samples | Steady_State | Electricity | Fieldmap.pde". The tricky part about this is the boundary conditions for the adjoint equation. | ||
| Archie Campbell (amc1) Junior Member Username: amc1 Post Number: 3 Registered: 03-2008 |
I have often needed to do this but it usually involves solving the problem again with a different potential. In the electrical example quoted the vector D is the one with zero divergence which means it is possible to draw its lines of force. To do this we use a new potential Q defined by curl(Q)=D, this can be done for any vector with zero divergence. Since D=-eps*grad(V) then curl(curl(Q)/eps)=curl(grad(V))=0. In principle Q is a vector potential but, since D is in the xy plane, it only has a z component. Q is therefore in effect a scalar so this equation can be solved by FlexPDE. Contour(Q) gives the lines of D. Natural(Q)=0 puts the lines of force normal to the boundary which makes it an equipotential of V. Putting a differing value(Q) on the other facing sides of the boundary defines the flux of D as the difference between the values and therefore the charge on the plates. For the case of current flow with the current density as the vector it defines the total current flowing. If it is the potential V to be defined it is necessary to solve for V as well as Q as is done in the example. Although this leads to differential boundary conditions, so that Q is not unique, FlexPDE seems to get the right solution without problems. Deriving the equation in the example, div(grad(Q)/eps)=0, requires three further steps, and can only be done if eps is piecewise constant, i.e does not depend on Q or position, apart from changing discontinuously at region boundaries. Firstly the eps can be taken outside the bracket, secondly curl (curl) written as -div (grad) and then the eps put back in. It is not clear why this has been done as solving the curl curl equation takes exactly the same time as with div(grad). The rules about moving a material parameter across brackets are not made clear in FlexPDE documentation and, judging from the explanation in the Comsol notes, involve some rather heavy mathematics. I think you can put them where you want, provided that they stay with the group of second order differentials where they originated and are piecewise constant. In the thermal diffusion equation it seems possible to move the specific heat from the d/dt term to the div grad term and combine it with the conductivity so that only a single diffusivity is needed, but this can affect the value of the normal boundary condition. | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 1110 Registered: 06-2003 |
I don't think moving factors in and out of operators should be that obscure. We have the vector identities (a and b scalar, B a vector): grad(a*b)= a*grad(b)+ b*grad(a) div(a*B) = a*div(B) + B<dot>grad(a) curl(a*B) = a*curl(B) - B<cross>grad(a) If grad(a) is zero, you can pull it out at will. Otherwise, you need to add the grad(a) term. If the value of a is discontinuous at material boundaries, the discontinuity is handled by integration by parts in each region. You need the "a" inside the outermost div or curl to get the Divergence Theorem or Curl Theorem surface terms right. FlexPDE will convert a*div(grad(u)) to div(a*grad(u)) - grad(u)<dot>grad(a) The div() will be integrated by parts, giving a Natural BC of a*grad(u), but you will have a spurious grad(u)<dot>grad(a) term inserted in the equation. | ||
| DukeVN (dukevn) New member Username: dukevn Post Number: 2 Registered: 04-2008 |
Hi Archie and Robert, Thanks for your replies. To my understand, if div(D) is not zero (which is the case I am following), then the above method will not work? Also about the adjoint equations, it sounds interesting and perspective for solving my problem, but I have been searching around for days and I almost found nothing about adjoint equations related to E&M problems. Do you know any source for that? Thanks, |