Constraints on dependent variables in PDE Toolbox

Can one impose constrains on a dependent variable in using the PDE Toolbox? That's the general question. My specific problem is as follows. Originally the problem is the 3D Laplace equation in variable T(x,y,z). Because the geometry is a plate-like object, I finite difference in the z direction and get a system of elliptic equation in dependent variables T1(x,y), T2(x,y),... Tn(x,y). Because of additional physics going on at the top and bottom z's, T1 and Tn are described by differential equations, BUT only on part of the uppper and lower z surfaces. On the other parts, the original boundary condition (Robin) turns into contraints for T1 and Tn. Bottom line: I have n elliptic equations in T1,...,Tn. BC's on each don't matter (say zero flux). T1 and Tn have contraints in some subdomain of (x,y) of the form (T2-T1)/dz = a*T1 and similar for Tn, where dz is the distance between z levels. I tried to handle this by simply making the c coefficients for the T1 and Tn equations very small in the subset of triangles corresponding to where the contraints are. This seemingly works well, and is independent of my choice of small value for c, as long as it's small enough,, however the result is dependent on my choice of dz (i.e., number of discrete z levels I choose). Perhaps I've done the math wrong, but this in any case seems related to the more general question of imposing constraints. It would be better if there were simply a way of doing this.

[Bill Greene]

Before trying to answer this I want to make sure I understand the question. Do I understand right that you are trying to apply either a Dirichlet or Neumann constraint at points that are not on the boundary edges of the 2D model-- i.e. are at some points in the interior?

[Richard Garner]

Not quite. The original 3D equation has some boundary condition at the z surfaces (doesn't matter which really, for the purposes of this argument, but they are usually of the Robin type). When I finite difference to change this 3D equation into a system of 2D equations, the result is that the 2D equations for the variables corresponding to z=0 and z=zmax now have these boundary condtions for the original 3D problem embedded in them. However, that condition only applies to a part of their 2D geometry. The bulk of that 2D geometry is a regular 2D differential equation. To be a little more explicit, my system of n equations looks as follows: -del2(T_k) + avec*Tvec =fvec for k=1,2,...,n and where avec is an nxn vector and fvec is an nx1 vector and Tvec=[T_1,...T_n]' (nx1). However, for T_1 and T_n, the above equation is applicable to only a part of the domain. For the other part the del2 term is absent. As I mentioned in my first post, I tried to handle this by having a "c" coefficient (in the parlance of the manual) be a small number (10^-8) for the the T_1 and T_n equations in the domain where the del2 term is absent. The boundary conditions fo rthe above 2D system (at the boundary of the 2D domain in which they are defined) is simply zero flux, although I don't think this is germane to the discussion.

Actually, in retrospect, I suppose employing a constraint, which is what I think I am doing, is alot like setting internal "boundary" conditions. So in this sense the answer to your original question is "yes". Now, I think the answer to my original quesiton will be "no". (Another way to look at this is that I have a system of differntial / algebraic equations.)

Answer 1

[Bill Greene]

Hi Richard,

I'm still not sure I'm understanding you so am not sure my "solution" will help but I'll explain it anyway.

It is possible include some additional constraints in your model using the following approach. (This assumes your equations are linear and requires a bit of programming.)

You can use the matrix form of assempde:

[K,M,F,Q,G,H,R]=assempde(b,p,e,t,c,a,f)

to return the global finite element matrices for the model. The H and R matrices define a constraint equation for every boundary point where you have specifed a Dirichlet constraint. If you want to constrain a dependent degree-of-freedom at some node not on the boundary, you can simply append an equation to the H and R matrices defining this constraint. Then you pass your modified H and R matrices along with the original K, M, F, Q, and G matrices back to assempde for the solution:

u=assempde(K,M,F,Q,G,H,R)

If you decide to give this a try, I suggest trying it on a very simple model first so you can print the original and new H and R matrices to verify their correctness.

Regards,

Bill

[Richard Garner]

Bill, thanks for this information. It is very useful, in spite of my figuring out what I did wrong in my original formulation. In my original formulation, I had normalized all the equations in the system by dividing out the multipliers of all the del2 terms. Thus, all c coefficients became unity and therefore it assumed flux was simply -del(T). As you know the flux is made to be continuous internally in the domain. But, in my case the coefficients actually change internally (the coefficients are essentially thermal conductivities), but I had divided them out and incorporated them into the a and f coefficients. When I restored them into the c coefficients, the results made more sense (they agree with complete 3D finite element calculations done with Comsol Multiphyics). And, results are only minimally effected by my choice of number of z levels. There are still some little differences with Comsol, but I don't know which one is actually correct. But, thanks for your information on how to get into the interstices of (PDE Toolbox) method to set contraints. I will definitely give that a try when I get the time.

[xl z]

Here, What is the composition of the elements in the R and H matrices ? In other words, how should I modify the R and H matrices? Thank