Soluciones singulares del problema elástico antiplano en esquinas con condiciones de contorno mixtas
Full text
Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Ingeniería Aeroespacial Soluciones singulares del problema elástico antiplano en esquinas con condiciones de contorno mixtas Autor: Víctor M. Villalba Corbacho Tutores: Vladislav Mantic Lescisin Departamento de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 14 de julio de 2015
Trabajo Fin de Grado Ingeniería Aeroespacial Soluciones singulares del problema elástico antiplano en esquinas con condiciones de contorno mixtas Autor: Víctor M. Villalba Corbacho Tutores: Vladislav Mantic Lescisin Catedrático de Universidad Departamento de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 14 de julio de 2015
Trabajo Fin de Grado: Soluciones singulares del problema elástico antiplano en esquinas con condiciones de contorno mixtas Autor: Víctor M. Villalba Corbacho Tutores: Vladislav Mantic Lescisin El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha:
Índice 1. Introducción 1 1.1. La ecuación de Laplace 1 1.2. Problema Antiplano en 2-D 2 1.3. El problema de Dirichlet-Robin 3 2. Objetivos 5 3. Problema clásico de Dirichlet-Neumann 7 4. Procedimiento de Mghazli para condiciones de contorno de Robin 11 4.1. Primera descomposición del problema 11 4.2. Generalización 13 4.3. Cálculo de los términos de sombra 15 4.4. Comentarios sobre la implementación computacional 22 5. Resultados y Discusión 25 5.1. Esquina a 240◦26 5.2. Unión a tope 28 5.3. Modelo de grieta con condición de Robin 30 6. Modificación del método 33 6.1. Unión a tope 33 6.2. Modelo de grieta con condición de Robin 34 7. Conclusiones y trabajos futuros 39 Apéndice A.Códigos de Mathematica 41 Bibliografía 53 I
1 Introducción Si entendiéramos lo que estamos haciendo, no lo llamaríamos investigación. 1.1 La ecuación de Laplace Seguramente el lector esté ya familiarizado con la llamada ‘Ecuación de Laplace’, un problema en derivadas parciales, elíptico, que aparece a la hora de modelar el comportamiento de fluidos, flujos de calor, campos electromagnéticos, o, como es de mayor interés para el presente trabajo, desplazamientos en un sólido deformable. ∇2u= 0 Es bien sabido que todos los problemas en mecánica del sólido deformable son tridimensionales. Sin embargo, es posible introducir simplificaciones en los modelos matemáticos, de manera que las soluciones a problemas tridimensionales se aproximen de forma muy precisa mediante problemas bidimensionales. Explicaciones extensivas de cómo obtener casos de elasticidad plana a partir de una formulación completa en tres dimensiones para el caso de la elasticidad lineal se pueden encontrar en Teoría de la Elasticidad[3] o en Elasticity [1]. Dentro de este dominio de problemas “bidimensionales”, un caso de particular interés en ingeniería es el cálculo de las tensiones y los desplazamientos que sufre un cuerpo sometido a carga que presenta una esquina, y más concretamente una esquina reentrante, con un ángulo entre los bordes superior a 180 º , puesto que se observa en la práctica que estas geometrías suponen un decremento importante de la resistencia del sólido, al ser lo que se denomina ‘concentradores de tensiones’. Si bien a lo largo del trabajo que se desarrollará a continuación trabajaremos casi en exclusiva con un problema puramente matemático, se considera de interés el probar que el problema antiplano en elasticidad 2-D se puede reducir, en el caso de fuerzas másicas despreciables, a una ecuación de Laplace. 1
8Capítulo 3. Problema clásico de Dirichlet-Neumann Sin embargo, lo que a nosotros nos interesa aquí es desarrollar una expansión en serie de términos que cumplan, exacta o aproximadamente, las condiciones de contorno expuestas, de manera que se pueda aproximar con creciente exactitud la solución variacional del problema, a medida que aumentamos la m. Para aplicar su teorema, en un caso totalmente general, con condiciones inhomogéneas en ambas esquinas, Grisvard exige una serie de propiedades a las funciones que aparecen tanto en el lugar de nuestra γ como en los términos independientes de dichas condiciones. Se pueden consultar dichas propiedades en el artículo citado de Mghazli, y claramente se cumplen en nuestro caso, al ser γ∈C∞ , por ser una constante en nuestro caso, tomando m≥ 2, y sabiendo que la constante 0pertenece a los espacios que dice el teorema. Lo que nosotros esperamos que se pueda cumplir es la siguiente relación: u≈X 0<λj<m−1 αju(0) j Siendo αj un conjunto de constantes que habría que determinar. En la teoría de mecánica de la fractura, estos coeficientes serían los llamados "factores intensficadores de tensión generalizados". Nótese que, para un problema de fractura, los que habría que calcular aquí son factores asociados al modo III, ya que distintos coeficientes aquí no se corresponden con ningún otro modo. Para construir las funciones u(0) j a partir de los autovalores λj , hay que distinguir los casos en los que las condiciones de contorno sean del mismo tipo en ambos bordes o no. En todo momento consideramos tanto los casos mixtos como los que no lo son, por lo que podríamos tener, para este problema, condiciones de Dirichlet-Dirichlet, Dirichlet-Neumann, o Neumann-Neumann. Distinguiremos aquí las distintas situaciones por las condiciones en θ = ω , puesto que así lo hace Mghazli. Si la condición en este borde es de Dirichlet, diremos que es D, si es de Neumann o de Robin diremos que es N o R respectivamente. A su vez, si el problema tiene condiciones mixtas, diremos que es M (Mixed), mientras que en caso contrario diremos que es S (Same). Los autovalores λjse definen de la siguiente manera: λj= jπ/ω si condiciones S. (j−1 2)π/ω si condiciones M (3.1) Debemos destacar, aunque pueda parecer evidente, que si el ángulo ω se hace más pequeño, los autovalores están más dispersos en la recta real, y son mayores. Esto quiere decir que para un mismo valor de m , que sería para nuestros propósitos una especie de parámetro de calidad de la solución, tendremos menos autovalores con
9 ángulos pequeños que con ángulos grandes. Por otro lado, el aumentar excesivamente la m puede incrementar de forma muy notable el coste computacional de calcular las soluciones que se plantean en apartados posteriores. Se pueden dar dos casos que generan a su vez otra división, según lo que valga el autovalor. Si λj/∈N: u(0) j= rλjsinλjθsi condiciones D-S o N-M. rλjcosλjθsi condiciones D-M o N-S. (3.2) Si λj∈N: u(0) j= rλj[log(r)sinλjθ+θcosλjθ]si condiciones D-S o N-M. rλj[log(r)cosλjθ−θsinλjθ]si condiciones D-M o N-S. (3.3) Tendremos que comprobar a continuación que, en efecto, estas soluciones cumplen las condiciones de contorno impuestas. Los casos en que λj/∈Nson fáciles de estudiar: •Caso D-S u|θ=0 = sin(0) = 0 u|θ=ω= sin(kπ)=0 •Caso N-S 1 r ∂u ∂θ θ=0 =−sin(0) = 0 1 r ∂u ∂θ θ=ω =−sin(kπ)=0 •Caso D-M 1 r ∂u ∂θ θ=0 =−sin(0) = 0 u|θ=ω= cos((k−1/2)π)=0 Lógicamente el caso mixto en que la condición de Dirichlet se encuentra en θ = 0 es totalmente análogo, si uno considera el cambio de variable τ = ω−θ sale igualmente que se cumplen las condiciones. El caso de que un cierto autovalor sea entero es mucho más problemático, puesto que para este caso la condición de Neumann no se cumple. Este comportamiento resulta preocupante si queremos cumplir la condición de Neumann homogénea, aunque no tendría que ser preocupante si, eligiendo los coeficientes de forma adecuada, podemos aproximar otro tipo de funciones. En cualquier caso, esta problemática no la podemos abordar en este momento, de modo que seguiremos hacia adelante, puesto que si al final conseguimos cumplir bien
10 Capítulo 3. Problema clásico de Dirichlet-Neumann 0.2 0.4 0.6 0.8 1.0 r -0.4 -0.2 0.2 0.4 0.6 Neumann Homogeneous Condition Figura 3.1 Se observan en la figura tres funciones de las que propone Grisvard para autovalores enteros. Podríamos decir que a medida que la potencia aumenta, la condición se cumple mejor en un entorno cercano a la esquina, pero no es el caso para el primer término, que es puramente lineal. las condiciones de contorno habremos resuelto el problema desde un punto de vista práctico.
4 Procedimiento de Mghazli para condiciones de contorno de Robin 4.1 Primera descomposición del problema El trabajo de Mghazli supone una extensión del de Grisvard, hacia el problema de Dirichlet-Robin que recordamos aquí: ∇2u= 0 u|θ=0 = 0 1 r ∂u ∂θ +γuθ=ω = 0 Ella consigue demostrar que este problema posee la misma propiedad que el de Dirichlet-Neumann, de manera que se puede obtener una descomposición que permita afirmar que la solución variacional está en un espacio de Sobolev de orden m. Pensemos que podemos aproximar la solución al problema de Dirichlet-Robin como la correspondiente al de Dirichlet-Neumann según el método de Grisvard: u≈X 0<λj<m−1 αju(0) j ∂u ∂θ ≈X 0<λj<m−1 αj ∂u(0) j ∂θ Entonces, la condición de Robin se puede escribir como: ∂u ∂θ θ=0 =−γr X 0<λj<m−1 αju(0) j|θ=0 Si ahora asumimos que esta relación se cumple término a término, para cada valor de j , se puede plantear una serie de j problemas en derivadas parciales con condiciones 11
12 Capítulo 4. Procedimiento de Mghazli para condiciones de contorno de Robin de Dirichlet-Neumann, pero no homogéneas: ∇2u(1) j= 0 u(1) j|θ=0 = 0 ∂u(1) j ∂θ θ=ω =−γrλj+1 sinλjω Nos ocuparemos de la solución de este problema en derivadas parciales en la siguiente sección del documento, con el propósito de tratar en detalle un caso más general, que es el correspondiente a la solución para cualquier problema de este tipo en que el coeficiente multiplicando a la potencia de rsea una constante. Lo que esperamos es que el término solución a este problema mejore las propiedades de la solución, como manera de tener en cuenta que la primera solución presenta una cierta diferencia con la real. . Figura 4.1 La línea azul representa el error en el cumplimiento de la condición de Robin en la solución de Grisvard para el primer autovalor calculado con un ángulo de 240 º con condiciones de Dirichlet-Robin sobre el borde de Robin. La amarilla representa lo mismo cuando se le suma la solución al problema anterior. Lo que se ha hecho para obtener la gráfica anterior es evaluar numéricamente la condición de Robin en el borde correspondiente, para una γ= 1. La expresión concreta del término, que deduciremos a posteriori, es: u(1) j= a(1) jrλj+1 sin(λj+1)θsi ω /∈π,2π a(1) jrλj+1[sin(λj+ 1)θlogr+θcos(λj+1)θ]si ω∈π,2π(4.1)
4.2 Generalización 13 Con a(1) jdefinido como: a(1) j= −γ (λj+1)sinωsi ω /∈π,2π γ ω(λj+1)cosωsi ω∈π,2π(4.2) Lo que nos dice esta gráfica, es que el añadir este término de superíndice (1) "mejora" el cumplimiento de la condición de Robin en el borde. En lo que sigue denominaremos a estos términos u(l) j como ’términos de sombra de orden ló ’shadow terms’, y cálculo constituirá una gran parte del objetivo del trabajo. A pesar del nombre, los términos de sombra no son, en general un solo término, siendo en general un sumatorio de funciones que se calculan a continuación. La idea es que tomar suficientes términos de este tipo consigue adaptar la solución para que, en un cierto rango de valores de r, la condición de Robin se cumpa razonablemente bien, y asumimos que la solución no es válida cuando estamos lejos del origen. Esto coincide con el enfoque clásico en tensiones que se utiliza ampliamente en mecánica de la fractura, con soluciones que presentan una singularidad en el labio de la grieta, y se consideran inválidas en distancias grandes. Un ejemplo clásico de este tipo de solución sería el campo de Westergaard en Mecánica de la Fractura Elástico-Lineal, si bien esta es válida para el problema plano, regido por una ecuación biharmónica. La importancia de obtener buenos resultados en el caso de ω = 180 ◦ radica en que es uno de los que introduce términos en logaritmos de r. A la vista de los resultados, al menos podemos afirmar que podemos aplicar este método para aproximar la solución al problema. Falta extenderlo para ver si podemos obtener soluciones que sean cada vez mejores a costa solamente de un mayor tiempo de computación. 4.2 Generalización Puede que lo primero que pensase uno a la vista de los resultados del apartado anterior fuese plantear inmediatamente otro problema en derivadas parciales por el mismo procedimiento para todos los autovalores λj para tratar de adaptar la solución aún más. Sin embargo, no es eso exactamente lo que hace Mghazli. Se plantea ahora la siguiente igualdad: u≈X 0<λj<m−1 αjΦj
14 Capítulo 4. Procedimiento de Mghazli para condiciones de contorno de Robin . Figura 4.2 La línea azul representa la condición de Robin para el caso en que ω = 180 ◦ usando solamente el término de Grisvard, añadiendo para la amarilla el primer término de sombra. Los términos Φson sumatorios de funciones con la forma siguiente: Φj= m−1−hj X l=0 hj−1<λj≤hj λj6=m−1 u(l) j Entender este sumatorio es extremadamente importante a la hora de aplicar el método. Cada uno de los términos Φ j es un sumatorio sobre l con el límite superior delimitado por la variable hj. La hj es el número natural más pequeño tal que λj≤hj . Esto es, el valor de λj determina el valor de hj. La introducción de esta variable en el sumatorio tiene el efecto de forzar la convergencia de la serie, puesto que a medida que aumenta λj , y con ello hj el límite superior va reduciéndose hasta 0, acabando con la suma. Desde el punto de vista de nuestros propósitos, posiblemente esto sea restarle "calidad" a las soluciones que queremos encontrar, al quitar términos que en principio nos ayudarían a aproximar, pero puesto que Mghazli lo definió así, elegimos dejarlo. Para entender por qué, tendremos que ver cual es la forma de los términos de sombra: Se considerarán dos casos, en que las condiciones sean de Dirichlet-Robin, y de Robin-Robin, puesto que el problema con condiciones de Dirichlet-Dirichlet estaría
4.3 Cálculo de los términos de sombra 15 cubierto por el teorema de Grisvard y el de Robin-Dirichlet sería igual al de DirichletRobin: •Condiciones de Dirichlet-Neumann ∇2u(l+1) j= 0 u(l+1) j|θ=0 = 0 ∂u(l+1) j ∂θ θ=ω =−γru(l) j|θ=ω •Condiciones de Neumann-Neumann ∇2u(l+1) j= 0 ∂u(l+1) j ∂θ θ=0 =γ1ru(l) j|θ=0 ∂u(l+1) j ∂θ θ=ω =−γ2ru(l) j|θ=ω Donde se ha tenido en cuenta que para la condición de contorno en θ = ω la normal exterior es opuesta al avance de la coordenada θ. El cálculo de los términos de sombra se plantea por tanto de forma recursiva: Cada uno de ellos se calcula a partir del anterior, siendo el primero de todos el correspondiente al teorema de Grisvard. 4.3 Cálculo de los términos de sombra Empezaremos aquí por analizar el caso más simple, que correspondería al cálculo del primer término de sombra y todos aquellos que no presentasen comportamiento logarítmico, y demostraría que pueden aparecer términos logaritmicos incluso si el autovalor del que emanan no es entero. Empezaremos definiendo la constante G1 como el coeficiente que multiplica a la potencia de r en la condición de Neumann, de momento. Veremos inmediatamente que las constantes en realidad son coeficientes de un polinomio en logaritmos. Estas Gi se definen solo para el cálculo de cada término de sombra correspondiente para simplificar la notación, y no aparecen explícitamente en las soluciones finales.
16 Capítulo 4. Procedimiento de Mghazli para condiciones de contorno de Robin ∇2u(l+1) j= 0 u(l+1) j|θ=0 = 0 ∂u(l+1) j ∂θ θ=ω =G1rλj+l+1 Donde se tiene una r de la imposición de la condición de Neumann y una rλj+l de particularizar u(l) jen el borde correspondiente. Es fácil ver que, para el caso del primer término de sombra que dejamos pendiente en el apartado anterior: rλj+1G1=−γu(0) j|θ=ω La solución más inmediata consistiría en buscar soluciones de la forma rλj+l+1 Γ( θ ), pero para ciertos valores de ω no existe solución de esta forma. Esto pasa si cos(λj+l+1)ω= 0. Las soluciones serán pues del tipo: u(l+1) j=rλj+l+1[Γ1(θ)logr+Γ2(θ)] Aplicando el operador de Laplace en coordenadas polares (definido en la introducción), se obtienen dos sistemas de ecuaciones diferenciales ordinarias: Γ00 1+(λj+l+ 1)2Γ1= 0 Γ1(0) = 0 Γ0 1(ω)=0 (4.3) Γ00 2+(λj+l+ 1)2Γ2=−2(λj+l+1)Γ1 Γ2(0) = 0 Γ0 1(ω) = G1 (4.4) La resolución a este sistema distingue los casos en que cos ( λj + l + 1) ω = 0 o cos(λj+l+1)ω6= 0. Finalmente obtenemos los siguientes resultados: Si cos(λj+l+1)ω6= 0: u(l+1) j=rλj+l+1a(l+1) jsin(λj+l+1)θ Si cos(λj+l+1)ω= 0: u(l+1) j=a(l+1) j[sin(λj+l+1)θlog(r)+θcos(λj+l+1)θ]
4.3 Cálculo de los términos de sombra 17 Con las constantes a(l+1) jdadas por: a(l+1) j= G1 (λj+l+1)cos(λj+l+1)ωsi cos(λj+l+1)ω6= 0 −G1 ω(λj+l+1)sin(λj+l+1)ωsi cos(λj+l+1)ω= 0 (4.5) Donde se ha tomado la constante que no queda fijada por las condiciones de contorno como nula. Si se particularizan los valores correspondientes al cálculo del primer término de sombra, se obtienen los resultados ya presentados anteriormente. El resultado implica que, incluso aunque en el término de sombra anterior no aparezca explícitamente un logaritmo, es posible que haya uno en el siguiente. En general, se puede comprobar que para el caso en que aparezca un logaritmo la condición de Neumann del problema que hay que resolver, si no consideramos por un momento el término no logarítmico, es de la forma: ∂f ∂θ θ=ω =G2rβ+1log(r) Las soluciones que se obtienen para este caso, que es también el caso en que el término de Grisvard tenga tabién un logaritmo, realizando una separación de variables similar a la anterior, son del tipo: f=rβ+1[Γ1(θ)(log(r))2+Γ2(θ)log(r)+Γ3(θ)] No se va a resolver explícitamente este caso, puesto que en general tendremos un término logarítmico y un término que no lo es, multiplicando a la potencia de r . El propósito es ilustrar que la presencia de logaritmos puede inducir la aparición de otros elevados a potencias superiores a medida que uno calcula más términos de sombra. Esto ocurre en general para cualquier orden del término de sombra, como puede comprobar el lector si ejecuta el algoritmo de Wolfram Mathematica correspondiente. Esto nos lleva al caso más general de problema en derivadas parciales para calcular un término de sombra: ∇2u(l+1) j= 0 u(l+1) j|θ=0 = 0 ∂u(l+1) j ∂θ θ=ω =rλj+l+1"n X i=0 Gi(log(r))i# Donde el índice n representa el máximo exponente al que aparece elevado un término del tipo log ( r ), que en general no tiene por qué tener una relación directa con el orden l más allá de que esta es la cota superior de lo que puede valer. No podemos expresar n mediante una fórmula directamente, pero podemos calcularla observando
5 Resultados y Discusión En esta sección se aplicará el método presentado con anterioridad a la resolución de dos problemas elásticos con interés particular: la unión a tope o "Butt-joint" y la grieta con borde cohesivo. También presentaremos, como preámbulo, la solución a un ángulo de 240 ◦ con condiciones de Dirichlet-Robin, para comentar efectos y propiedades que no son evidentes a simple vista, y que ayudarán al lector, esperemos, a entender mejor cómo se presentan los resultados. Se podría haber elegido otro ángulo, se consideró este al ser mayor que 180 ◦ , puesto que la esquina reentrante es la que, en la práctica, genera problemas estructurales por concentración de tensiones. El primer problema "especial" es una placa con un ángulo de 90 ◦ unida a otra. Esto quiere decir que se presenta un borde libre de tensiones y otro unido a un sólido con propiedades cualesquiera. Se supone una interfase compuesta por una distribución continua de muelles ideales de rigidez k y longitud natural nula. Traducido al método que presentamos aquí, esto nos daría una esquina con condiciones de Robin-Robin, con γ1=kyγ2= 0. El segundo problema realmente sería un modelo de "mitad de grieta", en que nos quedamos con la parte de arriba, y suponemos que la grieta está libre de tensiones, mientras que el borde sigue siendo un sólido continuo, pero dañado. La ley cohesiva sigue siendo la misma que antes, lo cual nos daría condiciones Robin-Robin con un ángulo de 180◦y una de las γnula, con la otra teniendo valor k. Curiosamente, ambos casos son de los que podemos llamar "problemáticos". No solo porque aparecen términos logarítmicos, además son ángulos que, según los valores de m , pueden quedar excluidos de la descomposición de Mghazli. En la práctica, vamos a ignorar esos casos, y a estudiar si cumplen la condición de Robin, aunque puedan no estar incluidos en los teoremas de Grisvard o Mghazli. Sin embargo, esto es un motivo para probar también la esquina a 240 ◦ con condiciones de RobinRobin, para verificar que el método es bueno aunque pueda fallar en casos particulares. Se informa al lector de que, a la hora de comprobar el cumplimiento o no de la condición de Robin, se presentará el error evaluado en el borde correspondiente, junto con el valor de cada aproximación, Φ. El motivo es que si nuestro método no 25
26 Capítulo 5. Resultados y Discusión va a dar resultados buenos, el error de la condición de Robin será: 1 r ∂u ∂θ +γuθ=ω|0 ≈O(u|θ=ω|0) De forma que, si la gráfica del error está muy por debajo en valor absoluto de la de la Φcorrespondiente, tendremos una buena aproximación a la solución y daremos el método por válido. 5.1 Esquina a 240◦ Para este caso se ha utilizado m = 8, aunque muy posiblemente habríamos sido capaces de poner un grado más alto. Se obtienen pues los autovalores siguientes: 3 8,9 8,15 8,21 8,27 8,33 8,39 8,45 8,51 8 Vamos ahora a estudiar si en efecto las funciones solución cumplen las condiciones de contorno, aunque no las vamos a escribir explícitamente. Empezando para la de Dirichlet: ω= 0. Las funciones Φ, para este caso, están todas compuestas de funciones seno multiplicadas por funciones de r, de modo que, al ser sin (0) = 0 la condición se cumple de forma exacta para todas ellas. Figura 5.1 Se observa aquí el error de la condición de Robin para tres de las funciones Φcalculadas con m=8, con condiciones de Dirichlet-Robin. Se puede ver que el error es pequeño, puesto que en el eje de ordenadas aparecen magnitudes del orden de 10−6. Por ser el más importante, vamos a estudiar en detalle la primera de las funciones, que es suma de un término singular más los términos de sombra.
5.1 Esquina a 240◦27 Figura 5.2 Aquí vemos el comportamiento esperado, el sumar términos de sombra mejora el cumplimiento de la condición de Robin hasta hacer el error prácticamente nulo. La línea azul representa solo el primer término, la amarilla aparece al sumar un término de sombra más, y ya con cuatro términos parece que cumplimos de forma prácticamente exacta.. .
28 Capítulo 5. Resultados y Discusión Vamos a probar ahora lo mismo, con la misma esquina, para las condiciones de Robin-Robin, y veremos el cumplimiento de las condiciones en ambos bordes de la esquina: Figura 5.3 Representamos el error en la esquina θ = 0 calculada con m = 8. Tenemos una serie de potencias que van aumentando a medida que se alejan del origen.. Figura 5.4 Aquí también tenemos el error en el cumplimiento de la condición de Robin. Como se ve obtenemos un numérico al particularizar para θ = ω , por lo que se puede afirmar que cumplimos. 5.2 Unión a tope Para este cálculo se utilizará un valor m = 8, por ser este el límite que nos impone nuestra potencia de cálculo, y se calcularán las funciones Phi correspondientes. A continuación, se calculará el error del cumplimiento en la condición de Robin. Los autovalores correspondientes a este problema resultan ser: λ={2,4,6,8}
5.2 Unión a tope 29 Y se obtienen las siguientes funciones Φ: Φ1=r4(12πcos(4θ)log(r)−12πθsin(4θ)) 144π+r31 3(cos(3θ)log(r)−θsin(3θ))−1 9cos(3θ) +r2(cos(2θ)log(r)−θsin(2θ))+O(r5) (5.1) Φ2=r6(30πcos(6θ)log(r)−30πθsin(6θ)) 900π+r51 5(cos(5θ)log(r)−θsin(5θ))−1 25 cos(5θ) +r4(cos(4θ)log(r)−θsin(4θ))+O(r7) (5.2) Φ3=r71 7(cos(7θ)log(r)−θsin(7θ))−1 49 cos(7θ)+r6(cos(6θ)log(r)−θsin(6θ)) (5.3) Existen más términos, que el código de Mathematica extrae sin problemas. Aunque no se vea a simple vista, la condición homogénea de Neumann se cumple estrictamente al derivar con respecto de θy particularizar en θ= 0. Vamos pues con la condición de Robin. Si realizamos la misma operación, las expresiones para el error sobre la condición de Robin son: e1=π r7 5040 +r!(5.4) e2=1 420πr3r4−840(5.5) e3=1 14πr5r2+42(5.6) Que obviamente no da 0. De hecho, la condición de Robin no se cumple bien para la primera de estas funciones, ya que presenta un término lineal que se desprende inmediatamente del eje de abscisas. Vamos a ver más de cerca el resultado:
30 Capítulo 5. Resultados y Discusión Figura 5.5 Si bien todas las funciones correspondientes a autovalores distintos de 2 cumplen bien las condiciones en un cierto intervalo, el primer autovalor de la serie incumple mediante una función lineal. 5.3 Modelo de grieta con condición de Robin Para un m = 6 obtenemos los siguientes resultados de error en el cumplimiento de la condición de Robin: Figura 5.6 Se representa el error el la condición de Robin para las Φque salen con m = 6 en la cara cohesiva ( θ = 0). El comportamiento aquí se adapta bastante bien en un pequeño intervalo, pero se acaba despegando del eje, y no se observa una buena adaptación como en otros casos.. Se observa pues que el método, al menos tal como está planteado, no da buenos resultados, puesto que no nos cumple las condiciones que queremos que cumpla. Esto no quiere decir que Mghazli se equivocase, seguramente consiguiese sus objetivos de estudio de regularidad pero sí que intentamos extender su trabajo directamente a un dominio para el que no estaba pensado, y no ha dado resultado. Sin embargo, una modificación del método podría compensar esto.
5.3 Modelo de grieta con condición de Robin 31 Figura 5.7 Aquí se representa el error en la cara libre. Se vuelve a incumplir, al igual que en el problema anterior, una de las condiciones de contorno de forma fuerte por la presencia de una recta asociada al autovalor 2.. Escribimos a continuación las funciones encontradas. Φ1=−r9(9(π−2θ)cos(9θ) +2sin(9θ)(1−9log(r))) 3265920 +r8((π−2θ)sin(8θ)+2cos(8θ)log(r)) 40320 −r7(7(π−2θ)cos(7θ) +2sin(7θ)(1−7log(r))) 35280 +1 720r6((π−2θ)sin(6θ)+2cos(6θ)log(r)) −1 600r5(5(π−2θ)cos(5θ)+2sin(5θ)(1−5log(r)))+ 1 24r4((π−2θ)sin(4θ)+2cos(4θ)log(r)) −1 18r3(3(π−2θ)cos(3θ)+ 2sin(3θ)(1 −3log(r)))+r2(cos(2θ)log(r)−θsin(2θ)) Φ2=−r9(9(π−2θ)cos(9θ) +2sin(9θ)(1−9log(r))) 272160 +r8((π−2θ)sin(8θ)+2cos(8θ)log(r)) 3360 −r7(7(π−2θ)cos(7θ) +2sin(7θ)(1−7log(r))) 2940 +1 60r6((π−2θ)sin(6θ)+2cos(6θ)log(r)) −1 50r5(5(π−2θ)cos(5θ)+ 2sin(5θ)(1 −5log(r)))+r4(cos(4θ)log(r)−θsin(4θ)) Φ3=−r9(9(π−2θ)cos(9θ) +2sin(9θ)(1−9log(r))) 9072 +1 112r8((π−2θ)sin(8θ)+2cos(8θ)log(r)) −1 98r7(7(π−2θ)cos(7θ)+ 2sin(7θ)(1 −7log(r)))+r6(cos(6θ)log(r)−θsin(6θ)) Φ4=1 162r8(2(−81θsin(8θ)+rsin(9θ)(9log(r)−1)+81cos(8θ)log(r))−9(π−2θ)rcos(9θ))
6 Modificación del método Para poder cumplir las condiciones de contorno, el profesor Mantic propuso una modificación de las ecuaciones de Mghazli, basada en que los términos del desarrollo de Grisvard para autovalores fraccionarios cumple también las condiciones en el caso de que sean enteros. El cambio consiste en ignorar la diferencia entre autovalores enteros y fraccionarios, y montar así un conjunto de primeros términos para el desarrollo, a partir de los cuales se calcularán los términos de sombra según el algoritmo de Mghazli. Esto quiere decir que según los valores de ω tendremos o no términos logarítmicos, pero lo que es seguro es que estos no aparecerán en el primer término. Así pues, los primeros términos del desarrollo serían: u(0) j= rλjsinλjθsi condiciones D-S o N-M. rλjcosλjθsi condiciones D-M o N-S. (6.1) Y no distinguimos casos según los valores λj. Lo único que se hará mediante este método es calcular las soluciones para los dos casos particulares que se han planteado: butt-crack joint y modelo de grieta. 6.1 Unión a tope Las funciones obtenidas en el primer caso para m= 8 son: Φ1=1 3r3cos(3θ)+ r2cos(2θ)(6.2) Φ2=1 5r5cos(5θ)+ r4cos(4θ)(6.3) Φ3=1 7r7cos(7θ)+ r6cos(6θ)(6.4) 33
40 Capítulo 7. Conclusiones y trabajos futuros Para concluir, si bien el problema antiplano que se ha estudiado aquí tiene su interés en la práctica, es indiscutible que el problema que verdaderamente uno querría resolver, tratando elasticidad bidimensional, es el problema plano, gobernado por la Ecuación Biarmónica. ∇4φ= 0 (7.1) No se ha hecho ningún estudio sobre la posibilidad de exportar este método a este problema, pero sería muy interesante obtener funciones aproximantes que permitiesen obtener cualquier grado de precisión para esquinas trabajando en este régimen.
Apéndice A Códigos de Mathematica 41
42 Apéndice A. Códigos de Mathematica (*This is the unmodified code for the DirichletRobin problem. It needs a .txt file as input and with the exact structure presented in the example file, although some of the inputs are unused and can be removed.*) ClearAll[m, θi, k1, k2, bbcc, ω, nterms, λ, j, a1, u0, u1, r, θ, u, a, h, A, G, B, Φ,Ω,β,γ1, γ2, Γ] (*Reading block, in this section of the code we should be able to extract the data for use: the regularity of the solution m, a vector with the corner's angles θi(two components), three constants with the material's properties and another vector containing the boundary conditions, bbcc.*) SetDirectory["C:\Cosas de Ingenieros\TFG\Ejecutar"]; infile =InputString["Input file?"]; input =ReadList[infile, Number, RecordLists →True]; m=input[[1, 1]]; θi=input[[2]] * π/180; (*γ1=input[[3,1]];γ2=input[[3,2]];*) (*Activate this line if you wish to use γ1 and γ2 as datafile input.*) bbcc =input[[4]]; ω=θi[[1]] - θi[[2]]; ωaux =ω; If[bbcc[[1]] ⩵bbcc[[2]], nterms =IntegerPart[ω/π* (m-1)], nterms =IntegerPart[ω/π* (m-1)+1/2] ] (*θi is called like that to avoid confusion with the polar angle*) (*Definition of the problem module, here we need to identify the kind of problem we have introduced with the datafile *) (*WARNING: THESE EQUATIONS ARE WRITTEN UNDER MGHAZLI'S POLAR COORDINATES, WITH ORIGIN ON ONE OF THE BOUNDARIES OF THE WEDGES.*) (*Table of groups: M: Mixed boundary conditions. S: Same boundary conditions R:Robin condition D: Dirichlett condition N: Neumann condition *) (*The philosophy now is as follows: Once the datafile has been read and the bbcc vector has been created, it is compared in succession with different sets of boundary conditions. If the bbcc vector doesn't fulfill the conditions, it is compared with the next until the correct one is found. *) j; If[bbcc[[1]] ⩵bbcc[[2]], λ=Table[j*Pi /ω,{j, nterms}], λ=Table[(j-1/2)*Pi /ω,{j, nterms}]] h=Table[Ceiling[λ[[j]]],{j, nterms}] If[bbcc[[1]] ⩵2 && bbcc[[2]] == 0, (*For this situation, Mghazli's notation would give j contained in M,R*) (*As for "our" indexes, j runs from 1 to the dimension of λ, while l runs from 0 to m-h-2, in order to give shadow terms up to m-h-1. Thus, l needs an offset in order to be put in Mathematica's lists, while j doesn't.*)
43 Φ=Table[0, {j, Dimensions[λ][[1]]}]; u=Table[0, {i, 1, Dimensions[λ][[1]]},{k, 1, m}]; Do[ If[! Element[λ[[j]], Integers], u[[j, 1]]=r^(λ[[j]]) * Sin[λ[[j]]*θ], u[[j, 1]] = r^λ[[j]] * (Log[r]Sin[λ[[j]] θ]+θCos[λ[[j]] θ])] ,{j, 1, Dimensions[λ][[1]]}] Do[ Do[ (*Within this loop, we are supposed to find each one of the singular functions or "shadow terms" associated with the problem for a "frozen" λ_k. The process is as follows: Each PDE results in the shadow term u_l, and is in turn given as a particularization for θ=ωof the prior shadow term u_(l-1)The parameters we need to define for the problem are β, n, and G. The G coefficients are those that multiply each one of the Log[r]terms in the prior shadow term.*) (*The general expression for the Neumann nonhomogeneous boundary condition is - r*γ*u_(l-1). Obviously each shadow term can be regarded as a polynomial in Log[r]multiplied by a power of r, the power being r^(λ_k+l-1).*) (*Thus, it is trivial to see that the shadow term u_l has a power superior to the one in u_(l-1)by one.*) (*A very important remark is that the Log[r]terms won't even exist in the first place if λ_k is not an integer, since the shadow terms are derived in succession, if the first shadow term doesn't have a logarithm, those that follow will not present logarithmic behaviour.*) (*Remark: The proof of the global theorem is valid for q=0, 1,2...,m-h-2. I don't have to worry about the upper limit, since q=m-h-2 would give a PDE for which the solution is the shadow term of order m-h-1. *) β=λ[[k1]] + l+1-1; (*This expression matches the exponent of r in the complete problem with that of Mghazli's appendix.*) θ=ω; Ω= -r*γ1*u[[k1, l +1]]; ClearAll[θ]; (*We know the power of the u_l shadow term is λ_k+l, so if we divide Ωby r^(λ_k+l+1)we have a pure polynomial in Log[r] (if there are logarithms at all)*) G=CoefficientList[Ω/r^(β+1), Log[r]]; n=Dimensions[G, 1][[1]]; Print["G=", G]; Print["DimensionG=", Dimensions[G]]; Print["n=", n]; Print["β=", β]; Print["ω=", ω]; Print["Ω=", Ω]; A=Table[0, {i, 1, n +1}]; Ct = Table[If[k⩵1, 1, Product[(i-1+l),{l, 1, k -1}] / (k-1)! / (β+1)^(k-1)], {i, 1, n +1},{k, 1, n +1}]; Print[Ct]; 2
44 Apéndice A. Códigos de Mathematica Print[ω]; ClearAll[ω]; If[Cos[(β+1)*ωaux]≠0 && Dimensions[G][[1]] ≠0, A[[n+1]]=0; A[[n+1-1]] = G[[n+1-1]] / (β+1)/Cos[(β+1)ω]; Do[ A[[i+1]]=1/ (β+1)/Cos[(β+1)ω] * (G[[i+1]] - Sum[A[[i+1+k]] * Ct[[i+1, k +1]] * ω^(k-1)*(k*D[Sin[(β+1) * ω], {ω, k}] + ω*D[Sin[(β+1) * ω],{ω, k +1}]),{k, 1, n -i}]); ,{i, n -2, 0, -1}]; ω=ωaux; u[[k1, l +1+1]] = r^(β+1)*(Sum[A[[j+1]] * Sum[Ct[[j+1-k, k +1]]*θ^k * D[Sin[(β+1) * θ] * Log[r]^(j-k),{θ, k}],{k, 0, j}],{j, 0, n}]); ]; If[Cos[(β+1)*ωaux] == 0 && Dimensions[G][[1]] ≠0, A[[1]]=0; A[[n+1]] = -G[[n+1-1]]/ω/ (β+1)/n/Sin[(β+1) * ω]; Do[A[[i+1]] = 1/ω/i/ (β+1) / Sin[(β+1)*ω] * (-G[[i+1-1]] + Sum[A[[i+1-1+k]] * Ct[[i+1-1, k +1]] * ω^(k-1) * (k*D[Sin[(β+1) * ω],{ω, k}]+ω*D[Sin[(β+1) * ω],{ω, k +1}]), {k, 2, n -i+1}]),{i, n -1, 1, -1}]; ω=ωaux; u[[k1, l +1+1]] = r^(β+1)*(Sum[A[[j+1]] * Sum[Ct[[j+1-k, k +1]]*θ^k * D[Sin[(β+1) * θ] * Log[r]^(j-k),{θ, k}],{k, 0, j}],{j, 0, n}]); ]; If[Dimensions[G][[1]] ⩵0, u[[k1, l +1+1]] = 0]; Print["ωdeb=", ω]; Print["A=", A]; Print["u_kl=", u[[k1, l +1+1]]]; Print["k l=", k l] Print["Iteration complete"]; ,{l, 0, m -2-h[[k1]]}]; Φ[[k1]] = Sum[u[[k1, j1]],{j1, 1, Dimensions[λ][[1]]}] ,{k1, 1, Dimensions[λ][[1]]}] ] If[bbcc[[1]] ⩵2 && bbcc[[2]] ⩵2, (*For this situation, Mghazli's notation would give j contained in M,R*) (*As for "our" indexes, j runs from 1 to the dimension of λ, while l runs from 0 to m-h-2, in order to give shadow terms up to m-h-1. Thus, l needs an offset in order to be put in Mathematica's lists, while j doesn't.*) Φ=Table[0, {j, Dimensions[λ][[1]]}]; u=Table[0, {i, 1, Dimensions[λ][[1]]},{k, 1, m}]; Do[ If[! Element[λ[[j]], Integers], u[[j, 1]]=r^(λ[[j]]) * Cos[λ[[j]]*θ], u[[j, 1]] = r^λ[[j]] * (Log[r]Cos[λ[[j]] θ]-θSin[λ[[j]] θ])] ,{j, 1, Dimensions[λ][[1]]}]; Do[ Do[ 3
45 (*Within this loop, we are supposed to find each one of the singular functions or "shadow terms" associated with the problem for a "frozen" λ_k. The process is as follows: Each PDE results in the shadow term u_l, and is in turn given as a particularization for θ=ωof the prior shadow term u_(l-1)The parameters we need to define for the problem are β, n, and G. The G coefficients are those that multiply each one of the Log[r]terms in the prior shadow term.*) (*The general expression for the Neumann nonhomogeneous boundary condition is - r*γ*u_(l-1). Obviously each shadow term can be regarded as a polynomial in Log[r]multiplied by a power of r, the power being r^(λ_k+l-1).*) (*Thus, it is trivial to see that the shadow term u_l has a power superior to the one in u_(l-1)by one.*) (*A very important remark is that the Log[r]terms won't even exist in the first place if λ_k is not an integer, since the shadow terms are derived in succession, if the first shadow term doesn't have a logarithm, those that follow will not present logarithmic behaviour.*) (*Remark: The proof of the global theorem is valid for q=0, 1,2...,m-h-2. I don't have to worry about the upper limit, since q=m-h-2 would give a PDE for which the solution is the shadow term of order m-h-1. *) β=λ[[k1]] + l+1-1; (*This expression matches the exponent of r in the complete problem with that of Mghazli's appendix.*) θ=0; Ω=r*γ1*u[[k1, l +1]]; (*Both Ωand Γrepresent the value of the Neumann condition, each in its corresponding edge. We can then extract the H coefficients, which determine the PDE we solve.*) ClearAll[θ]; θ=ω; Γ= -r*γ2*u[[k1, l +1]]; ClearAll[θ]; (*Print["Ω="Ω]; Print["Γ="Γ];*) (*We know the power of the u_l shadow term is λ_k+l, so if we divide Ωby r^(λ_k+l+1)we have a pure polynomial in Log[r] (if there are logarithms at all)*) H= {CoefficientList[Ω/r^(β+1), Log[r]], CoefficientList[Γ/r^(β+1), Log[r]]}; n=Dimensions[H][[2]]; Print["H=", H]; Print["DimensionH=", Dimensions[H]]; Print["n=", n]; Print["β=", β]; Print["ω=", ω]; Print["Ω=", Ω]; Print["Sin ", Sin[(β+1) * ω]]; A=Table[0, {i, 1, n +1}]; B=Table[0, {i, 1, n +1}]; (*The column [[:,1]] of H represents the Robin condition for θ=0, while [[:,2]] does likewise for θ=ω*) 4
46 Apéndice A. Códigos de Mathematica Ct =Table[If[k⩵1, 1, Product[(i-1+l),{l, 1, k -1}] / (k-1)! / (β+1)^(k-1)],{i, 1, n +1},{k, 1, n +1}]; (*We calculate the Ct constants in the inner loop because we need to guarantee they reach n+1 in both dimensions, and only within the loop we know n.*) i=1; Print["Ct", Ct]; ClearAll[ω]; If[Sin[(β+1)*ωaux]≠0 && Dimensions[H][[2]] ≠0, A[[n+1]]=0; Do[ A[[i+1]]=H[[1, i +1]] / (β+1) - (i+1) / (β+1) * A[[i+1+1]], {i, n -1, 0, -1}]; Print[A]; B[[n+1]]=0; B[[n+1-1]] = -H[[2, n +1-1]] / (β+1) / Sin[(β+1)*ω] + A[[n+1-1]] * Cot[(β+1) * ω]; Do[ B[[i+1]] = 1/ (β+1)/Sin[(β+1)*ω] * (-H[[2, i +1]] + Sum[Ct[[i+1, k +1]] * (k*ω^ (k-1) * (A[[i+1+k]]*D[Sin[(β+1) * ω],{ω, k}]+B[[i+1+k]] * D[Cos[(β+1) * ω],{ω, k}]) + ω^k * (A[[i+1+k]]*D[ Sin[(β+1) * ω],{ω, k +1}] + B[[i+1+k]] * D[Cos[(β+1)*ω], {ω, k +1}])),{k, 1, n -i}]) + A[[i+1]]*Cot[(β+1)*ω]; ,{i, n -2, 0, -1}]; ω=ωaux; Print[B]; u[[k1, l +1+1]] = r^(β+1) * Sum[A[[j+1]]*Sum[Ct[[j+1-k, k +1]] * θ^k *D[Sin[(β+1)*θ],{θ, k}] * Log[r]^(j-k),{k, 0, j}]+B[[j+1]] * Sum[Ct[[j+1-k, k +1]] * θ^k * D[Cos[(β+1) * θ],{θ, k}]*Log[r]^(j-k),{k, 0, j}],{j, 0, n}]]; If[Sin[(β+1)*ωaux]⩵0 && Dimensions[H][[2]] ≠0, (*This applies to angles for which Sin[(β+1)*ω]=0*) A[[n+1]]=0; Do[ A[[i+1]] = -H[[1, i +1]] / (β+1)-(i+1)/(β+1) * A[[i+1+1]], {i, n -1, 0, -1}]; Print[A]; B[[n+1]] = -H[[2, n +1-1]] / (β+1) / Cos[(β+1)*ω]/n/ω+ A[[n+1-1]] / n/ω;(*Warning Hn isn't defined, since the boundary conditions only have H till order n-1, so Mghazli has a mistake*) B[[1]] 0; Print[B]; Do[B[[i+1]] = 1/i/ω/ (β+1)/Cos[(β+1) * ω] * (-H[[2, i +1-1]] + Sum[Ct[[i+1-1, k +1]] * (k*ω^(k-1) * (A[[i+1-1+k]]*D[Sin[(β+1) * ω],{ω, k}] + B[[i+1-1+k]]*D[Cos[(β+1)*ω],{ω, k}]) + ω^k * (A[[i+1-1+k]] * D[Sin[(β+1) * ω],{ω, k +1}]+B[[i+ 1-1+k]]*D[Cos[(β+1)*ω],{ω, k +1}])),{k, 2, n -i+1}]) + A[[i+1]]/ω/ (β+1) + A[[i+1-1]]/i/ω,{i, n -1, 1, -1}]; ω=ωaux; 5
47 u[[k1, l +1+1]] = r^(β+1) * Sum[A[[j+1]]*Sum[Ct[[j+1-k, k +1]] * θ^k *D[Sin[(β+1)*θ],{θ, k}] * Log[r]^(j-k),{k, 0, j}]+B[[j+1]] * Sum[Ct[[j+1-k, k +1]] * θ^k * D[Cos[(β+1) * θ],{θ, k}]*Log[r]^(j-k),{k, 0, j}],{j, 0, n}]; u[[k1, l +1+1]] = Simplify[u[[k1, l +1+1]]] ]; Print["u_kl=" u[[k1, l +1+1]]]; If[Dimensions[H][[2]] ⩵0, u[[k1, l +1+1]] = 0]; Print["Lower Iteration complete"]; ,{l, 0, m -2-h[[k1]]}]; Print["Higher Iteration Complete"] Print["Dimensions=" Dimensions[u][[2]]] Φ[[k1]] = Sum[u[[k1, j1]],{j1, 1, Dimensions[u][[2]]}]; ,{k1, 1, Dimensions[λ][[1]]}]; ] Do[Φ[[i]] = Sum[u[[i, j]],{j, 1, Dimensions[u][[2]]}], {i, 1, Dimensions[Φ][[1]]}] If[λ[[Dimensions[λ][[1]]]] ⩵m-1, Φ[[Dimensions[λ][[1]]]] = 0]; ω=ωaux; Print["ω"] (*This piece of code runs the modified version of the method. Otherwise, it is identical, and accepts exactly the same datafile.*) ClearAll[m, θi, k1, k2, bbcc, ω, nterms, λ, j, a1, u0, u1, r, θ, u, a, h, A, G, B, Φ,Ω,β,γ1, γ2, Γ] (*Reading block, in this section of the code we should be able to extract the data for use: the regularity of the solution m, a vector with the corner's angles θi(two components), three constants with the material's properties and another vector containing the boundary conditions, bbcc.*) SetDirectory["C:\Cosas de Ingenieros\TFG\Ejecutar"]; infile =InputString["Input file?"]; input =ReadList[infile, Number, RecordLists →True]; m=input[[1, 1]]; θi=input[[2]] * π/180; (*γ1=input[[3,1]];γ2=input[[3,2]];*) bbcc =input[[4]]; ω=θi[[1]] - θi[[2]]; ωaux =ω; If[bbcc[[1]] ⩵bbcc[[2]], nterms =IntegerPart[ω/π* (m-1)], nterms =IntegerPart[ω/π* (m-1)+1/2] ] (*θi is called like that to avoid confusion with the polar angle*) (*Definition of the problem module, here we need to identify the kind of problem we have introduced with the datafile *) (*WARNING: THESE EQUATIONS ARE WRITTEN UNDER MGHAZLI'S POLAR COORDINATES, WITH ORIGIN ON ONE OF THE BOUNDARIES OF THE WEDGES.*) (*Table of groups: M: Mixed boundary conditions. 6
48 Apéndice A. Códigos de Mathematica S: Same boundary conditions R:Robin condition D: Dirichlett condition N: Neumann condition *) (*The philosophy now is as follows: Once the datafile has been read and the bbcc vector has been created, it is compared in succession with different sets of boundary conditions. If the bbcc vector doesn't fulfill the conditions, it is compared with the next until the correct one is found. *) j; If[bbcc[[1]] ⩵bbcc[[2]], λ=Table[j*Pi /ω,{j, nterms}], λ=Table[(j-1/2)*Pi /ω,{j, nterms}]] h=Table[Ceiling[λ[[j]]],{j, nterms}] If[bbcc[[1]] ⩵2 && bbcc[[2]] == 0, (*For this situation, Mghazli's notation would give j contained in M,R*) (*As for "our" indexes, j runs from 1 to the dimension of λ, while l runs from 0 to m-h-2, in order to give shadow terms up to m-h-1. Thus, l needs an offset in order to be put in Mathematica's lists, while j doesn't.*) Φ=Table[0, {j, Dimensions[λ][[1]]}]; u=Table[0, {i, 1, Dimensions[λ][[1]]},{k, 1, m}]; Do[ If[! Element[λ[[j]], Integers], u[[j, 1]]=r^(λ[[j]]) * Sin[λ[[j]]*θ], u[[j, 1]] = r^(λ[[j]]) * Sin[λ[[j]] * θ]] ,{j, 1, Dimensions[λ][[1]]}] Do[ Do[ (*Within this loop, we are supposed to find each one of the singular functions or "shadow terms" associated with the problem for a "frozen" λ_k. The process is as follows: Each PDE results in the shadow term u_l, and is in turn given as a particularization for θ=ωof the prior shadow term u_(l-1)The parameters we need to define for the problem are β, n, and G. The G coefficients are those that multiply each one of the Log[r]terms in the prior shadow term.*) (*The general expression for the Neumann nonhomogeneous boundary condition is - r*γ*u_(l-1). Obviously each shadow term can be regarded as a polynomial in Log[r]multiplied by a power of r, the power being r^(λ_k+l-1).*) (*Thus, it is trivial to see that the shadow term u_l has a power superior to the one in u_(l-1)by one.*) (*A very important remark is that the Log[r]terms won't even exist in the first place if λ_k is not an integer, since the shadow terms are derived in succession, if the first shadow term doesn't have a logarithm, those that follow will not present logarithmic behaviour.*) (*Remark: The proof of the global theorem is valid for q=0, 1,2...,m-h-2. I don't have to worry about the upper limit, since q=m-h-2 would give a PDE for which the solution is the shadow term of order m-h-1. *) β=λ[[k1]] + l+1-1; (*This expression matches the exponent of 7
49 r in the complete problem with that of Mghazli's appendix.*) θ=ω; Ω= -r*γ1*u[[k1, l +1]]; ClearAll[θ]; (*We know the power of the u_l shadow term is λ_k+l, so if we divide Ωby r^(λ_k+l+1)we have a pure polynomial in Log[r] (if there are logarithms at all)*) G=CoefficientList[Ω/r^(β+1), Log[r]]; n=Dimensions[G, 1][[1]]; Print["G=", G]; Print["DimensionG=", Dimensions[G]]; Print["n=", n]; Print["β=", β]; Print["ω=", ω]; Print["Ω=", Ω]; A=Table[0, {i, 1, n +1}]; Ct = Table[If[k⩵1, 1, Product[(i-1+l),{l, 1, k -1}] / (k-1)! / (β+1)^(k-1)], {i, 1, n +1},{k, 1, n +1}]; Print[Ct]; Print[ω]; ClearAll[ω]; If[Cos[(β+1)*ωaux]≠0 && Dimensions[G][[1]] ≠0, A[[n+1]]=0; A[[n+1-1]] = G[[n+1-1]] / (β+1)/Cos[(β+1)ω]; Do[ A[[i+1]]=1/ (β+1)/Cos[(β+1)ω] * (G[[i+1]] - Sum[A[[i+1+k]] * Ct[[i+1, k +1]] * ω^(k-1)*(k*D[Sin[(β+1) * ω], {ω, k}] + ω*D[Sin[(β+1) * ω],{ω, k +1}]),{k, 1, n -i}]); ,{i, n -2, 0, -1}]; ω=ωaux; u[[k1, l +1+1]] = r^(β+1)*(Sum[A[[j+1]] * Sum[Ct[[j+1-k, k +1]]*θ^k * D[Sin[(β+1) * θ] * Log[r]^(j-k),{θ, k}],{k, 0, j}],{j, 0, n}]); ]; If[Cos[(β+1)*ωaux] == 0 && Dimensions[G][[1]] ≠0, A[[1]]=0; A[[n+1]] = -G[[n+1-1]]/ω/ (β+1)/n/Sin[(β+1) * ω]; Do[A[[i+1]] = 1/ω/i/ (β+1) / Sin[(β+1)*ω] * (-G[[i+1-1]] + Sum[A[[i+1-1+k]] * Ct[[i+1-1, k +1]] * ω^(k-1) * (k*D[Sin[(β+1) * ω],{ω, k}]+ω*D[Sin[(β+1) * ω],{ω, k +1}]), {k, 2, n -i+1}]),{i, n -1, 1, -1}]; ω=ωaux; u[[k1, l +1+1]] = r^(β+1)*(Sum[A[[j+1]] * Sum[Ct[[j+1-k, k +1]]*θ^k * D[Sin[(β+1) * θ] * Log[r]^(j-k),{θ, k}],{k, 0, j}],{j, 0, n}]); ]; If[Dimensions[G][[1]] ⩵0, u[[k1, l +1+1]] = 0]; Print["ωdeb=", ω]; Print["A=", A]; Print["u_kl=", u[[k1, l +1+1]]]; Print["k l=", k l] 8